Merge pull request #515 from FZJ-INM1-BDA/fix_region_get_boundingbox

`provider.boundingbox` -> `provider.get_boundingbox()`

Author: AhmetNSimsek

e2e/core/test_region.py

@@ -106,3 +106,20 @@ def test_related_region_hemisphere():
     all_related_reg = [reg for reg in reg.get_related_regions()]
     assert any("left" in ass.assigned_structure.name for ass in all_related_reg)
     assert any("right" in ass.assigned_structure.name for ass in all_related_reg)
+
+
+spaces = ['mni152', 'colin27']
+
+
+@pytest.mark.parametrize("space", spaces)
+def test_boundingbox(space):
+    hoc1_l = siibra.get_region('julich', 'hoc1 left')
+    hoc1_r = siibra.get_region('julich', 'hoc1 right')
+    bbox_l = hoc1_l.get_bounding_box(space)
+    bbox_r = hoc1_r.get_bounding_box(space)
+    assert bbox_l != bbox_r, "Left and right hoc1 should not have the same bounding boxes"
+
+    v_l = hoc1_l.get_regional_map(space)
+    v_r = hoc1_r.get_regional_map(space)
+    assert bbox_l == v_l.get_boundingbox(clip=True, background=0.0), "Boundingbox of regional mask should be the same as bouding mask of the regions"
+    assert bbox_r == v_r.get_boundingbox(clip=True, background=0.0), "Boundingbox of regional mask should be the same as bouding mask of the regions"

examples/02_maps_and_templates/004_access_bigbrain.py

@@ -117,7 +117,7 @@
 # %%
 # Let's fetch a crop inside hoc5 at full resolution. We intersect the bounding
 # box of hoc5l and the section.
-hoc5_bbox = hoc5l.get_bounding_box('bigbrain').intersection(section1402.boundingbox)
+hoc5_bbox = hoc5l.get_bounding_box('bigbrain').intersection(section1402)
 print(f"Size of the bounding box: {hoc5_bbox.shape}")
 
 # this is quite large, so we shrink it

siibra/core/region.py

@@ -17,7 +17,7 @@
 from . import concept, structure, space as _space, parcellation as _parcellation
 from .assignment import Qualification, AnatomicalAssignment
 
