diff --git a/setup.py b/setup.py
index 8fc8b2932d3a8d755dfe7f0c4ef4e332dce2941b..3624595a7dca53aacced312e24aee190cfdaf7dc 100644
--- a/setup.py
+++ b/setup.py
@@ -42,13 +42,17 @@ 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'
+ 'imageSurfaceMesh = gnomon_package_tissueimagemesh.algorithm.meshFromImage.imageSurfaceMesh',
+ 'segmentationSurfaceMesh = gnomon_package_tissueimagemesh.algorithm.meshFromImage.segmentationSurfaceMesh'
],
'pointCloudQuantification': [
'nucleiCurvatureNukem = gnomon_package_tissueimagemesh.algorithm.pointCloudQuantification.nucleiCurvatureNukem',
- 'nucleiLayerEstimation = gnomon_package_tissueimagemesh.algorithm.pointCloudQuantification.nucleiLayerEstimation'
+ 'nucleiLayerEstimation = gnomon_package_tissueimagemesh.algorithm.pointCloudQuantification.nucleiLayerEstimation',
+ 'surfaceMeshProperty = gnomon_package_tissueimagemesh.algorithm.pointCloudQuantification.surfaceMeshProperty'
],
},
keywords='',
diff --git a/src/gnomon_package_tissueimagemesh/algorithm/cellImageQuantification/surfaceMeshCellCurvature.py b/src/gnomon_package_tissueimagemesh/algorithm/cellImageQuantification/surfaceMeshCellCurvature.py
index b291c79760cb291a430b2b4109a01e85a840c0a1..c76a31128b623cd6e7975c89186c1d5169c63a47 100644
--- a/src/gnomon_package_tissueimagemesh/algorithm/cellImageQuantification/surfaceMeshCellCurvature.py
+++ b/src/gnomon_package_tissueimagemesh/algorithm/cellImageQuantification/surfaceMeshCellCurvature.py
@@ -7,7 +7,7 @@ import pandas as pd
from dtkcore import d_real, d_inliststring
from gnomon.utils import algorithmPlugin
-from gnomon.utils.decorators import cellImageInput, cellImageOutput, dataFrameOutput
+from gnomon.utils.decorators import cellImageInput, meshInput, cellImageOutput, dataFrameOutput
from gnomon.core import gnomonAbstractCellImageQuantification
from timagetk import TissueImage3D
@@ -18,6 +18,7 @@ from timagetk_geometry.features.tissue_mesh import compute_cell_curvature_featur
@algorithmPlugin(version="0.3.0", coreversion="0.81.0", name="Surface Mesh Curvature")
+@meshInput('surface_topomesh', data_plugin="gnomonMeshDataPropertyTopomesh")
@cellImageInput("tissue", data_plugin="gnomonCellImageDataTissueImage")
@cellImageOutput("out_tissue", data_plugin="gnomonCellImageDataTissueImage")
@dataFrameOutput('df', data_plugin="gnomonDataFrameDataPandas")
@@ -34,6 +35,7 @@ class surfaceMeshCellCurvature(gnomonAbstractCellImageQuantification):
super().__init__()
self.tissue = {}
+ self.surface_topomesh = {}
self.out_tissue = {}
self.df = {}
@@ -65,13 +67,16 @@ class surfaceMeshCellCurvature(gnomonAbstractCellImageQuantification):
out_tissue.cells = deepcopy(tissue.cells)
out_tissue.cells.image = out_tissue
- 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']
- )
+ 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_curvature_features(out_tissue, surface_topomesh=surface_topomesh)
self.out_tissue[time] = out_tissue
diff --git a/src/gnomon_package_tissueimagemesh/algorithm/cellImageQuantification/surfaceMeshCellLayer.py b/src/gnomon_package_tissueimagemesh/algorithm/cellImageQuantification/surfaceMeshCellLayer.py
new file mode 100644
index 0000000000000000000000000000000000000000..d65dfa410b80d80c6f6aa19e33721a0b5f4f879c
--- /dev/null
+++ b/src/gnomon_package_tissueimagemesh/algorithm/cellImageQuantification/surfaceMeshCellLayer.py
@@ -0,0 +1,89 @@
+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()])}
diff --git a/src/gnomon_package_tissueimagemesh/algorithm/cellImageQuantification/surfaceMeshCellProperty.py b/src/gnomon_package_tissueimagemesh/algorithm/cellImageQuantification/surfaceMeshCellProperty.py
new file mode 100644
index 0000000000000000000000000000000000000000..2371dc1fba6506139e13da87055350611d41cc79
--- /dev/null
+++ b/src/gnomon_package_tissueimagemesh/algorithm/cellImageQuantification/surfaceMeshCellProperty.py
@@ -0,0 +1,112 @@
+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()])}
diff --git a/src/gnomon_package_tissueimagemesh/algorithm/meshFromImage/imageSurfaceMesh.py b/src/gnomon_package_tissueimagemesh/algorithm/meshFromImage/imageSurfaceMesh.py
index 0cc29b48e19410b77f8cd5b2d45fc7bd98e20887..dd3874814773d1c70beaf4724478effcb22df377 100644
--- a/src/gnomon_package_tissueimagemesh/algorithm/meshFromImage/imageSurfaceMesh.py
+++ b/src/gnomon_package_tissueimagemesh/algorithm/meshFromImage/imageSurfaceMesh.py
@@ -11,10 +11,9 @@ from gnomon.utils import corePlugin
from gnomon.utils.decorators import imageInput, meshOutput
-from timagetk_geometry.image_surface.image_mesh import image_surface_topomesh
+from timagetk_geometry.image_surface.image_mesh import image_surface_topomesh, up_facing_surface_topomesh
from cellcomplex.property_topomesh.analysis import compute_topomesh_property
-from cellcomplex.property_topomesh.extraction import cut_surface_topomesh
from cellcomplex.property_topomesh.creation import triangle_topomesh
from cellcomplex.utils import array_dict
@@ -29,8 +28,10 @@ class imageSurfaceMesh(gnomonAbstractMeshFromImage):
filter of standard deviation gaussian_sigma and the result is thresholded
by the value of the threshold parameter. The obtained binary mask is
meshed using a Marching Cubes algorithm and decimated to reach the target
- edge_length. Only the top or bottom part of the mesh is kept depending on
- the orientation parameter, up to a given cut_level.
+ edge_length. Only the largest connected component o the top or bottom part
+ of the mesh is kept depending on the orientation parameter, after removing
+ the parts beyond a given threshold of vertical orientation of the surface
+ normal vectors.
"""
def __init__(self):
@@ -46,13 +47,13 @@ class imageSurfaceMesh(gnomonAbstractMeshFromImage):
self._parameters['gaussian_sigma'] = d_real('Sigma', 2, 0., 20., 2, "Standard deviation of the Gaussian kernel used to smooth signal")
self._parameters['orientation'] = d_inliststring('Orientation', "up", ["up", "down"], "Whether to keep upper or lower part")
- self._parameters['cut_level'] = d_real('Cut level', 0, 0, 50, 1, "Level of the z cut")
+ self._parameters['down_threshold'] = d_real('Orientation threshold', 0, -1, 1, 1, "The z coordinate below which a normal vector is considered to face downwards")
self._parameters['edge_length'] = d_real('Edge length', 10, 0, 20, 1, "Maximal length of mesh edges")
self._parameter_groups = {}
for parameter_name in ['resampling_voxelsize', 'gaussian_sigma']:
self._parameter_groups[parameter_name] = 'image_preprocessing'
- for parameter_name in ['orientation', 'cut_level']:
+ for parameter_name in ['orientation', 'down_threshold']:
self._parameter_groups[parameter_name] = 'surface_cutting'
def refreshParameters(self):
@@ -79,16 +80,12 @@ class imageSurfaceMesh(gnomonAbstractMeshFromImage):
self._parameters['threshold'].setMax(255)
elif img.dtype == np.uint16:
self._parameters['threshold'].setMax(65535)
- self._parameters['threshold'].setValue(1.2*threshold_otsu(img.get_array()))
+ self._parameters['threshold'].setValue(int(np.round(1.2*threshold_otsu(img.get_array()))))
self._parameters['resampling_voxelsize'].setMin(0)
self._parameters['resampling_voxelsize'].setMax(10*np.max(img.voxelsize))
self._parameters['resampling_voxelsize'].setValue(4*np.mean(img.voxelsize))
- self._parameters['cut_level'].setMin(0.)
- self._parameters['cut_level'].setMax(img.shape[0]*img.voxelsize[0])
- self._parameters['cut_level'].setValue(img.shape[0]*img.voxelsize[0]/2.)
-
def run(self):
self.topomesh = {}
@@ -109,9 +106,13 @@ class imageSurfaceMesh(gnomonAbstractMeshFromImage):
target_length=self['edge_length']
)
- topomesh = cut_surface_topomesh(topomesh, z_cut=self['cut_level'], below=self['orientation']=="up")
- # topomesh = up_facing_surface_topomesh(topomesh, normal_method="orientation")
-
+ topomesh = up_facing_surface_topomesh(
+ topomesh,
+ upwards=self['orientation']=="up",
+ down_facing_threshold=self['down_threshold'],
+ normal_method='orientation',
+ connected=True
+ )
compute_topomesh_property(topomesh, 'vertices', 2)
points = topomesh.wisp_property('barycenter', 0).values(list(topomesh.wisps(0)))
@@ -120,9 +121,6 @@ class imageSurfaceMesh(gnomonAbstractMeshFromImage):
triangles = point_mapping.values(triangles)
topomesh = triangle_topomesh(triangles, points)
- #def cut_surface_topomesh(input_topomesh, z_cut=0, below=True):
- #nuclei_sigma=2., density_voxelsize=1., intensity_threshold=2000., microscope_orientation=1, maximal_length=10., remeshing_iterations=10,
- #erosion_radius=0.0, smoothing=50, decimation=100):
for degree in [1, 2, 3]:
compute_topomesh_property(topomesh, 'barycenter', degree)
diff --git a/src/gnomon_package_tissueimagemesh/algorithm/meshFromImage/segmentationSurfaceMesh.py b/src/gnomon_package_tissueimagemesh/algorithm/meshFromImage/segmentationSurfaceMesh.py
new file mode 100644
index 0000000000000000000000000000000000000000..3d4cd4fe43df392ff33cd635407638dc404e3dad
--- /dev/null
+++ b/src/gnomon_package_tissueimagemesh/algorithm/meshFromImage/segmentationSurfaceMesh.py
@@ -0,0 +1,69 @@
+import numpy as np
+
+from skimage.filters.thresholding import threshold_otsu
+
+from dtkcore import d_int, d_inliststring, 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.tissue_image_mesh import tissue_image_surface_topomesh
+
+from cellcomplex.property_topomesh.analysis import compute_topomesh_property
+
+
+@corePlugin(version="0.3.0", coreversion="0.81.0", name="Segmentation Surface Mesh")
+@cellImageInput("seg_img", data_plugin='gnomonCellImageDataTissueImage')
+@meshOutput("topomesh", data_plugin="gnomonMeshDataPropertyTopomesh")
+class segmentationSurfaceMesh(gnomonAbstractMeshFromImage):
+ """Reconstruct the surface on an segmented image as a triangular mesh.
+ """
+
+ 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['orientation'] = d_inliststring('Orientation', "up", ["up", "down"], "Whether to keep upper or lower part")
+
+ 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_surface_topomesh(
+ seg_img=seg_img,
+ resampling_voxelsize=resampling_voxelsize,
+ orientation=(1 if self['orientation']=="up" else -1),
+ surface_matching='cell',
+ padding=False,
+ )
+
+ 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
diff --git a/src/gnomon_package_tissueimagemesh/algorithm/pointCloudQuantification/nucleiCurvatureNukem.py b/src/gnomon_package_tissueimagemesh/algorithm/pointCloudQuantification/nucleiCurvatureNukem.py
index 1570e37eb7be202182ed2ddc11ff32c6ab00cfaa..00e6d39588ae635006537b849d73d47aa6ec3c12 100644
--- a/src/gnomon_package_tissueimagemesh/algorithm/pointCloudQuantification/nucleiCurvatureNukem.py
+++ b/src/gnomon_package_tissueimagemesh/algorithm/pointCloudQuantification/nucleiCurvatureNukem.py
@@ -2,15 +2,17 @@ from copy import deepcopy
import numpy as np
import pandas as pd
-from cellcomplex.property_topomesh.utils.pandas_tools import topomesh_to_dataframe
-from dtkcore import d_inliststring
-from dtkcore import d_int
+
+from dtkcore import d_inliststring, d_int
+
from gnomon.utils import algorithmPlugin
-from gnomon.utils.decorators import dataFrameOutput
-from gnomon.utils.decorators import imageInput
-from gnomon.utils.decorators import pointCloudInput
-from gnomon.utils.decorators import pointCloudOutput
+from gnomon.utils.decorators import imageInput, meshInput, pointCloudInput
+from gnomon.utils.decorators import dataFrameOutput, pointCloudOutput
+
from gnomon.core import gnomonAbstractPointCloudQuantification
+
+from cellcomplex.property_topomesh.utils.pandas_tools import topomesh_to_dataframe
+
from tissue_nukem_3d.nuclei_mesh_tools import nuclei_surface_topomesh_curvature
diff --git a/src/gnomon_package_tissueimagemesh/algorithm/pointCloudQuantification/nucleiLayerEstimation.py b/src/gnomon_package_tissueimagemesh/algorithm/pointCloudQuantification/nucleiLayerEstimation.py
index e3b851c3d589261895972fab3df129556336ff6e..76314ff0387585b3fe1c76c48f174be6cf66410c 100644
--- a/src/gnomon_package_tissueimagemesh/algorithm/pointCloudQuantification/nucleiLayerEstimation.py
+++ b/src/gnomon_package_tissueimagemesh/algorithm/pointCloudQuantification/nucleiLayerEstimation.py
@@ -6,7 +6,7 @@ import pandas as pd
from dtkcore import d_inliststring, d_int
from gnomon.utils import algorithmPlugin
-from gnomon.utils.decorators import pointCloudInput, imageInput, pointCloudOutput, dataFrameOutput
+from gnomon.utils.decorators import pointCloudInput, imageInput, meshInput, pointCloudOutput, dataFrameOutput
from gnomon.core import gnomonAbstractPointCloudQuantification
from timagetk.components.spatial_image import SpatialImage
@@ -15,6 +15,7 @@ from timagetk_geometry.features.nuclei import nuclei_layer
@algorithmPlugin(version="0.3.0", coreversion="0.81.0", name="Nuclei L1 Estimation")
@imageInput('img', data_plugin='gnomonImageDataMultiChannelImage')
+@meshInput('surface_topomesh', data_plugin="gnomonMeshDataPropertyTopomesh")
@pointCloudInput("df", data_plugin="gnomonPointCloudDataPandas")
@pointCloudOutput('out_df', data_plugin="gnomonPointCloudDataPandas")
@dataFrameOutput('data_df', data_plugin="gnomonDataFrameDataPandas")
@@ -31,6 +32,7 @@ class nucleiLayerEstimation(gnomonAbstractPointCloudQuantification):
super().__init__()
self.img = {}
+ self.surface_topomesh = {}
self.df = {}
self.out_df = {}
self.data_df = {}
@@ -71,6 +73,11 @@ class nucleiLayerEstimation(gnomonAbstractPointCloudQuantification):
nuclei_points = df[['center_'+dim for dim in 'xyz']].values
nuclei_positions = dict(zip(df['label'].values, nuclei_points))
+ if (time in self.surface_topomesh) and self.surface_topomesh[time] is not None:
+ surface_topomesh = self.surface_topomesh[time]
+ else:
+ surface_topomesh = None
+
if (time in self.img) and self.img[time] is not None:
img = self.img[time]
if 'surface_channel' in self._parameters.keys():
@@ -83,8 +90,10 @@ class nucleiLayerEstimation(gnomonAbstractPointCloudQuantification):
img = SpatialImage(np.zeros(shape), voxelsize=(1., 1., 1.))
surface_mode = "density"
+ print(surface_topomesh)
cell_layer, surface_topomesh = nuclei_layer(
nuclei_positions,
+ surface_topomesh=surface_topomesh,
nuclei_img=img,
upwards=(self['orientation'] == 'upright'),
surface_mode=surface_mode,
diff --git a/src/gnomon_package_tissueimagemesh/algorithm/pointCloudQuantification/surfaceMeshProperty.py b/src/gnomon_package_tissueimagemesh/algorithm/pointCloudQuantification/surfaceMeshProperty.py
new file mode 100644
index 0000000000000000000000000000000000000000..5a61becb6bbe3ae3d5e553de5557a694c53fbcfb
--- /dev/null
+++ b/src/gnomon_package_tissueimagemesh/algorithm/pointCloudQuantification/surfaceMeshProperty.py
@@ -0,0 +1,99 @@
+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 pointCloudInput, imageInput, meshInput, pointCloudOutput, dataFrameOutput
+from gnomon.core import gnomonAbstractPointCloudQuantification
+
+
+@algorithmPlugin(version="0.3.0", coreversion="0.81.0", name="Nuclei Mesh Property")
+@meshInput('surface_topomesh', data_plugin="gnomonMeshDataPropertyTopomesh")
+@pointCloudInput("df", data_plugin="gnomonPointCloudDataPandas")
+@pointCloudOutput('out_df', data_plugin="gnomonPointCloudDataPandas")
+@dataFrameOutput('data_df', data_plugin="gnomonDataFrameDataPandas")
+class surfaceMeshProperty(gnomonAbstractPointCloudQuantification):
+ """Transfer vertex attributes of a surface mesh to the nearest points.
+ """
+
+ def __init__(self):
+ super().__init__()
+
+ self.surface_topomesh = {}
+ self.df = {}
+ self.out_df = {}
+ self.data_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.data_df = {}
+ self.out_df = {}
+
+ for time in self.df.keys():
+ df = self.df[time]
+
+ surface_topomesh = self.surface_topomesh[time]
+
+ out_df = deepcopy(df)
+
+ cell_labels = df['label'].values
+ 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:
+ cell_centers = df[['center_'+dim for dim in 'xyz']].values
+
+ 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)
+ if property_is_binary:
+ out_df[property_name] = (nd.sum(vertex_property, vertex_cells, index=cell_labels) > 0).astype(int)
+ else:
+ if vertex_property.ndim == 1:
+ out_df[property_name] = nd.mean(vertex_property, vertex_cells, index=cell_labels)
+ elif vertex_property.ndim == 2:
+ out_df[property_name] = 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:
+ out_df[property_name] = 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)))
+
+ self.out_df[time] = out_df
+
+ self.data_df[time] = deepcopy(out_df)
+ self.data_df[time]['time'] = time
+
+ if len(self.data_df)>0:
+ self.data_df = {0: pd.concat([self.data_df[time] for time in self.data_df.keys()])}