Merge branch 'feature/wall-mesh-quantification' into 'develop'

Towards PIN1 quantification

Closes gnomon/gnomon#652

See merge request !6

setup.py

@@ -43,11 +43,14 @@ 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'
+            'surfaceMeshCellProperty = gnomon_package_tissueimagemesh.algorithm.cellImageQuantification.surfaceMeshCellProperty',
+            'wallMeshImageSignal = gnomon_package_tissueimagemesh.algorithm.cellImageQuantification.wallMeshImageSignal',
+            'wallMeshSignalPolarities = gnomon_package_tissueimagemesh.algorithm.cellImageQuantification.wallMeshSignalPolarities',
         ],
         'meshFromImage': [
             'imageSurfaceMesh = gnomon_package_tissueimagemesh.algorithm.meshFromImage.imageSurfaceMesh',
-            'segmentationSurfaceMesh = gnomon_package_tissueimagemesh.algorithm.meshFromImage.segmentationSurfaceMesh'
+            'segmentationSurfaceMesh = gnomon_package_tissueimagemesh.algorithm.meshFromImage.segmentationSurfaceMesh',
+            'segmentationWallMesh = gnomon_package_tissueimagemesh.algorithm.meshFromImage.segmentationWallMesh'
         ],
         'pointCloudQuantification': [
             'nucleiCurvatureNukem = gnomon_package_tissueimagemesh.algorithm.pointCloudQuantification.nucleiCurvatureNukem',

src/gnomon_package_tissueimagemesh/algorithm/cellImageQuantification/wallMeshImageSignal.py

+import logging
+from copy import deepcopy
+
+import numpy as np
+import pandas as pd
+
+from dtkcore import d_real, d_inliststring, d_inliststringlist
+
+from gnomon.utils import algorithmPlugin
+from gnomon.utils.decorators import cellImageInput, imageInput, meshInput, cellImageOutput, dataFrameOutput
+from gnomon.core import gnomonAbstractCellImageQuantification
+
+from timagetk import TissueImage3D
+from timagetk.features.pandas_tools import cells_to_dataframe
+
+from timagetk_geometry.signal_quantification.wall_mesh import quantify_wall_topomesh_signal_intensity
+
+
+@algorithmPlugin(version="0.3.0", coreversion="0.81.0", name="Wall Mesh Signal")
+@meshInput('wall_topomesh', data_plugin="gnomonMeshDataPropertyTopomesh")
+@imageInput('img', data_plugin='gnomonImageDataMultiChannelImage')
+@cellImageInput("tissue", data_plugin="gnomonCellImageDataTissueImage")
+@cellImageOutput("out_tissue", data_plugin="gnomonCellImageDataTissueImage")
+@dataFrameOutput('df', data_plugin="gnomonDataFrameDataPandas")
+class wallMeshImageSignal(gnomonAbstractCellImageQuantification):
+    """Measure cell-wall image signal intensity using a triangle mesh.
+
+    The method computes a value of signal for each cell wall represented in a
+    wall triangle mesh by averaging the image intensity at the level of the
+    mesh vertices.
+
+    """
+
+    def __init__(self):
+        super().__init__()
+
+        self.tissue = {}
+        self.img = {}
+        self.wall_topomesh = {}
+
+        self.out_tissue = {}
+        self.df = {}
+
+        self._parameters = {}
+        self._parameters['channels'] = d_inliststringlist("Channels", [""],[""], "Channels on which to quantify signal intensity at the cell walls")
+        self._parameters['gaussian_sigma'] = d_real("Sigma", 0.5, 0., 50., 2, "Standard deviation of the Gaussian kernel used to smooth signal")
+
+        self._parameter_groups = {}
+
+    def __del__(self):
+        pass
+
+    def refreshParameters(self):
+        if len(self.img)>0:
+            img = list(self.img.values())[0]
+
+            if len(img.channel_names) == 1:
+                if 'channels' in self._parameters.keys():
+                    del self._parameters['channels']
+            else:
+                if not 'channels' in self._parameters.keys():
+                    self._parameters['channels'] = d_inliststringlist("channels", [""], [""], "Channels on which to quantify signal intensity at the cell walls")
+                self._parameters['channels'].setValues(img.channel_names)
+                self._parameters['channels'].setValue(img.channel_names)
+
+    def run(self):
+        self.out_tissue = {}
+        self.df = {}
+
+        for time in self.tissue.keys():
+            if not time in self.img.keys():
+                logging.error("Impossible to quantify! No signal image provided")
+                self.df = {}
+                return
+            if not time in self.wall_topomesh.keys():
+                logging.error("Impossible to quantify! No wall mesh provided")
+                self.df = {}
+                return
+            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
+            out_tissue.walls = deepcopy(tissue.walls)
+            out_tissue.walls.image = out_tissue
+
+            img = self.img[time]
+            wall_topomesh = self.wall_topomesh[time]
+
+            wall_signals = quantify_wall_topomesh_signal_intensity(
+                signal_img=img,
+                wall_topomesh=wall_topomesh,
+                wall_sigma=self['gaussian_sigma']
+            )
+
+            for channel in wall_signals:
+                signal_name = channel if channel != "" else "image_signal"
+                if ('channels' not in self._parameters) or (signal_name in self['channels']):
+                    out_tissue.walls.set_feature(signal_name, wall_signals[channel])
+
+            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/wallMeshSignalPolarities.py

+import logging
+from copy import deepcopy
+
+import numpy as np
+import pandas as pd
+
+from dtkcore import d_real, d_inliststring, d_inliststringlist
+
+from gnomon.utils import algorithmPlugin
+from gnomon.utils.decorators import cellImageInput, imageInput, meshInput, cellImageOutput, dataFrameOutput
+from gnomon.core import gnomonAbstractCellImageQuantification
+
+from timagetk import TissueImage3D
+from timagetk.features.pandas_tools import cell_walls_to_dataframe
+
+from timagetk_geometry.signal_quantification.wall_mesh import quantify_wall_topomesh_membrane_signal_intensity, compute_wall_property
+
+
+@algorithmPlugin(version="0.3.0", coreversion="0.81.0", name="Wall Signal Polarities")
+@meshInput('wall_topomesh', data_plugin="gnomonMeshDataPropertyTopomesh")
+@imageInput('img', data_plugin='gnomonImageDataMultiChannelImage')
+@cellImageInput("tissue", data_plugin="gnomonCellImageDataTissueImage")
+@cellImageOutput("out_tissue", data_plugin="gnomonCellImageDataTissueImage")
+@dataFrameOutput('df', data_plugin="gnomonDataFrameDataPandas")
+class wallMeshSignalPolarities(gnomonAbstractCellImageQuantification):
+    """Measure image signal polarities using cell wall meshes.
+
+    The method computes a value of signal on both sides of each cell wall
+    represented in a wall triangle mesh and assess the polarity of the signal
+    at the level of each cell-cell interface.
+
+    """
+
+    def __init__(self):
+        super().__init__()
+
+        self.tissue = {}
+        self.img = {}
+        self.wall_topomesh = {}
+
+        self.out_tissue = {}
+        self.df = {}
+
+        self._parameters = {}
+        self._parameters['membrane_channel'] = d_inliststring("Reference channel", "",[""], "Channel to use as a reference for positioning the cell wall")
+        self._parameters['channels'] = d_inliststringlist("Channels", [""],[""], "Channels on which to quantify signal polarity at the cell walls")
+        self._parameters['wall_distance'] = d_real("Wall distance", 0.6, 0., 10., 2, "Distance used to average the signal at each wall vertex")
+        self._parameters['max_distance'] = d_real("Max distance", 1.5, 0., 20., 2, "Maximal distance when projecting image signal on normal axes")
+        self._parameters['max_radius'] = d_real("Max radius", 1, 0., 20., 2, "Maximal radius when projecting image signal on normal axes")
+
+        self._parameter_groups = {}
+
+    def __del__(self):
+        pass
+
+    def refreshParameters(self):
+        if len(self.img)>0:
+            img = list(self.img.values())[0]
+
+            if len(img.channel_names) == 1:
+                if 'membrane_channel' in self._parameters.keys():
+                    del self._parameters['membrane_channel']
+                if 'channels' in self._parameters.keys():
+                    del self._parameters['channels']
+            else:
+                if 'membrane_channel' not in self._parameters.keys():
+                    self._parameters['membrane_channel'] = d_inliststring("Reference channel", "",[""], "Channel to use as a reference for positioning the cell wall")
+
+                if 'channels' not in self._parameters.keys():
+                    self._parameters['channels'] = d_inliststringlist("channels", [""], [""], "Channels on which to quantify signal polarity at the cell walls")
+
+            membrane_channel = self['membrane_channel']
+            self._parameters['membrane_channel'].setValues(img.channel_names)
+            membrane_channel = membrane_channel if membrane_channel in img.channel_names else img.channel_names[0]
+            self._parameters['membrane_channel'].setValue(membrane_channel)
+
+            channel_names = self['channels']
+            self._parameters['channels'].setValues(img.channel_names)
+            channel_names = [c for c in channel_names if c in img.channel_names]
+            if len(channel_names) == 0:
+                channel_names = img.channel_names
+            self._parameters['channels'].setValue(channel_names)
+
+    def run(self):
+        self.out_tissue = {}
+        self.df = {}
+
+        for time in self.tissue.keys():
+            if not time in self.img.keys():
+                logging.error("Impossible to quantify! No signal image provided")
+                self.df = {}
+                return
+            if not time in self.wall_topomesh.keys():
+                logging.error("Impossible to quantify! No wall mesh provided")
+                self.df = {}
+                return
+            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
+            out_tissue.walls = deepcopy(tissue.walls)
+            out_tissue.walls.image = out_tissue
+
+            img = self.img[time]
+            wall_topomesh = self.wall_topomesh[time]
+
+            cell_centers = out_tissue.cells.barycenter()
+
+            membrane_channel = "" if 'channels' not in self._parameters else self['membrane_channel']
+            channel_names = [""] if 'channels' not in self._parameters else self['channels']
+
+            wall_areas = compute_wall_property(wall_topomesh, 'area')
+            wall_centers = compute_wall_property(wall_topomesh, 'barycenter')
+            wall_normals = compute_wall_property(wall_topomesh, 'normal')
+
+            out_tissue.walls.set_feature('area', wall_areas)
+            for k, dim in enumerate('xyz'):
+                out_tissue.walls.set_feature(f'center_{dim}', {w: c[k] for w, c in wall_centers.items()})
+            for k, dim in enumerate('xyz'):
+                out_tissue.walls.set_feature(f'normal_{dim}', {w: n[k] for w, n in wall_normals.items()})
+
+            wall_signals = quantify_wall_topomesh_membrane_signal_intensity(
+                signal_img=img,
+                all_wall_topomesh=wall_topomesh,
+                cell_centers=cell_centers,
+                membrane_channel=membrane_channel,
+                channel_names=channel_names,
+                compute_polarities=True,
+                wall_distance=self['wall_distance'],
+                wall_sigma=self['max_distance'],
+                line_sigma=self['max_radius']
+            )
+
+            for channel_name in channel_names:
+                signal_name = channel_name if channel_name != "" else "image_signal"
+                for suffix in ["", "_left", "_right", "_polarity"]:
+                    out_tissue.walls.set_feature(f"{signal_name}{suffix}", wall_signals[f"{channel_name}{suffix}"])
+
+            self.out_tissue[time] = out_tissue
+
+            self.df[time] = cell_walls_to_dataframe(out_tissue.walls, 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/meshFromImage/segmentationWallMesh.py

+import numpy as np
+
+
+from dtkcore import d_bool, d_int, d_inliststringlist, d_real
+
+import gnomon.core
+from gnomon.core import gnomonAbstractMeshFromImage
+
+from gnomon.utils import corePlugin
+from gnomon.utils.decorators import cellImageInput, meshOutput
+
+from timagetk_geometry.image_surface.wall_mesh import tissue_image_wall_topomesh
+
+from cellcomplex.property_topomesh.analysis import compute_topomesh_property
+
+
+@corePlugin(version="0.3.0", coreversion="0.81.0", name="Cell Wall Meshes")
+@cellImageInput("seg_img", data_plugin='gnomonCellImageDataTissueImage')
+@meshOutput("topomesh", data_plugin="gnomonMeshDataPropertyTopomesh")
+class segmentationWallMesh(gnomonAbstractMeshFromImage):
+    """Extract the cell walls of a segmented image as triangular meshes.
+    """
+
+    def __init__(self):
+        super().__init__()
+
+        self.seg_img = {}
+        self.topomesh = {}
+
+        self._parameters = {}
+        self._parameters['resampling_voxelsize'] = d_real("Resampling voxelsize", 0, 0, 1, 2, "Cubic voxelsize in which to resample the image before mesh extraction")
+        self._parameters['smoothing'] = d_bool("Smoothing", True, "Whether to smooth and decimate the wall meshes")
+        self._parameters['wall_types'] = d_inliststringlist('Walls', ["anticlinal_L1"], ["anticlinal_L1", "epidermis_L1", "all"], "Which cell walls to extract in the resulting mesh")
+        self._parameters['edge_length'] = d_real('Edge length', 1, 0, 10, 1, "Target length of mesh edges")
+
+        self._parameter_groups = {}
+
+    def refreshParameters(self):
+        if len(self.seg_img)>0:
+            seg_img = list(self.seg_img.values())[0]
+
+            self._parameters['resampling_voxelsize'].setMin(0)
+            self._parameters['resampling_voxelsize'].setMax(10*np.max(seg_img.voxelsize))
+            self._parameters['resampling_voxelsize'].setValue(4*np.mean(seg_img.voxelsize))
+
+    def run(self):
+        self.topomesh = {}
+
+        for time in self.seg_img.keys():
+            seg_img = self.seg_img[time]
+
+            resampling_voxelsize = None if self['resampling_voxelsize'] == 0 else self['resampling_voxelsize']
+            topomesh = tissue_image_wall_topomesh(
+                tissue=seg_img,
+                resampling_voxelsize=resampling_voxelsize,
+                wall_types=self['wall_types'],
+                smoothing=self['smoothing'],
+                target_edge_length=self['edge_length']
+            )
+
+            for degree in [1, 2, 3]:
+                compute_topomesh_property(topomesh, 'barycenter', degree)
+
+            for degree in [0, 1, 2, 3]:
+                positions = topomesh.wisp_property('barycenter', degree)
+                for k, dim in enumerate(['x', 'y', 'z']):
+                    topomesh.update_wisp_property('barycenter_'+dim, degree, dict(zip(positions.keys(), positions.values()[:, k])))
+
+            self.topomesh[time] = topomesh