-from ..locations import location, boundingbox, point, pointset
+from ..locations import location, point, pointset
 from ..volumes import parcellationmap, volume
 from ..commons import (
     logger,
@@ -535,7 +535,7 @@ def get_regional_map(
             # provided by 'via_space'. Now transform the affine to match the
             # desired target space.
             intermediary_result = result
-            transform = intermediary_result.boundingbox.estimate_affine(space)
+            transform = intermediary_result.get_boundingbox(clip=True, background=0.0).estimate_affine(space)
             result = volume.from_array(
                 imgdata,
                 np.dot(transform, region_img.affine),
@@ -682,6 +682,7 @@ def get_bounding_box(
         space: _space.Space,
         maptype: MapType = MapType.LABELLED,
         threshold_statistical=None,
+        **fetch_kwargs
     ):
         """Compute the bounding box of this region in the given space.
 
@@ -705,7 +706,7 @@ def get_bounding_box(
             mask = self.get_regional_map(
                 spaceobj, maptype=maptype, threshold=threshold_statistical
             )
-            return mask.boundingbox
+            return mask.get_boundingbox(clip=True, background=0.0, **fetch_kwargs)
         except (RuntimeError, ValueError):
             for other_space in self.parcellation.spaces - spaceobj:
                 try:
@@ -718,7 +719,7 @@ def get_bounding_box(
                         f"No bounding box for {self.name} defined in {spaceobj.name}, "
                         f"will warp the bounding box from {other_space.name} instead."
                     )
-                    bbox = mask.boundingbox
+                    bbox = mask.get_boundingbox(clip=True, background=0.0, **fetch_kwargs)
                     if bbox is not None:
                         return bbox.warp(spaceobj)
                 except RuntimeError:

siibra/features/image/image.py

@@ -42,7 +42,7 @@ def location(self):
         Loads the bounding box only if required, since it demands image data access.
         """
         if self._location_cached is None:
-            self._location_cached = self.volume.boundingbox
+            self._location_cached = self.volume.get_boundingbox(clip=False)  # use unclipped to preseve exisiting behaviour
         return self._location_cached
 
     @property

siibra/locations/boundingbox.py

@@ -21,10 +21,11 @@
 
 import hashlib
 import numpy as np
-from typing import TYPE_CHECKING
+from typing import TYPE_CHECKING, Union
 if TYPE_CHECKING:
     from ..core.structure import BrainStructure
     from nibabel import Nifti1Image
+    from ..core.space import Space
 
 
 class BoundingBox(location.Location):
@@ -35,7 +36,14 @@ class BoundingBox(location.Location):
     from the two corner points.
     """
 
-    def __init__(self, point1, point2, space=None, minsize: float = None, sigma_mm=None):
+    def __init__(
+        self,
+        point1,
+        point2,
+        space: Union[str, 'Space'] = None,
+        minsize: float = None,
+        sigma_mm=None
+    ):
         """
         Construct a new bounding box spanned by two 3D coordinates
         in the given reference space.
@@ -99,24 +107,6 @@ def shape(self) -> float:
     def is_planar(self) -> bool:
         return any(d == 0 for d in self.shape)
 
-    @staticmethod
-    def _determine_bounds(A, threshold=0):
-        """
-        Bounding box of nonzero values in a 3D array.
-        https://stackoverflow.com/questions/31400769/bounding-box-of-numpy-array
-        """
-        x = np.any(A > threshold, axis=(1, 2))
-        y = np.any(A > threshold, axis=(0, 2))
-        z = np.any(A > threshold, axis=(0, 1))
-        nzx, nzy, nzz = [np.where(v) for v in (x, y, z)]
-        if any(len(nz[0]) == 0 for nz in [nzx, nzy, nzz]):
-            # empty array
-            return None
-        xmin, xmax = nzx[0][[0, -1]]
-        ymin, ymax = nzy[0][[0, -1]]
-        zmin, zmax = nzz[0][[0, -1]]
-        return np.array([[xmin, xmax + 1], [ymin, ymax + 1], [zmin, zmax + 1], [1, 1]])
-
     def __str__(self):
         if self.space is None:
             return (
@@ -390,3 +380,35 @@ def __eq__(self, other: 'BoundingBox'):
 
     def __hash__(self):
         return super().__hash__()
+
+
+def _determine_bounds(array: np.ndarray, threshold=0):
+    """
+    Bounding box of nonzero values in a 3D array.
+    https://stackoverflow.com/questions/31400769/bounding-box-of-numpy-array
+    """
+    x = np.any(array > threshold, axis=(1, 2))
+    y = np.any(array > threshold, axis=(0, 2))
+    z = np.any(array > threshold, axis=(0, 1))
+    nzx, nzy, nzz = [np.where(v) for v in (x, y, z)]
+    if any(len(nz[0]) == 0 for nz in [nzx, nzy, nzz]):
+        # empty array
+        return None
+    xmin, xmax = nzx[0][[0, -1]]
+    ymin, ymax = nzy[0][[0, -1]]
+    zmin, zmax = nzz[0][[0, -1]]
+    return np.array([[xmin, xmax + 1], [ymin, ymax + 1], [zmin, zmax + 1], [1, 1]])
+
+
+def from_array(array: np.ndarray, threshold=0, space: "Space" = None) -> "BoundingBox":
+    """
+    Find the bounding box of an array.
+
+    Parameters
+    ----------
+    array : np.ndarray
+    threshold : int, default: 0
+    space : Space, default: None
+    """
+    bounds = _determine_bounds(array, threshold)
+    return BoundingBox(bounds[:3, 0], bounds[:3, 1], space=space)

siibra/volumes/providers/gifti.py

@@ -44,8 +44,10 @@ def __init__(self, url: Union[str, Dict[str, str]], volume=None):
     def _url(self) -> Union[str, Dict[str, str]]:
         return self._init_url
 
-    @property
-    def boundingbox(self) -> _boundingbox.BoundingBox:
+    def get_boundingbox(self, clip=False, background=0.0, **fetch_kwargs) -> '_boundingbox.BoundingBox':
+        """
+        Bounding box calculation is not yet implemented for meshes.
+        """
         raise NotImplementedError(
             f"Bounding box access to {self.__class__.__name__} objects not yet implemented."
         )
@@ -143,10 +145,9 @@ def fetch(self, fragment: str = None, **kwargs):
 
         return {"labels": np.hstack(labels)}
 
-    @property
-    def boundingbox(self) -> _boundingbox.BoundingBox:
+    def get_boundingbox(self, clip=False, background=0.0) -> '_boundingbox.BoundingBox':
         raise NotImplementedError(
-            f"Bounding boxes of {self.__class__.__name__} objects not defined."
+            f"Bounding box access to {self.__class__.__name__} objects not yet implemented."
         )
 
     @property

siibra/volumes/providers/neuroglancer.py

@@ -118,22 +118,47 @@ def fetch(
 
         return result
 
-    @property
-    def boundingbox(self):
+    def get_boundingbox(self, clip=False, background=0, **fetch_kwargs) -> "_boundingbox.BoundingBox":
         """
-        Return the bounding box in physical coordinates
-        of the union of fragments in this neuroglancer volume.
+        Return the bounding box in physical coordinates of the union of
+        fragments in this neuroglancer volume.
+
+        Parameters
+        ----------
+        clip: bool, default: True
+            Whether to clip the background of the volume.
+        background: float, default: 0.0
+            The background value to clip.
+            Note
+            ----
+            To use it, clip must be True.
+        fetch_kwargs:
+            key word arguments that are used for fetchin volumes,
+            such as voi or resolution_mm.
         """
         bbox = None
         for frag in self._fragments.values():
-            next_bbox = _boundingbox.BoundingBox((0, 0, 0), frag.shape, space=None) \
-                .transform(frag.affine)
+            if len(frag.shape) > 3:
+                logger.warning(
+                    f"N-D Neuroglancer volume has shape {frag.shape}, but "
+                    f"bounding box considers only {frag.shape[:3]}"
+                )
+            if clip:
+                img = frag.fetch(**fetch_kwargs)
+                next_bbox = _boundingbox.from_array(
+                    np.asanyarray(img.dataobj), threshold=background, space=None
+                ).transform(img.affine)  # use the affine of the image matching fetch_kwargs
+            else:
+                shape = frag.shape[:3]
+                next_bbox = _boundingbox.BoundingBox(
+                    (0, 0, 0), shape, space=None
+                ).transform(frag.affine)
             bbox = next_bbox if bbox is None else bbox.union(next_bbox)
         return bbox
 
     def _merge_fragments(self) -> nib.Nifti1Image:
         # TODO this only performs nearest neighbor interpolation, optimized for float types.
-        bbox = self.boundingbox
+        bbox = self.get_boundingbox(clip=False, background=0.0)
         num_conflicts = 0
         result = None
 
@@ -486,10 +511,12 @@ def __init__(self, resource: Union[str, dict], volume=None):
     def _url(self) -> Union[str, Dict[str, str]]:
         return self._init_url
 
-    @property
-    def boundingbox(self) -> _boundingbox.BoundingBox:
+    def get_boundingbox(self, clip=False, background=0.0, **fetch_kwargs) -> '_boundingbox.BoundingBox':
+        """
+        Bounding box calculation is not yet implemented for meshes.
+        """
         raise NotImplementedError(
-            f"Fast bounding box access to {self.__class__.__name__} objects not yet implemented."
+            f"Bounding box access to {self.__class__.__name__} objects not yet implemented."
         )
 
     def _get_fragment_info(self, meshindex: int) -> Dict[str, Tuple[str, ]]:

siibra/volumes/providers/nifti.py

@@ -65,12 +65,28 @@ def _url(self) -> Union[str, Dict[str, str]]:
     def fragments(self):
         return [k for k in self._img_loaders if k is not None]
 
-    @property
-    def boundingbox(self):
+    def get_boundingbox(self, clip=True, background=0, **fetch_kwargs) -> "_boundingbox.BoundingBox":
         """
-        Return the bounding box in physical coordinates
-        of the union of fragments in this nifti volume.
+        Return the bounding box in physical coordinates of the union of
+        fragments in this nifti volume.
+
+        Parameters
+        ----------
+        clip : bool, default: True
+            Whether to clip the background of the volume.
+        background : float, default: 0.0
+            The background value to clip.
+            Note
+            ----
+            To use it, clip must be True.
+        fetch_kwargs:
+            Not used
         """
+        if fetch_kwargs:
+            logger.warning(
+                "`volume.fetch()` keyword arguments supplied. Nifti volumes"
+                " cannot pass them for bounding box calculation."
+            )
         bbox = None
         for loader in self._img_loaders.values():
             img = loader()
@@ -79,15 +95,21 @@ def boundingbox(self):
                     f"N-D NIfTI volume has shape {img.shape}, but "
                     f"bounding box considers only {img.shape[:3]}"
                 )
-            shape = img.shape[:3]
-            next_bbox = _boundingbox.BoundingBox((0, 0, 0), shape, space=None) \
-                .transform(img.affine)
+            if clip:
+                next_bbox = _boundingbox.from_array(
+                    np.asanyarray(img.dataobj), threshold=background, space=None
+                ).transform(img.affine)
+            else:
+                shape = img.shape[:3]
+                next_bbox = _boundingbox.BoundingBox(
+                    (0, 0, 0), shape, space=None
+                ).transform(img.affine)
             bbox = next_bbox if bbox is None else bbox.union(next_bbox)
         return bbox
 
     def _merge_fragments(self) -> nib.Nifti1Image:
         # TODO this only performs nearest neighbor interpolation, optimized for float types.
-        bbox = self.boundingbox
+        bbox = self.get_boundingbox(clip=False, background=0.0)
         num_conflicts = 0
         result = None
 

siibra/volumes/providers/provider.py

@@ -35,9 +35,8 @@ def __init_subclass__(cls, srctype: str) -> None:
         VolumeProvider._SUBCLASSES.append(cls)
         return super().__init_subclass__()
 
-    @property
     @abstractmethod
-    def boundingbox(self) -> "BoundingBox":
+    def get_boundingbox(self, clip=True, background=0.0) -> "BoundingBox":
         raise NotImplementedError
 
     @property
@@ -84,9 +83,8 @@ def __init__(self, parent_provider: VolumeProvider, z: int):
         self.srctype = parent_provider.srctype
         self.z = z
 
-    @property
-    def boundingbox(self) -> "BoundingBox":
-        return self.provider.boundingbox
+    def get_boundingbox(self, clip=True, background=0.0, **fetch_kwargs) -> "BoundingBox":
+        return self.provider.get_boundingbox(clip=clip, background=background, **fetch_kwargs)
 
     def fetch(self, **kwargs):
         # activate caching at the caller using "with SubvolumeProvider.UseCaching():""