Merge pull request #523 from FZJ-INM1-BDA/docs_add_location_examples

Add the first location example

examples/04_locations/000_employing_locations_of_interest.py

@@ -0,0 +1,101 @@
+# 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.
+
+"""
+Utilizing locations of interest
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+`siibra` provides common locations of interests on reference spaces as objects.
+These location objects can be used to query features or assignment to the maps,
+see :ref:`sphx_glr_examples_05_anatomical_assignment_001_coordinates.py` and
+:ref:`sphx_glr_examples_05_anatomical_assignment_002_activation_maps.py`.
+The conversion of these locations to other spaces are done in the background
+but one can also be invoked when needed.
+"""
+
+# %%
+import siibra
+from nilearn import plotting
+import numpy as np
+
+# %%
+# The simplest location type is a point. It can have a uncertainity variable
+# or can be exact.
+point = siibra.Point((27.75, -32.0, 63.725), space='mni152')
+point_uncertain = siibra.Point((27.75, -32.0, 63.725), space='mni152', sigma_mm=3.)
+point_uncertain
+
+# %%
+# From several points or a set of coordinates we can create a PointSet.
+siibra.PointSet(
+    [(27.75, -32.0, 63.725), (27.75, -32.0, 63.725)],
+    space='mni152',
+    sigma_mm=[0.0, 3.]
+)
+
+# %%
+# There are several helper proeprties and methods for locations, some specific
+# to the location type. For example, we can warp the points to another space
+# (currently limited), to COLIN27 or BigBrain space
+print(point.warp('bigbrain'))
+print(point.warp('colin27'))
+
+# %%
+# To explore further, let us first create a random pointset and get the box
+# that contains these points, which is another location type.
+ptset = siibra.PointSet(
+    np.concatenate([
+        np.random.randn(1000, 3) * 5 + (-27.75, -32.0, 63.725),
+        np.random.randn(1000, 3) * 5 + (27.75, -32.0, 63.725)
+    ]),
+    space='mni152'
+)
+ptset.boundingbox
+
+# %%
+# We can display these points as a kernel density estimated volume
+ptset.labels = np.ones(len(ptset), dtype=int)
+kde_volume = siibra.volumes.from_pointset(ptset)
+plotting.view_img(kde_volume.fetch())
+
+# %%
+# `siibra` can find the clusters (using HDBSCAN) and label the points.
+ptset.find_clusters()
+ptset.labels += (1 - ptset.labels.min())  # offset the labels to be able to display as nifti
+clusters_kde_volume = siibra.volumes.from_pointset(ptset)
+plotting.view_img(clusters_kde_volume.fetch())
+
+# %%
+# Moreover, a location objects can be used to query features. For illusration,
+# we first crate a BoundingBox
+bbox = siibra.locations.BoundingBox(
+    point1=(-29.75, -33.0, 63.725),
+    point2=(-25.75, -30.0, 60.725),
+    space='mni152'
+)
+features_of_interest = siibra.features.get(bbox, 'image')  # let us search for images
+# and print the comparisons of the anatomical anchors to our BoundingBox
+for f in features_of_interest:
+    print(f.last_match_description)
+
+# %%
+# And now let us simply select the features that overlaps with our BoundingBox
+selected_features = [
+    f
+    for f in features_of_interest
+    if "overlaps" in str(f.last_match_result[0].qualification)
+]
+for f in selected_features:
+    print(f.name)

examples/04_locations/README.txt

@@ -3,6 +3,10 @@
 Locations in reference spaces
 -----------------------------
 
-coming soon
+``siibra`` supports defining several types of locations in reference spaces for the analysis
+of point or region of interests. Using the location objects, you can query for features
+and employ their utility functions for detailed analysis. Additionally, ``siibra`` supports
+warping coordinates between most popular spaces to enable users to study these locations of
+interests in different spaces and parcellation maps. 
 
 

siibra/volumes/volume.py

@@ -28,7 +28,8 @@
 from typing import List, Dict, Union, Set, TYPE_CHECKING
 from time import sleep
 import json
-from skimage import filters, feature as skimage_feature
+from skimage import feature as skimage_feature
+from skimage.filters import gaussian
 
 if TYPE_CHECKING:
     from ..retrieval.datasets import EbrainsDataset
