Introduce image volumes as an input for feature query

CITATION.cff

@@ -15,9 +15,9 @@ license: Apache License, Version 2.0
 
 repository-code: https://github.com/FZJ-INM1-BDA/siibra-python
 
-version: v0.5a02
+version: v1.0a01
 
-date-released: 2023-11-15
+date-released: 2023-11-16
 
 doi: https://zenodo.org/records/7885729
 

docs/developer.rst

@@ -157,9 +157,14 @@ Basic definitions and notes
 Fetching volumes
 ----------------
 
+To access the actual image or mesh data from a volume requires explicit fetching.
+The process of fetching typically involves finding an appropriate volume provider, loading the data, and returning the image or mesh data.
+In simple volume objects however, this might only require to return already loaded image data. 
+Images are typially returned as Nifti1Image, while meshes are returned as a dictionary including the vertices, faces and sometimes label arrays.
+
 Fetching volumes occurs in two main stages:
 
-1. The determination of the volume by the user.
+1. Request of a volume object, such as a template, parcellation map, or other plain image volume.
 
     * The user sets the object they would like to fetch a volume from:
 
@@ -171,7 +176,7 @@ Fetching volumes occurs in two main stages:
         * template directly accesses to ``volume.fetch()``
         * ``fetch()`` first goes through ``map.fetch()`` to determine the associated volume.
 
-2. Actual retrieval of the volume object by siibra after the user asks for the
+2. Actual retrieval of the volume data by siibra after the user asks for the
    volume via ``fetch()`` method. When ``fetch()`` is invoked it accesses to
    corresponding volume provider based on the specifications given by volume
    index, fragment, z, label, variant, and format. According to the source type,

e2e/core/test_region.py

@@ -4,7 +4,8 @@
 from itertools import repeat
 
 import siibra
