diff --git a/siibra/features/tabular/cell_density_profile.py b/siibra/features/tabular/cell_density_profile.py index fea5af7a..257755dd 100644 --- a/siibra/features/tabular/cell_density_profile.py +++ b/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,6 +250,7 @@ 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: @@ -249,16 +258,5 @@ def _values(self): 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 - )) diff --git a/siibra/features/tabular/layerwise_cell_density.py b/siibra/features/tabular/layerwise_cell_density.py index 1c71636c..d5f2e3c8 100644 --- a/siibra/features/tabular/layerwise_cell_density.py +++ b/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 - ))