@@ -457,7 +458,7 @@ def __init__(self, parent_volume: Volume, z: int):
         )
 
 
-def from_file(filename: str, space: str, name: str):
+def from_file(filename: str, space: str, name: str) -> Volume:
     """ Builds a nifti volume from a filename. """
     from ..core.concept import get_registry
     from .providers.nifti import NiftiProvider
@@ -469,7 +470,7 @@ def from_file(filename: str, space: str, name: str):
     )
 
 
-def from_nifti(nifti: Nifti1Image, space: str, name: str):
+def from_nifti(nifti: Nifti1Image, space: str, name: str) -> Volume:
     """Builds a nifti volume from a Nifti image."""
     from ..core.concept import get_registry
     from .providers.nifti import NiftiProvider
@@ -486,7 +487,7 @@ def from_array(
     affine: np.ndarray,
     space: Union[str, Dict[str, str]],
     name: str
-):
+) -> Volume:
     """Builds a siibra volume from an array and an affine matrix."""
     if len(name) == 0:
         raise ValueError("Please provide a non-empty string for `name`")
@@ -503,25 +504,63 @@ def from_array(
 
 def from_pointset(
     points: pointset.PointSet,
-    label: int,
-    target: Volume,
+    labels: List[int] = [],
+    target: Volume = None,
     min_num_points=10,
     normalize=True,
     **kwargs
-):
+) -> Volume:
+    """
+    Get the kernel density estimate as a volume from the points using their
+    average uncertainty on target volume.
+    Parameters
+    ----------
+    points: pointset.PointSet
+    target: Volume, default: None
+        If no volumes supplied, the template of the space points are defined on
+        will be used.
+    labels: List[int], default: []
+    min_num_points: int, default 10
+    normalize: bool, default: True
+
+    Raises
+    ------
+    RuntimeError
+        If no points with labels found
+    """
+    if target is None:
+        target = points.space.get_template()
     targetimg = target.fetch(**kwargs)
     voxels = points.transform(np.linalg.inv(targetimg.affine), space=None)
-    selection = [_ == label for _ in points.labels]
-    if np.count_nonzero(selection) == 0:
-        raise RuntimeError(f"No points with label {label} in the set: {', '.join(map(str, points.labels))}")
-    X, Y, Z = np.split(
-        np.array(voxels.as_list()).astype('int')[selection, :],
-        3, axis=1
-    )
-    if len(X) < min_num_points:
-        return None
     cimg = np.zeros_like(targetimg.get_fdata())
-    cimg[X, Y, Z] += 1
+
+    if len(labels) == 0:
+        if points.labels is None:
+            labels = [1]
+            logger.info("Found no point labels. Labelling all with 1.")
+            points.labels = np.ones(len(points), dtype=int)
+        else:
+            labels = points.labels
+    found_enough_points = False
+    for label in labels:
+        selection = [_ == label for _ in points.labels]
+        if selection.count(True) == 0:
+            logger.warning(f"No points with label {label} in the set: {set(points.labels)}")
+            continue
+        X, Y, Z = np.split(
+            np.array(voxels.as_list()).astype('int')[selection, :],
+            3,
+            axis=1
+        )
+        if len(X) < min_num_points:
+            logger.warning(f"Not enough points (<{min_num_points}) with label {label}, skipping.")
+            continue
+        found_enough_points = True
+        cimg[X, Y, Z] = label
+
+    if not found_enough_points:
+        raise RuntimeError(f"Not enough poinsts with labels {labels} found in {points}.")
+
     if isinstance(points.sigma_mm, (int, float)):
         bandwidth = points.sigma_mm
     elif isinstance(points.sigma_mm, list):
@@ -532,14 +571,14 @@ def from_pointset(
     else:
         logger.warn("Poinset has no uncertainty, using bandwith=1mm for kernel density estimate.")
         bandwidth = 1
-    data = filters.gaussian(cimg, bandwidth)
+    data = gaussian(cimg, bandwidth)
     if normalize:
         data /= data.sum()
     return from_array(
         data=data,
         affine=targetimg.affine,
         space=target.space,
-        name=f'KDE map of {sum(selection)} points with label={label}'
+        name=f'KDE map of {points} with labels={labels}'
     )