add mesh-based layer and generic property plugins

setup.py

@@ -42,6 +42,8 @@ setup_kwds = dict(
         ],
         'cellImageQuantification': [
             'surfaceMeshCellCurvature = gnomon_package_tissueimagemesh.algorithm.cellImageQuantification.surfaceMeshCellCurvature',
+            'surfaceMeshCellLayer = gnomon_package_tissueimagemesh.algorithm.cellImageQuantification.surfaceMeshCellLayer',
+            'surfaceMeshCellProperty = gnomon_package_tissueimagemesh.algorithm.cellImageQuantification.surfaceMeshCellProperty'
         ],
         'meshFromImage': [
             'imageSurfaceMesh = gnomon_package_tissueimagemesh.algorithm.meshFromImage.imageSurfaceMesh',

src/gnomon_package_tissueimagemesh/algorithm/cellImageQuantification/surfaceMeshCellLayer.py

+import logging
+from copy import deepcopy
+
+import numpy as np
+import pandas as pd
+
+from dtkcore import d_real, d_inliststring
+
+from gnomon.utils import algorithmPlugin
+from gnomon.utils.decorators import cellImageInput, meshInput, cellImageOutput, dataFrameOutput
+from gnomon.core import gnomonAbstractCellImageQuantification
+
+from timagetk import TissueImage3D
+from timagetk.features.pandas_tools import cells_to_dataframe
+
+from timagetk_geometry.image_surface.tissue_image_mesh import tissue_image_surface_topomesh
+from timagetk_geometry.features.tissue_mesh import compute_cell_layer_feature
+
+
+@algorithmPlugin(version="0.3.0", coreversion="0.81.0", name="Surface Mesh Layer")
+@meshInput('surface_topomesh', data_plugin="gnomonMeshDataPropertyTopomesh")
+@cellImageInput("tissue", data_plugin="gnomonCellImageDataTissueImage")
+@cellImageOutput("out_tissue", data_plugin="gnomonCellImageDataTissueImage")
+@dataFrameOutput('df', data_plugin="gnomonDataFrameDataPandas")
+class surfaceMeshCellLayer(gnomonAbstractCellImageQuantification):
+    """Estimate the cell layer to which each cell in the image belong.
+
+    The method estimates the first layer of cells using a triangle mesh of the
+    tissue surface by selecting the nearest cells to at least one vertex of
+    the mesh. The layer information is then propagated to the rest of the
+    tissue through cell adjacencies.
+
+    """
+
+    def __init__(self):
+        super().__init__()
+
+        self.tissue = {}
+        self.surface_topomesh = {}
+
+        self.out_tissue = {}
+        self.df = {}
+
+        self._parameters = {}
+        self._parameters['orientation'] = d_inliststring('Orientation', "up", ["up", "down"], "Whether to keep upper or lower part")
+        self._parameters['resampling_voxelsize'] = d_real("Resampling voxelsize", 1., 0.1, 5., 1, "Voxelsize to use for the resampling allowing to estimate the surface cell triangulation")
+        self._parameters['down_facing_threshold'] = d_real("Down threshold", 0., -1, 1., 2, "Angle threshold under with the surface (pointing downwards) will be removed")
+
+        self._parameter_groups = {}
+        for parameter_name in ['orientation','resampling_voxelsize', 'down_facing_threshold']:
+            self._parameter_groups[parameter_name] = 'surface_extraction'
+
+    def __del__(self):
+        pass
+
+    def run(self):
+        self.out_tissue = {}
+        self.df = {}
+
+        if len(self.tissue) == 0:
+            logging.error("Impossible to quantify! No segmented image provided")
+            return
+
+        for time in self.tissue.keys():
+            tissue = self.tissue[time]
+
+            out_tissue = TissueImage3D(tissue, not_a_label=0, background=1)
+            out_tissue.cells = deepcopy(tissue.cells)
+            out_tissue.cells.image = out_tissue
+
+            if time in self.surface_topomesh:
+                surface_topomesh = self.surface_topomesh[time]
+            else:
+                surface_topomesh = tissue_image_surface_topomesh(
+                    out_tissue,
+                    resampling_voxelsize=self['resampling_voxelsize'],
+                    surface_matching='cell',
+                    orientation=(1 if self['orientation'] == "up" else -1),
+                    down_facing_threshold=self['down_facing_threshold']
+                )
+            compute_cell_layer_feature(out_tissue, method='surface', surface_topomesh=surface_topomesh)
+
+            self.out_tissue[time] = out_tissue
+
+            self.df[time] = cells_to_dataframe(out_tissue.cells, cell_ids=out_tissue.cell_ids())
+            self.df[time]['time'] = time
+
+        if len(self.df) > 0:
+            self.df = {0: pd.concat([self.df[time] for time in self.df.keys()])}

src/gnomon_package_tissueimagemesh/algorithm/cellImageQuantification/surfaceMeshCellProperty.py

+from copy import deepcopy
+
+import numpy as np
+import scipy.ndimage as nd
+import pandas as pd
+
+from dtkcore import d_inliststringlist, d_int
+
+import gnomon.core
+from gnomon.utils import algorithmPlugin
+from gnomon.utils.decorators import cellImageInput, meshInput, cellImageOutput, dataFrameOutput
+from gnomon.core import gnomonAbstractCellImageQuantification
+
+from timagetk import TissueImage3D
+from timagetk.features.pandas_tools import cells_to_dataframe
+
+from timagetk_geometry.features.cells import compute_cell_surface_center_feature
+
+@algorithmPlugin(version="0.3.0", coreversion="0.81.0", name="Surface Mesh Property")
+@meshInput('surface_topomesh', data_plugin="gnomonMeshDataPropertyTopomesh")
+@cellImageInput("tissue", data_plugin="gnomonCellImageDataTissueImage")
+@cellImageOutput('out_tissue', data_plugin="gnomonCellImageDataTissueImage")
+@dataFrameOutput('df', data_plugin="gnomonDataFrameDataPandas")
+class surfaceMeshCellProperty(gnomonAbstractCellImageQuantification):
+    """Transfer vertex attributes of a surface mesh to the nearest cells.
+    """
+
+    def __init__(self):
+        super().__init__()
+
+        self.surface_topomesh = {}
+        self.tissue = {}
+        self.out_tissue = {}
+        self.df = {}
+
+        self._parameters = {}
+        self._parameters['attribute_names'] = d_inliststringlist("Attributes", [""], [""], "List of vertex attributes to transfer from the surface mesh")
+
+    def refreshParameters(self):
+        if len(self.surface_topomesh)>0:
+            surface_topomesh = list(self.surface_topomesh.values())[0]
+
+            attribute_names = self['attribute_names']
+            property_names = list(surface_topomesh.wisp_property_names(0))
+            if len(property_names) == 0:
+                property_names = [""]
+            attribute_names = [p for p in attribute_names if p in property_names]
+            if len(attribute_names) == 0:
+                attribute_names = property_names
+            self._parameters['attribute_names'].setValues(property_names)
+            self._parameters['attribute_names'].setValue(attribute_names)
+
+    def run(self):
+        self.out_tissue = {}
+        self.df = {}
+
+        for time in self.tissue.keys():
+            tissue = self.tissue[time]
+
+            surface_topomesh = self.surface_topomesh[time]
+
+            out_tissue = TissueImage3D(tissue, not_a_label=0, background=1)
+            out_tissue.cells = deepcopy(tissue.cells)
+            out_tissue.cells.image = out_tissue
+
+            cell_labels = out_tissue.cell_ids()
+            if surface_topomesh.has_wisp_property('cell', 0, is_computed=True):
+                vertex_cells = surface_topomesh.wisp_property('cell', 0).values(list(surface_topomesh.wisps(0)))
+            else:
+                if 'surface_center' in out_tissue.cells.feature_names():
+                    cell_surface_centers = out_tissue.cells.feature('surface_center')
+                else:
+                    cell_surface_centers = compute_cell_surface_center_feature(out_tissue)
+                cell_centers = np.array([cell_surface_centers[c] for c in cell_labels])
+
+                vertex_points = surface_topomesh.wisp_property('barycenter', 0).values(list(surface_topomesh.wisps(0)))
+                vertex_cell_distances = np.linalg.norm(vertex_points[:, np.newaxis] - cell_centers[np.newaxis], axis=-1)
+                vertex_cells = cell_labels[np.argmin(vertex_cell_distances, axis=-1)]
+
+                surface_topomesh.update_wisp_property('cell', 0, dict(zip(surface_topomesh.wisps(0), vertex_cells)))
+                surface_topomesh.update_wisp_property('label', 0, dict(zip(surface_topomesh.wisps(0), vertex_cells%256)))
+
+            for property_name in self['attribute_names']:
+                vertex_property = surface_topomesh.wisp_property(property_name, 0).values(list(surface_topomesh.wisps(0)))
+                property_values = np.unique(vertex_property)
+                property_is_binary = vertex_property.ndim == 1 and len(property_values) <= 2 and (0 in property_values or 1 in property_values)
+                cell_property = None
+                if property_is_binary:
+                    cell_property = (nd.sum(vertex_property, vertex_cells, index=cell_labels) > 0).astype(int)
+                else:
+                    if vertex_property.ndim == 1:
+                        cell_property = nd.mean(vertex_property, vertex_cells, index=cell_labels)
+                    elif vertex_property.ndim == 2:
+                        cell_property = list(np.transpose([
+                            nd.mean(vertex_property[:, k], vertex_cells, index=cell_labels)
+                            for k in range(vertex_property.shape[1])
+                        ]))
+                    elif vertex_property.ndim == 3:
+                        cell_property = list(np.transpose([[
+                            nd.mean(vertex_property[:, j, k], vertex_cells, index=cell_labels)
+                            for k in range(vertex_property.shape[2])
+                        ] for j in range(vertex_property.shape[1])], (2, 1, 0)))
+                if cell_property is not None:
+                    out_tissue.cells.set_feature(property_name, dict(zip(cell_labels, cell_property)))
+
+            self.out_tissue[time] = out_tissue
+
+            self.df[time] = cells_to_dataframe(out_tissue.cells, cell_ids=out_tissue.cell_ids())
+            self.df[time]['time'] = time
+
+        if len(self.df)>0:
+            self.df = {0: pd.concat([self.df[time] for time in self.df.keys()])}

src/gnomon_package_tissueimagemesh/algorithm/pointCloudQuantification/surfaceMeshProperty.py

@@ -90,9 +90,6 @@ class surfaceMeshProperty(gnomonAbstractPointCloudQuantification):
                             for k in range(vertex_property.shape[2])
                         ] for j in range(vertex_property.shape[1])], (2, 1, 0)))
 
-            for i_dim, dim in enumerate(['x', 'y', 'z']):
-                surface_topomesh.update_wisp_property('barycenter_' + dim, 0, surface_topomesh.wisp_property('barycenter', 0).values()[:, i_dim])
-
             self.out_df[time] = out_df
 
             self.data_df[time] = deepcopy(out_df)