-from siibra.core.relation_qualification import Qualification
+from siibra.core.assignment import Qualification
+from siibra.core.region import Region
 
 
 regions = [
@@ -14,7 +15,7 @@
 
 
 @pytest.mark.parametrize("region", regions)
-def test_region_spatial_props(region: siibra.core.parcellation.region.Region):
+def test_region_spatial_props(region: Region):
     props = region.spatial_props("mni152")
     for idx, cmp in enumerate(props.components, start=1):
         assert cmp.volume >= props.components[idx - 1].volume
@@ -76,7 +77,7 @@ def test_find(parc_spec, region_spec, result_len, check_regions):
 
 
 @pytest.mark.parametrize("parc, reg_spec, has_related, has_homology, has_related_ebrains_reg", [
-    ("2.9", "PGa", True, True, False),
+    ("julich 2.9", "PGa", True, True, False),
     ("monkey", "PG", False, True, False),
     ("waxholm v3", "cornu ammonis 1", True, False, True),
 ])
@@ -99,6 +100,7 @@ def test_homologies_related_regions(parc, reg_spec, has_related, has_homology, h
             )
         assert len([f for f in features]) > 0
 
+
 def test_related_region_hemisphere():
     reg = siibra.get_region("2.9", "PGa")
     all_related_reg = [reg for reg in reg.get_related_regions()]

e2e/features/external/test_ebrains.py

@@ -1,6 +1,6 @@
 import siibra
 import pytest
-from siibra.features.anchor import AssignmentQualification
+from siibra.features.anchor import Qualification
 
 concepts = [
     (
@@ -25,4 +25,4 @@ def features_fixture():
 def test_each_feature(features_fixture):
     for f in features_fixture:
         for qual in f.anchor.regions.values():
-            assert isinstance(qual, AssignmentQualification)
+            assert isinstance(qual, Qualification)

e2e/features/image/test_image.py

@@ -2,7 +2,7 @@
 import pytest
 
 features = siibra.features.get(
-    siibra.spaces["big brain"], "CellBodyStainedVolumeOfInterest"
+    siibra.get_template("big brain"), "CellBodyStainedVolumeOfInterest"
 )
 
 

e2e/features/test_get.py

@@ -11,7 +11,7 @@ def test_get_instances(Cls: siibra.features.Feature):
     assert isinstance(instances, list)
 
 
-ids = [
+selected_ids = [
     "lq0::EbrainsDataFeature::p:minds/core/parcellationatlas/v1.0.0/94c1125b-b87e-45e4-901c-00daee7f2579-290::r:Area hOc1 (V1, 17, CalcS) left::https://nexus.humanbrainproject.org/v0/data/minds/core/dataset/v1.0.0/3ff328fa-f48f-474b-bd81-b5ee7ca230b6",
     "cf0::BigBrainIntensityProfile::p:minds/core/parcellationatlas/v1.0.0/94c1125b-b87e-45e4-901c-00daee7f2579-290::r:Area hOc1 (V1, 17, CalcS) left::nodsid::6ba1f5c180a63705d84c25db9ff7efa7",  # CompoundFeature of 1579 BigBrainIntensityProfile features grouped by (Modified silver staining modality) anchored at Area hOc1 (V1, 17, CalcS) left with Set of 1579 points in the Bounding box from (-63.69,-59.94,-29.09) mm to (0.91,77.90,54.03)mm in BigBrain microscopic template (histology) space
     "lq0::BigBrainIntensityProfile::p:minds/core/parcellationatlas/v1.0.0/94c1125b-b87e-45e4-901c-00daee7f2579-290::r:Area hOc1 (V1, 17, CalcS) left::a48879314d93d15c74e6512e9af9c91d",  # BigBrainIntensityProfile (Modified silver staining) anchored at Area hOc1 (V1, 17, CalcS) left with Point in BigBrain microscopic template (histology) [0.4248340129852295,50.589298248291016,-14.839699745178223]
@@ -27,7 +27,7 @@ def test_get_instances(Cls: siibra.features.Feature):
 ]
 
 
-@pytest.mark.parametrize("fid", ids)
+@pytest.mark.parametrize("fid", selected_ids)
 def test_get_instance(fid):
     feat = siibra.features.Feature._get_instance_by_id(fid)
     assert feat.id == fid
@@ -36,7 +36,7 @@ def test_get_instance(fid):
 # this tests whether or not calling a live query caused proxy feature to be
 # added to subclasses. (It should not: causes memory leak and also increases
 # query time linearly)
-@pytest.mark.parametrize("fid", ids)
+@pytest.mark.parametrize("fid", selected_ids)
 def test_subclass_count(fid):
     len_before = len(siibra.features.Feature.SUBCLASSES[siibra.features.Feature])
     _ = siibra.features.Feature._get_instance_by_id(fid)

e2e/volumes/test_mapshape.py

@@ -39,7 +39,7 @@ def test_map_shape(space_shape, parc_id, map_endpoint, region_name):
     volume_data = None
     if region_name is not None:
         region = siibra.get_region(parc_id, region_name)
-        volume_data = region.fetch_regional_map(space_id, map_endpoint)
+        volume_data = region.get_regional_map(space_id, map_endpoint).fetch()
     else:
         labelled_map = siibra.get_map(parc_id, space_id, map_endpoint)
         assert labelled_map is not None
@@ -56,5 +56,5 @@ def test_template_resampling():
     mp_img = mp.fetch()
     template = mp.space.get_template().fetch()
     assert mp_img.shape != template.shape
-    resamp_template = mp.get_resampled_template()
-    assert mp_img.shape == resamp_template.shape
+    resamp_template = mp.get_resampled_template().fetch()
+    assert mp_img.shape == resamp_template.shape

e2e/volumes/test_parcellationmap.py

@@ -103,3 +103,10 @@ def test_region_1to1ness_in_parcellation():
                     }
                 )
     assert len(failed) == 0, print("The regions in maps that can't be gotten in respective parcellations\n", failed)
+
+
+def test_fetching_merged_volume():
+    # this also tests for volume.merge as it is used in Map.fetch
+    mp = siibra.get_map("julich 2.9", "bigbrain")
+    assert len(mp) > 1
+    _ = mp.fetch()

e2e/volumes/test_volume.py

@@ -0,0 +1,35 @@
+import siibra
+from siibra.volumes import Volume
+from siibra.core.region import Region
+import numpy as np
+import pytest
+# TODO write a test for the volume-region and volume-volume intersection
+
+
+# add more to the list when centroids can be calculated for non-convex regions as well
+selected_regions = [
+    (siibra.get_region('julich 2.9', 'CA2 (Hippocampus) right'), 'mni152'),
+    (siibra.get_region('julich 2.9', 'CA2 (Hippocampus) left'), 'colin27'),
+    (siibra.get_region('julich 2.9', 'hoc1 left'), 'mni152'),
+    (siibra.get_region('julich 3', 'sts'), 'mni152')
+]
+
+
+@pytest.mark.parametrize("region, space", selected_regions)
+def test_region_intersection_with_its_own_volume(region, space):
+    assert isinstance(region, Region)
+    volume = region.get_regional_map(space)
+    intersection = region.intersection(volume)
+    assert isinstance(intersection, Volume)
+    assert np.all(
+        np.equal(intersection.fetch().dataobj, volume.fetch().dataobj)
+    ), "Intersection of a regional map with its region object should be the same volume."
+
+
+@pytest.mark.parametrize("region, space", selected_regions)
+def test_region_intersection_with_its_centroid(region, space):
+    assert isinstance(region, Region)
+    centroids = region.compute_centroids(space)
+    assert centroids in region
+    assert region.intersection(centroids) == centroids
+    assert centroids.intersection(region) == centroids

examples/01_atlases_and_parcellations/005_brain_region_spatialprops.py

@@ -61,6 +61,6 @@
 # gives us a Nifti1Image object as provided by `nibabel <https://nipy.org/nibabel/>`_, 
 # and which we can directly visualize using plotting functions like the ones in 
 # `nilearn <https://nilearn.github.io/stable/index.html>`_:
-mask = v1_left.fetch_regional_map(space, maptype="labelled")
+mask = v1_left.get_regional_map(space, maptype="labelled")
 from nilearn import plotting
-plotting.plot_roi(mask, title=f"Mask of {v1_left.name} in {space.name}")
+plotting.plot_roi(mask.fetch(), title=f"Mask of {v1_left.name} in {space.name}")

examples/02_maps_and_templates/003_accessing_maps.py

@@ -42,9 +42,9 @@
 
 # %%
 # We select the maximum probability map of Julich-Brain in MNI152 space,
-# which is a parcellation map with discrete labels. 
+# which is a parcellation map with discrete labels.
 # `get_map` assumes maptype='labelled' by default.
-julich_mpm = siibra.get_map(space="icbm 2009c", parcellation="julich 2.9", maptype="labelled")
+julich_mpm = siibra.get_map(space="icbm 2009c", parcellation="julich 3", maptype="labelled")
 print(julich_mpm)
 
 # %%
@@ -120,6 +120,5 @@
 # %%
 # In addition to parcellation maps, `siibra` can produce binary masks of brain regions.
 hoc5L = siibra.get_region(parcellation='julich 2.9', region='hoc5 left')
-hoc5L_mask = hoc5L.fetch_regional_map(space="mni152", maptype="labelled")
-plotting.plot_roi(hoc5L_mask, title=f"Mask of {hoc5L.name}")
-
+hoc5L_mask = hoc5L.get_regional_map(space="mni152", maptype="labelled")
+plotting.plot_roi(hoc5L_mask.fetch(), title=f"Mask of {hoc5L.name}")

examples/02_maps_and_templates/004_access_bigbrain.py

@@ -35,9 +35,10 @@
 # %%
 # Per default, `siibra` will fetch the whole brain volume at a reasonably
 # reduced resolution.
-bigbrain = siibra.get_template('bigbrain')
-bigbrain_whole = bigbrain.fetch()
-plotting.view_img(bigbrain_whole, bg_img=None, cmap='gray')
+space = siibra.spaces['bigbrain']
+bigbrain_template = space.get_template()
+bigbrain_whole_img = bigbrain_template.fetch()
+plotting.view_img(bigbrain_whole_img, bg_img=None, cmap='gray')
 
 # %%
 # To see the full resolution, we may specify a bounding box in the physical
@@ -47,19 +48,24 @@
 # a string representation, which could be conveniently copy pasted from the
 # interactive viewer `siibra explorer <https://atlases.ebrains.eu/viewer>`_.
 # Note that the coordinates can be specified by 3-tuples, and in other ways.
-space = siibra.spaces.get('bigbrain')
-voi = space.get_bounding_box(
+voi = siibra.locations.BoundingBox(
     point1="-30.590mm, 3.270mm, 47.814mm",
-    point2="-26.557mm, 6.277mm, 50.631mm"
+    point2="-26.557mm, 6.277mm, 50.631mm",
+    space=space
 )
-bigbrain_chunk = bigbrain.fetch(voi=voi, resolution_mm=0.02)
+bigbrain_chunk = bigbrain_template.fetch(voi=voi, resolution_mm=0.02)
 plotting.view_img(bigbrain_chunk, bg_img=None, cmap='gray')
 
 # %%
 # Note that since both fetched image volumes are spatial images with a properly
 # defined transformation between their voxel and physical spaces, we can
 # directly plot them correctly superimposed on each other:
-plotting.view_img(bigbrain_chunk, bg_img=bigbrain_whole, cmap='magma', cut_coords=tuple(voi.center))
+plotting.view_img(
+    bigbrain_chunk,
+    bg_img=bigbrain_whole_img,
+    cmap='magma',
+    cut_coords=tuple(voi.center)
+)
 
 # %%
 # Next we select a parcellation which provides a map for BigBrain, and extract
@@ -103,7 +109,7 @@
 section1402 = features[3]
 plotting.plot_img(
     section1402.fetch(),
-    bg_img=bigbrain_whole,
+    bg_img=bigbrain_whole_img,
     title="#1402",
     cmap='gray'
 )

examples/02_maps_and_templates/006_bigbrain_cortical_layers.py

@@ -77,4 +77,3 @@
     surf_map=thicknesses[4],
     symmetric_cmap=False, cmap='magma', vmax=0.3
 )
-# %%

examples/02_maps_and_templates/007_adding_custom_parcellation.py

@@ -0,0 +1,83 @@
+# Copyright 2018-2023
+# Institute of Neuroscience and Medicine (INM-1), Forschungszentrum Jülich GmbH
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+"""
+Adding a custom parcellation map
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Sometimes you might want to use a custom parcellation map to perform feature
+queries in siibra. For this example, we retrieve a freely available AICHA -
+Atlas of Intrinsic Connectivity of Homotopic Areas (Tzourio-Mazoyer N, Landeau B,
+Papathanassiou D, Crivello F, Etard O, Delcroix N, Mazoyer B, Joliot M (2002)
+Automated anatomical labeling of activations in SPM using a macroscopic
+anatomical parcellation of the MNI MRI single-subject brain. Neuroimage
+15:273-289.). This atlas provided in the MNI ICBM 152 space.
+"""
+
+
+# %%
+import siibra
+from nilearn import plotting
+
+# %%
+# Load the custom parcellation map from the online resource. For the retrieval,
+# siibra's zipfile connector is very helpful. We can use the resulting NIfTI
+# to create a custom parcellation map inside siibra.
+
+# connect to the online zip file
+conn = siibra.retrieval.ZipfileConnector(
+    "http://www.gin.cnrs.fr/wp-content/uploads/aicha_v1.zip"
+)
+
+# the NIfTI file is easily retrieved:
+nifti = conn.get("AICHA/AICHA.nii")
+# and create a volume on space MNI152 (note that this assumes our
+# external knowledge that the map is in MNI152 space)
+volume = siibra.volumes.from_nifti(nifti, 'mni152', "AICHA")
+
+# The text file with label mappings has a custom format. We provide a tsv
+# decoder to extract the list of region/label pairs since the txt file is tab
+# seperated.
+volume_info = conn.get("AICHA/AICHA_vol1.txt", decode_func=siibra.retrieval.requests.DECODERS['.tsv'])
+volume_info
+
+# %%
+# Now we use this to add a custom map to siibra.
+regionnames = [
+    name.replace('-R', ' right').replace('-L', ' left')
+    for name in volume_info['nom_l']
+]
+labels = [int(label) for label in volume_info['color']]
+custom_map = siibra.volumes.parcellationmap.from_volume(
+    name="AICHA - Atlas of Intrinsic Connectivity of Homotopic Areas",
+    volume=volume,
+    regionnames=regionnames,
+    regionlabels=labels
+)
+
+# %%
+# let's plot the final map
+plotting.plot_roi(custom_map.fetch())
+
+
+# %%
+# We can already use this map to find spatial features, such as BigBrain
+# intensity profiles.
+region = custom_map.parcellation.get_region('S_Rolando-1 left')
+profiles = siibra.features.get(
+    region,
+    siibra.features.cellular.BigBrainIntensityProfile
+)[0]
+print(f"{len(profiles)} intensity profiles found.")

examples/03_data_features/002_colorize_map.py

@@ -23,10 +23,11 @@
 """
 
 # %%
-# We start by selecting an atlas and the Julich-Brain parcellation.
+from nilearn import plotting
 import siibra
+# We start by selecting an atlas and the Julich-Brain parcellation.
 atlas = siibra.atlases.MULTILEVEL_HUMAN_ATLAS
-jubrain = atlas.get_parcellation('julich 2.9')
+jubrain = atlas.get_parcellation('julich 3')
 # sphinx_gallery_thumbnail_path = '_static/example_thumbnails/03-002.png'
 
 # %%
@@ -53,8 +54,7 @@
 
 # %%
 # Now colorize the Julich-Brain maximum probability map and plot it.
-colorized_map = jubrain.get_map(space='mni152').colorize(mapping)
-from nilearn import plotting
+colorized_map = jubrain.get_map(space='mni152').colorize(mapping).fetch()
 plotting.view_img(
     colorized_map, cmap='magma',
     title=f"Average densities available for {receptor}", symmetric_cmap=False
@@ -64,6 +64,3 @@
 # Alternatively, we can display this map on a surface mesh using nilearn.
 # Note that, you can switch between the hemispheres or variants (inflated or pial) from the plot itself.
 plotting.view_img_on_surf(colorized_map, cmap='magma', symmetric_cmap=False, surf_mesh="fsaverage6")
-
-
-# %%