Introduce image volumes as an input for feature query (#480)

* remove slicing operator for building volumes of interest from Space objects

* remove setting __repr__ to name in Region

* introduce NIfTI FeatureMaps as Location subclass and move assignment logic from AnatomicalAnchor to Region/Location classes

* deal more consistently with points that could not be warped

* still accept Map.get_index() for multiple matches in case there is one exact match.

* introduce spatialmap as a location object; simplify matching in live queries

* WIP: introduce LocationFilter, make Volume a LocationFilter; use LocationFilters for feature queries; move providers to subfolder of Volume; rename fetch_regional_map to get_regional_map; return Nifti1Image objects only on fetch

* first attempt to re-enable feature query by region object - still testing and needs profiling

* introducde SpaceWapringFailedException; remove contains implementation in Region

* fix sparsemap.fetch, small bug fixes

* _assign_image deal with updated commons.split_components api

* Region.spatial_props accept modified get_regional_map function

* fix coordinate selection in Volume.points_inside

critical error: used voxels in wapred space instead original coordinates.

* fix collecting coordinates for feature anchor

* fix syntax error

* fix typo

* avoid redundant calls of child region maps in Region.get_regional_map

* raise proper exception in volume.warp

* Fix recursion error due to locationfilter.assign. Replace __contains__ with contains

* overwrite LocationFilter.assign() in Region to handle all cases

* fix order of "in" operator

* remove debug output

* allow quering receptor with a nifti volume

* accept volumes for intersection computations

* rename LocationFilter to BrainStructure

* some description what a brain structure is as opposed to a brain region

* fix issues caused by contains->__contains__

* Use tuple instead of subscripted generics in isinstance

As of python 3.7 Subscripted generics cannot be used with class and instance checks

* Fix import errors in tests

* allow sampling and clustering from volumes

- new method draw_samples() for (volumetric) volumes
- new methods find_clusters and label_colors for PointSet
- new method siibra.volumes.from_pointset

* add hdbscan to requirements, make label colors a property of pointset

* use sigma_mm of pointset for volume.from_points

* If volume intersection results in zeroslike, return None

* add find_peaks method to all volumes

* Fix image feature query (case: inverse assignment is None)

* fix from_pointset function in volume

* a few fixes

* a few fixes

* Fix small bugs and tests

* add matplotlib to test requirements

* More fixes

* Fix some bugs, examples, e2e tests, and all unit tests. Add hdbscan to setup requirements

* use scikit-learn HDBSCAN instead of hdbscan

* map.assign takes locations. _assign_image->_assign_volume

* Bug fix: connected_components has no element == 0

* Fix fetch cache for volumes

* Fix assignment of location to regions

If a region does not have a valid mask on a space, check if warp the locations to supported spaces of regional mask until the assignment qualification can be checked.

* Array conversion with `item` (Deprecation warning NumPy 1.25)

* Check child region masks for intersection if affine of the children does not match

* Use neuroglancer format if not format is given but resolution_mm is given

- Resolution selection is not implemented for nifti.
- This is important because there is a bug that alters which format is selected cyclically. (Could not pinpoint the source yet)

* Todo artifact and typo fix

* fix fetch format selection, format order. add _maintain_fetch_cache

* add notes to volume.fetch and protect against poisoning _FETCH_CACHE

* Array conversion with `item` (Deprecation warning NumPy 1.25)

* maint: parcmap

BUG report: _assign_volume does not work

* Maint: anchor

* sparsemap._assigne_image -> _assign_volume

* Clean up map._assing_volume. fetch_iter goes through fragments as well

* fix volume.points_inside

* feat: add maps from nifti and label information

This is meant to be a preparation for querying features with volumes. The previous iterations required computation of the relations to the defined maps. With ability top query with volumes, this computation is unnecessary.

In addition, this PR allows creating a map from multiple volumes and supplying parcellation specification.

TODO:
- [] a method for adding custom parcellations.
- [] a method for adding custom space.
- [] a method to extracting these custom atlas concepts as json preconfiguration

* Use AICHA atlas as an example

* Anchor: use a species specific match cache key

Solves the bug of fetching ebrains dataset features wrongly

* Use volume for querying image in 03.000 example

* carry from_volume to parcellationmap submodule

* Remove strict requirement of scikit-learn>=1.3.0 for finding clusters in pointsets

* remove hdbscan from requirements

* remove volume._maintain_fetch_cache

* hide volum.points_inside

* fix several bugs from the merge with 0.5a01

* Add unit to tqdm for building_objects

* Allow merging volumes from volume submodule

* Fix: Region.intersection(Location) if region map is not available on Locations.space but child regions are.

* Fix RegionRelationAssessments

* maint: location.union

* add some documentation to boundingbox

* Fix assignment.Qualification.inverse

* Fix: Volume.intersection(PointSet)

BrainStructure.intersects check for not None so return None if intersection has no points, not an empty PointSet

* Test: fetching a map with multiple volumes wo args and volume.merge

* Maint tests, fix volume.intersection, organize TODOs

* get_uuid -> generate_uuid

* Add tests: check region-volume intersection, centroid containednedness

* Make BigBrainIntensityProfile.location a property and refer to its anchor

* Remove gifit import from volumes package

* Remove redundant tests from e2e\feature\test_get

* Make Location a BrainStructure. Use this to ensure implementation of hash, eq, and repr

* update str and repr of core concepts

* Several bug fixes

-  Fix fetch cache
- fix region.fetch_cache -> getmap_cache
- Clean up get_regional_map
- Add names to volumes when extacting volumes (important for hashing)

* no name field bug fix volume.from_nifti

* further volume name adjustments for get_regional map

* Fix PointSet.intersection (rewrite)

also fix bounding box import

* Fix Point le and ge. Use other's intersection in intersection when not clear (this is established in all others)

* Fix boundingbox.intersection

* Remove non-convex region maps from volume tests

* fix try-except block of region assignment to locations

* Update compound feature _export

* Remove redundant `bigbrain._choose_space` and tidy bigbrain.py up

* add/correct type hints in anchor.py

* Use location's space for region assignment first. region.supported_spaces is now a sorted list by default

Fixes issues with reproducibility.

* Example 03-000: use region object for PLI feature query

* maint: allen query

* fix test_voi and enable more unit tests

* version bump v1.0a01

---------

Co-authored-by: Timo Dickscheid <[email protected]>

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")
-
-
-# %%