Add label smoothing (morphometrics#54)

* add label smoothing

* update docstring

* fix numpy dtype

* fix pandas iteritems deprecation

* remove repair mesh flag

src/morphometrics/_gui/models/annotation_model.py

@@ -269,7 +269,7 @@ def start_annotation(
 
     def _update_table(self, df: pd.DataFrame, new_column: pd.Series):
         column_name = new_column.name
-        for index, value in new_column.iteritems():
+        for index, value in new_column.items():
             df.at[index, column_name] = value
 
     def _get_group_by_keys(self, combo_widget=None) -> List[str]:

src/morphometrics/label/_tests/test_image_utils.py

@@ -161,7 +161,7 @@ def test_get_mask_bounding_box_3d():
 
     bounding_box = get_mask_bounding_box_3d(mask_image)
     np.testing.assert_allclose(bounding_box, expected_bounding_box)
-    assert bounding_box.dtype == np.int
+    assert bounding_box.dtype == int
 
 
 def test_expand_bounding_box():

src/morphometrics/measure/label.py

@@ -212,7 +212,7 @@ def measure_surface_properties_from_labels(
             continue
         label_mask = label_image == label_index
         smoothed_mesh = binary_mask_to_surface(
-            label_mask, n_mesh_smoothing_interations=n_mesh_smoothing_interations
+            label_mask, n_mesh_smoothing_iterations=n_mesh_smoothing_interations
         )
         object_measurements = measure_surface_properties(surface=smoothed_mesh)
 

src/morphometrics/utils/_tests/test_surface_utils.py

@@ -49,9 +49,7 @@ def test_voxelize_closed_surface():
     pitch = 0.5
     mesh = _make_cuboid_mesh(origin=origin, extents=cube_extents)
 
-    voxelized, image_origin = voxelize_closed_surface(
-        mesh, pitch=pitch, repair_mesh=True
-    )
+    voxelized, image_origin = voxelize_closed_surface(mesh, pitch=pitch)
 
     np.testing.assert_allclose(image_origin, [9.5, -0.5, 9.5])
     np.testing.assert_allclose([63, 63, 63], voxelized.shape)
@@ -67,7 +65,6 @@ def test_closed_surfaces_to_label_image_no_crop():
         [mesh_0, mesh_1],
         pitch=pitch,
         crop_around_mesh=False,
-        repair_mesh=True,
     )
 
     np.testing.assert_allclose(image_origin, [0, 0, 0])

src/morphometrics/utils/surface_utils.py

@@ -1,11 +1,11 @@
-from typing import List, Tuple
+from typing import List, Optional, Tuple
 
 import numpy as np
 import pymeshfix
 import trimesh.voxel.creation
 from skimage.measure import marching_cubes
 from trimesh import Trimesh
-from trimesh.smoothing import filter_taubin
+from trimesh.smoothing import filter_mut_dif_laplacian
 
 from ..types import BinaryImage, LabelImage
 
@@ -51,19 +51,25 @@ def repair_mesh(mesh: Trimesh) -> Trimesh:
 
 
 def binary_mask_to_surface(
-    object_mask: BinaryImage, n_mesh_smoothing_interations: int = 50
+    object_mask: BinaryImage,
+    n_mesh_smoothing_iterations: int = 10,
+    diffusion_coefficient: float = 0.5,
 ) -> Trimesh:
     """Convert surface of a 3D binary mask (segmented object) into a watertight mesh.
 
     Parameters
     ----------
     object_mask  : BinaryMask
         A 3D binary image corresponding to the object you want to mesh.
-    n_mesh_smoothing_interations : int
-        The number of interations of smooting to perform. Smoothing is
-        done by the trimesh taubin filter:
-        https://trimsh.org/trimesh.smoothing.html#trimesh.smoothing.filter_taubin
-        Default value is 50.
+    n_mesh_smoothing_iterations : int
+        The number of interations of smoothing to perform. Smoothing is
+        done by the trimesh mutable diffusion laplacian filter:
+        https://trimsh.org/trimesh.smoothing.html#trimesh.smoothing.filter_mut_dif_laplacian
+        Default value is 10.
+    diffusion_coefficient : float
+        The diffusion coefficient for smoothing. 0 is no diffusion.
+        Default value is 0.5.
+        https://trimsh.org/trimesh.smoothing.html#trimesh.smoothing.filter_mut_dif_laplacian
 
     Returns
     -------
@@ -79,14 +85,71 @@ def binary_mask_to_surface(
     mesh = Trimesh(vertices=vertices_clean, faces=faces_clean)
 
     # optionally clean up the mesh
-    if n_mesh_smoothing_interations > 0:
-        filter_taubin(mesh, iterations=n_mesh_smoothing_interations)
+    if n_mesh_smoothing_iterations > 0:
+        filter_mut_dif_laplacian(
+            mesh, iterations=n_mesh_smoothing_iterations, lamb=diffusion_coefficient
+        )
 
     return mesh
 
 
+def label_image_to_surface_collection(
+    label_image: LabelImage,
+    label_smoothing_radius: Optional[int] = None,
+    n_mesh_smoothing_iterations: int = 10,
+    diffusion_coefficient: float = 0.5,
+) -> List[Trimesh]:
+    """Convert a label image into a collection of surfaces.
+
+    Parameters
+    ----------
+    label_image : LabelImage
+        The image to transform into surfaces.
+    label_smoothing_radius : Optional[int]
+        The radius of the morphological filter to smooth the label
+        image before meshing. This is similar to performing a morphological
+        opening on the label image.
+        If None, no label smoothing is performed. Default value is None.
+     n_mesh_smoothing_iterations : int
+        The number of interations of smoothing to perform. Smoothing is
+        done by the trimesh mutable diffusion laplacian filter:
+        https://trimsh.org/trimesh.smoothing.html#trimesh.smoothing.filter_mut_dif_laplacian
+        Default value is 10.
+    diffusion_coefficient : float
+        The diffusion coefficient for smoothing. 0 is no diffusion.
+        Default value is 0.5.
+        https://trimsh.org/trimesh.smoothing.html#trimesh.smoothing.filter_mut_dif_laplacian
+
+    Returns
+    -------
+    mesh : trimesh.Trimesh
+        The resulting mesh as a trimesh.Trimesh object.
+        https://trimsh.org/trimesh.base.html#github-com-mikedh-trimesh
+    """
+    import pyclesperanto_prototype as cle
+
+    if label_smoothing_radius is not None:
+        # perform label smoothing
+        label_image = np.asarray(
+            cle.smooth_labels(label_image, radius=label_smoothing_radius)
+        )
+
+    label_values = np.unique(label_image)
+
+    return [
+        binary_mask_to_surface(
+            label_image == label_index,
+            n_mesh_smoothing_iterations=n_mesh_smoothing_iterations,
+            diffusion_coefficient=diffusion_coefficient,
+        )
+        for label_index in label_values
+        if label_index != 0
+    ]
+
+
 def voxelize_closed_surface(
-    mesh: Trimesh, pitch: float, repair_mesh: bool = True
+    mesh: Trimesh,
+    pitch: float,
 ) -> Tuple[BinaryImage, np.ndarray]:
     """Voxelize a closed surface mesh.
 
@@ -97,9 +160,6 @@ def voxelize_closed_surface(
     pitch : float
         The voxel width in mesh units. Voxels have the
         same width in each dimension (i.e., are cubes).
-    repair_mesh : bool
-        Flag to attept to repair the mesh if set to True.
-        Default value is True.
 
     Returns
     -------
@@ -115,7 +175,7 @@ def voxelize_closed_surface(
     # convert the centroid to the nearest integer multiple of the pitch
     rounded_centroid = _round_to_pitch(coordinate=centroid, pitch=pitch)
 
-    # find the minimum cube half-width that encompases the full mesh
+    # find the minimum cube half-width that encompasses the full mesh
     cube_half_width = np.max(bounding_box - rounded_centroid)
 
     # get the number of voxels for the cube half-width
@@ -208,9 +268,7 @@ def closed_surfaces_to_label_image(
     label_image = np.zeros(image_shape_voxels, dtype=np.uint16)
 
     for i, mesh in enumerate(meshes):
-        voxelized, origin = voxelize_closed_surface(
-            mesh, pitch=pitch, repair_mesh=repair_mesh
-        )
+        voxelized, origin = voxelize_closed_surface(mesh, pitch=pitch)
 
         # get the coordinates of the voxels inside of the mesh
         filled_voxel_coordinates = np.argwhere(voxelized)