fix: cell density profile

siibra/features/tabular/cell_density_profile.py

@@ -16,15 +16,35 @@
 from . import cortical_profile
 
 from .. import anchor as _anchor
-from ...commons import PolyLine, logger, create_key
+from ...commons import PolyLine, logger
 from ...retrieval import requests
 
 from skimage.draw import polygon
 from skimage.transform import resize
-from io import BytesIO
 import numpy as np
 import pandas as pd
 
+from io import BytesIO
+from typing_extensions import Tuple
+
+
+def cell_reader(bytes_buffer: bytes):
+    return pd.read_csv(BytesIO(bytes_buffer[2:]), delimiter=" ", header=0).astype(
+        {"layer": int, "label": int}
+    )
+
+
+def layer_reader(bytes_buffer: bytes):
+    return pd.read_csv(BytesIO(bytes_buffer[2:]), delimiter=" ", header=0, index_col=0)
+
+
+def poly_srt(poly):
+    return poly[poly[:, 0].argsort(), :]
+
+
+def poly_rev(poly):
+    return poly[poly[:, 0].argsort()[::-1], :]
+
 
 class CellDensityProfile(
     cortical_profile.CorticalProfile,
@@ -45,24 +65,6 @@ class CellDensityProfile(
 
     _filter_attrs = cortical_profile.CorticalProfile._filter_attrs + ["location"]
 
-    @classmethod
-    def CELL_READER(cls, b):
-        return pd.read_csv(BytesIO(b[2:]), delimiter=" ", header=0).astype(
-            {"layer": int, "label": int}
-        )
-
-    @classmethod
-    def LAYER_READER(cls, b):
-        return pd.read_csv(BytesIO(b[2:]), delimiter=" ", header=0, index_col=0)
-
-    @staticmethod
-    def poly_srt(poly):
-        return poly[poly[:, 0].argsort(), :]
-
-    @staticmethod
-    def poly_rev(poly):
-        return poly[poly[:, 0].argsort()[::-1], :]
-
     def __init__(
         self,
         section: int,
@@ -89,9 +91,9 @@ def __init__(
         )
         self._step = 0.01
         self._url = url
-        self._cell_loader = requests.HttpRequest(url, self.CELL_READER)
+        self._cell_loader = requests.HttpRequest(url, cell_reader)
         self._layer_loader = requests.HttpRequest(
-            url.replace("segments", "layerinfo"), self.LAYER_READER
+            url.replace("segments", "layerinfo"), layer_reader
         )
         self._density_image = None
         self._layer_mask = None
@@ -105,66 +107,74 @@ def location(self):
 
     @property
     def shape(self):
-        return tuple(self.cells[["y", "x"]].max().astype("int") + 1)
+        """(y,x)"""
+        return tuple(np.ceil(self.cells[["y", "x"]].max()).astype("int"))
 
-    def boundary_annotation(self, boundary):
+    def boundary_annotation(self, boundary: Tuple[int, int]) -> np.ndarray:
         """Returns the annotation of a specific layer boundary."""
-        y1, x1 = self.shape
+        shape_y, shape_x = self.shape
 
         # start of image patch
         if boundary == (-1, 0):
-            return np.array([[0, 0], [x1, 0]])
+            return np.array([[0, 0], [shape_x, 0]])
 
         # end of image patch
         if boundary == (7, 8):
-            return np.array([[0, y1], [x1, y1]])
+            return np.array([[0, shape_y], [shape_x, shape_y]])
 
         # retrieve polygon
         basename = "{}_{}.json".format(
             *(self.LAYERS[layer] for layer in boundary)
         ).replace("0_I", "0")
-        url = self._url.replace("segments.txt", basename)
-        poly = self.poly_srt(np.array(requests.HttpRequest(url).get()["segments"]))
+        poly_url = self._url.replace("segments.txt", basename)
+        poly = poly_srt(np.array(requests.HttpRequest(poly_url).get()["segments"]))
 
         # ensure full width
         poly[0, 0] = 0
-        poly[-1, 0] = x1
+        poly[-1, 0] = shape_x
 
         return poly
 
-    def layer_annotation(self, layer):
+    def layer_annotation(self, layer: int) -> np.ndarray:
         return np.vstack(
             (
                 self.boundary_annotation((layer - 1, layer)),
-                self.poly_rev(self.boundary_annotation((layer, layer + 1))),
+                poly_rev(self.boundary_annotation((layer, layer + 1))),
                 self.boundary_annotation((layer - 1, layer))[0, :],
             )
         )
 
     @property
-    def layer_mask(self):
+    def layer_mask(self) -> np.ndarray:
         """Generates a layer mask from boundary annotations."""
         if self._layer_mask is None:
-            self._layer_mask = np.zeros(np.array(self.shape).astype("int") + 1)
+            self._layer_mask = np.zeros(self.shape, dtype="int")
             for layer in range(1, 8):
                 pl = self.layer_annotation(layer)
                 X, Y = polygon(pl[:, 0], pl[:, 1])
-                self._layer_mask[Y, X] = layer
+                # remove pixels outside the image since the points may be delivered outside of bounds
+                valid_voxels = np.logical_and.reduce((
+                    X < self.shape[1],
+                    Y < self.shape[0],
+                    X >= 0,
+                    Y >= 0,
+                ))
+                self._layer_mask[Y[valid_voxels], X[valid_voxels]] = layer
         return self._layer_mask
 
     @property
-    def depth_image(self):
+    def depth_image(self) -> np.ndarray:
         """Cortical depth image from layer boundary polygons by equidistant sampling."""
 
         if self._depth_image is None:
-
+            logger.info("Calculating densities from cell and layer data...")
             # compute equidistant cortical depth image from inner and outer contour
             scale = 0.1
-            D = np.zeros((np.array(self.density_image.shape) * scale).astype("int") + 1)
+            depth_arr = np.zeros(np.ceil(np.array(self.shape) * scale + 1).astype("int"))
 
             # determine sufficient stepwidth for profile sampling
             # to match downscaled image resolution
-            vstep, hstep = 1.0 / np.array(D.shape) / 2.0
+            vstep, hstep = 1.0 / np.array(depth_arr.shape) / 2.0
             vsteps = np.arange(0, 1 + vstep, vstep)
             hsteps = np.arange(0, 1 + vstep, hstep)
 
@@ -176,15 +186,15 @@ def depth_image(self):
             # write sample depths to their location in the depth image
             for prof in profiles:
                 XY = prof.sample(vsteps).astype("int")
-                D[XY[:, 1], XY[:, 0]] = vsteps
+                depth_arr[XY[:, 1], XY[:, 0]] = vsteps
 
             # fix wm region, account for rounding error
-            XY = self.layer_annotation(7) * scale
-            D[polygon(XY[:, 1] - 1, XY[:, 0])] = 1
-            D[-1, :] = 1
+            XY = np.round(self.layer_annotation(7) * scale)
+            depth_arr[polygon(XY[:, 1] - 1, XY[:, 0])] = 1
+            depth_arr[-1, :] = 1
 
             # rescale depth image to original patch size
-            self._depth_image = resize(D, self.density_image.shape)
+            self._depth_image = resize(depth_arr, np.array(self.shape) + 1)
 
         return self._depth_image
 
@@ -200,7 +210,7 @@ def boundary_positions(self):
         return self._boundary_positions
 
     @property
-    def density_image(self):
+    def density_image(self) -> np.ndarray:
         if self._density_image is None:
             logger.debug("Computing density image for", self._url)
             # we integrate cell counts into 2D bins
@@ -209,9 +219,7 @@ def density_image(self):
             counts, xedges, yedges = np.histogram2d(
                 self.cells.y,
                 self.cells.x,
-                bins=(np.array(self.layer_mask.shape) / pixel_size_micron + 0.5).astype(
-                    "int"
-                ),
+                bins=np.ceil(np.array(self.shape) / pixel_size_micron).astype("int"),
             )
 
             # rescale the counts from count / pixel_size**2  to count / 0.1mm^3,
@@ -224,16 +232,16 @@ def density_image(self):
             counts /= np.cbrt(self.BIGBRAIN_VOLUMETRIC_SHRINKAGE_FACTOR) ** 2
 
             # to go to 0.1 millimeter cube, we multiply by 0.1 / 0.0002 = 500
-            self._density_image = resize(counts, self.layer_mask.shape, order=2)
+            self._density_image = resize(counts, np.array(self.shape) + 1, order=2)
 
         return self._density_image
 
     @property
-    def cells(self):
+    def cells(self) -> pd.DataFrame:
         return self._cell_loader.get()
 
     @property
-    def layers(self):
+    def layers(self) -> pd.DataFrame:
         return self._layer_loader.get()
 
     @property
@@ -242,23 +250,13 @@ def _depths(self):
 
     @property
     def _values(self):
+        # TODO: release a dataset update instead of on the fly computation
         densities = []
         delta = self._step / 2.0
         for d in self._depths:
             mask = (self.depth_image >= d - delta) & (self.depth_image < d + delta)
             if np.sum(mask) > 0:
                 densities.append(self.density_image[mask].mean())
             else:
-                densities.append(np.NaN)
+                densities.append(np.nan)
         return np.asanyarray(densities)
-
-    @property
-    def key(self):
-        assert len(self.species) == 1
-        return create_key("{}_{}_{}_{}_{}".format(
-            self.id,
-            self.species[0]['name'],
-            self.regionspec,
-            self.section,
-            self.patch
-        ))

siibra/features/tabular/layerwise_cell_density.py

@@ -16,13 +16,13 @@
 from . import cortical_profile
 from .. import anchor as _anchor
 from . import tabular
+from ..tabular.cell_density_profile import cell_reader, layer_reader
 
 from ... import commons
 from ...retrieval import requests
 
 import pandas as pd
 import numpy as np
-from io import BytesIO
 
 
 class LayerwiseCellDensity(
@@ -40,16 +40,6 @@ class LayerwiseCellDensity(
         "The cortical depth is estimated from the measured layer thicknesses."
     )
 
-    @classmethod
-    def CELL_READER(cls, b):
-        return pd.read_csv(BytesIO(b[2:]), delimiter=" ", header=0).astype(
-            {"layer": int, "label": int}
-        )
-
-    @classmethod
-    def LAYER_READER(cls, b):
-        return pd.read_csv(BytesIO(b[2:]), delimiter=" ", header=0, index_col=0)
-
     def __init__(
         self,
         segmentfiles: list,
@@ -77,8 +67,8 @@ def _load_densities(self):
         density_dict = {}
         for i, (cellfile, layerfile) in enumerate(self._filepairs):
             try:
-                cells = requests.HttpRequest(cellfile, func=self.CELL_READER).data
-                layers = requests.HttpRequest(layerfile, func=self.LAYER_READER).data
+                cells = requests.HttpRequest(cellfile, func=cell_reader).data
+                layers = requests.HttpRequest(layerfile, func=layer_reader).data
             except requests.SiibraHttpRequestError as e:
                 print(str(e))
                 commons.logger.error(f"Skipping to bootstrap a {self.__class__.__name__} feature, cannot access file resource.")
@@ -103,12 +93,3 @@ def data(self):
             )
             self._data_cached.index.name = 'layer'
         return self._data_cached
-
-    @property
-    def key(self):
-        assert len(self.species) == 1
-        return commons.create_key("{}_{}_{}".format(
-            self.dataset_id,
-            self.species[0]['name'],
-            self.regionspec
-        ))