Maint provider.get_boundingbox(), volume.boundingbox -> volume.get_boundingbox()

AhmetNSimsek · AhmetNSimsek · commit 72a3684f12c1 · 2023-11-21T12:48:34.000+01:00
diff --git a/e2e/core/test_region.py b/e2e/core/test_region.py
@@ -121,5 +121,5 @@ def test_boundingbox(space):
 
     v_l = hoc1_l.get_regional_map(space)
     v_r = hoc1_r.get_regional_map(space)
-    assert bbox_l == v_l.boundingbox, "Boundingbox of regional mask should be the same as bouding mask of the regions"
-    assert bbox_r == v_r.boundingbox, "Boundingbox of regional mask should be the same as bouding mask of the regions"
+    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"
diff --git a/siibra/core/region.py b/siibra/core/region.py
@@ -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:
diff --git a/siibra/features/image/image.py b/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
diff --git a/siibra/locations/boundingbox.py b/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)
diff --git a/siibra/volumes/providers/gifti.py b/siibra/volumes/providers/gifti.py
@@ -44,7 +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
 
-    def get_boundingbox(self, clip=False, background=0.0) -> '_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."
         )
diff --git a/siibra/volumes/providers/neuroglancer.py b/siibra/volumes/providers/neuroglancer.py
@@ -118,10 +118,23 @@ def fetch(
 
         return result
 
-    def get_boundingbox(self, clip=False, background=0) -> "_boundingbox.BoundingBox":
+    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():
@@ -131,11 +144,10 @@ def get_boundingbox(self, clip=False, background=0) -> "_boundingbox.BoundingBox
                     f"bounding box considers only {frag.shape[:3]}"
                 )
             if clip:
-                img = frag.fetch()
-                bounds = _boundingbox.BoundingBox._determine_bounds(np.asanyarray(img.dataobj), background)
-                next_bbox = _boundingbox.BoundingBox(
-                    bounds[:3, 0], bounds[:3, 1], space=None
-                )
+                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(
@@ -499,7 +511,10 @@ def __init__(self, resource: Union[str, dict], volume=None):
     def _url(self) -> Union[str, Dict[str, str]]:
         return self._init_url
 
-    def get_boundingbox(self, clip=False, background=0.0) -> '_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."
         )
diff --git a/siibra/volumes/providers/nifti.py b/siibra/volumes/providers/nifti.py
@@ -65,11 +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]
 
-    def get_boundingbox(self, clip=True, background=0) -> "_boundingbox.BoundingBox":
+    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,10 +96,9 @@ def get_boundingbox(self, clip=True, background=0) -> "_boundingbox.BoundingBox"
                     f"bounding box considers only {img.shape[:3]}"
                 )
             if clip:
-                bounds = _boundingbox.BoundingBox._determine_bounds(np.asanyarray(img.dataobj), background)
-                next_bbox = _boundingbox.BoundingBox(
-                    bounds[:3, 0], bounds[:3, 1], space=None
-                )
+                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(
@@ -93,7 +109,7 @@ def get_boundingbox(self, clip=True, background=0) -> "_boundingbox.BoundingBox"
 
     def _merge_fragments(self) -> nib.Nifti1Image:
         # TODO this only performs nearest neighbor interpolation, optimized for float types.
-        bbox = self.get_boundingbox(clip=True, background=0.0)
+        bbox = self.get_boundingbox(clip=False, background=0.0)
         num_conflicts = 0
         result = None
 
@@ -174,7 +190,7 @@ def fetch(
             result = loader()
 
         if voi is not None:
-            bb_vox = voi.transform_bbox(np.linalg.inv(result.affine))
+            bb_vox = voi.transform(np.linalg.inv(result.affine))
             (x0, y0, z0), (x1, y1, z1) = bb_vox.minpoint, bb_vox.maxpoint
             shift = np.identity(4)
             shift[:3, -1] = bb_vox.minpoint
diff --git a/siibra/volumes/providers/provider.py b/siibra/volumes/providers/provider.py
@@ -83,8 +83,8 @@ def __init__(self, parent_provider: VolumeProvider, z: int):
         self.srctype = parent_provider.srctype
         self.z = z
 
-    def get_boundingbox(self, clip=True, background=0.0) -> "BoundingBox":
-        return self.provider.get_boundingbox(clip=clip, background=background)
+    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():""
diff --git a/siibra/volumes/volume.py b/siibra/volumes/volume.py
@@ -115,11 +115,41 @@ def concat(url: Union[str, Dict[str, str]], concat: str):
             for srctype, prov in self._providers.items()
         }
 
-    @property
-    def boundingbox(self):
-        for provider in self._providers.values():
+    def get_boundingbox(self, clip: bool = True, background: float = 0.0, **fetch_kwargs) -> "boundingbox.BoundingBox":
+        """
+        Obtain the bounding box in physical coordinates of this 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. Currently, only possible for
+            Neuroglancer volumes except for `format`.
+
+        Raises
+        ------
+        RuntimeError
+            If the volume provider does not have a bounding box calculator.
+        """
+        fmt = fetch_kwargs.get("format")
+        if fmt in self._providers:
+            raise ValueError(
+                f"Requested format {fmt} is not available as provider of "
+                "this volume. See `volume.formats` for possible options."
+            )
+        providers = [self._providers[fmt]] if fmt else self._providers.values()
+        for provider in providers:
             try:
-                bbox = provider.get_boundingbox()
+                bbox = provider.get_boundingbox(
+                    clip=clip, background=background, **fetch_kwargs
+                )
                 if bbox.space is None:  # provider does usually not know the space!
                     bbox._space_cached = self.space
                     bbox.minpoint._space_cached = self.space
@@ -230,7 +260,7 @@ def intersection(self, other: structure.BrainStructure, **kwargs) -> structure.B
                 return None  # BrainStructure.intersects check for not None
             return result[0] if len(result) == 1 else result  # if PointSet has single point return as a Point
         elif isinstance(other, boundingbox.BoundingBox):
-            return self.boundingbox.intersection(other)
+            return self.get_boundingbox(clip=True, background=0.0, **kwargs).intersection(other)
         elif isinstance(other, Volume):
             format = kwargs.pop('format', 'image')
             v1 = self.fetch(format=format, **kwargs)
