add volume.isosurface_discrete()

docs/changes.md

@@ -58,6 +58,7 @@
 - improvements to `RayCastPlotter`
 - add `visual.scalar_range()` to control mesh coloring.
 - fix `shapes.Text3D.text()` by @gioda
+- add `volume.isosurface_discrete()` method
 
 
 ## Soft Breaking Changes
@@ -106,6 +107,8 @@ examples/volumetric/read_volume3.py
 examples/volumetric/rectl_grid1.py
 examples/volumetric/struc_grid1.py
 examples/volumetric/app_raycaster.py
+examples/volumetric/isosurfaces1.py
+examples/volumetric/isosurfaces2.py
 
 examples/simulations/mag_field1.py
 

docs/examples_db.js

@@ -828,7 +828,7 @@ vedo_example_db =
         imgsrc: 'images/volumetric/read_volume3.jpg',
     },
     {
-        pyname: 'isosurfaces',
+        pyname: 'isosurfaces1',
         kbd   : '',
         categ : 'volumetric',
         short : 'isosurface sets',

examples/volumetric/isosurfaces.py → examples/volumetric/isosurfaces1.py

@@ -1,11 +1,9 @@
 """Interactively cut a set of isosurfaces from a volumetric dataset"""
 from vedo import dataurl, show, BoxCutter, Volume
 
-# generate an isosurface the volume for each thresholds
-thresholds = [0.1, 0.25, 0.4, 0.6, 0.75, 0.9]
-
-# isos is of type Mesh
-isos = Volume(dataurl+'quadric.vti').isosurface(thresholds)
+# generate an isosurface the volume for each value
+values = [0.1, 0.25, 0.4, 0.6, 0.75, 0.9]
+isos = Volume(dataurl+'quadric.vti').isosurface(values) # Mesh
 
 plt = show(isos, __doc__, axes=1, interactive=False)
 

examples/volumetric/isosurfaces2.py

@@ -0,0 +1,29 @@
+"""Isosurface extraction from a volume dataset with discrete values."""
+from vedo import *
+
+lut = build_lut(
+    [
+        [0, "lightyellow"],
+        [1, "red"],
+        [2, "blue"],
+        [3, "yellow"],
+        [4, "orange"],
+        [5, "cyan"],
+        [6, "magenta"],
+        [7, "white"],
+        [8, "pink"],
+        [9, "brown"],
+        [10, "lightblue"],
+        [11, "lightgreen"],
+    ],
+    interpolate=False
+)
+
+# this dataset is a 3D volume of 64x64x64 voxels containing 12 "blobs"
+blobs = Volume(dataurl + "blobs.vti")
+
+isovalues = list(range(12))
+iso_discrete = blobs.isosurface_discrete(isovalues, nsmooth=10)
+iso_discrete.cmap(lut)
+
+show(iso_discrete, __doc__, axes=1).close()

vedo/core.py

@@ -1805,7 +1805,7 @@ def isosurface(self, value=None, flying_edges=False):
         Use flying_edges for faster results (but sometimes can interfere with `smooth()`).
 
         Examples:
-            - [isosurfaces.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/isosurfaces.py)
+            - [isosurfaces1.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/isosurfaces1.py)
 
                 ![](https://vedo.embl.es/images/volumetric/isosurfaces.png)
         """
@@ -1844,6 +1844,57 @@ def isosurface(self, value=None, flying_edges=False):
             c="#4cc9f0:#e9c46a",
         )
         return out
+
+    def isosurface_discrete(self, value=None, nsmooth=15):
+        """
+        Create boundary/isocontour surfaces from a label map (e.g., a segmented image) using a threaded,
+        3D version of the multiple objects/labels Surface Nets algorithm.
+        The input is a 3D image (i.e., volume) where each voxel is labeled
+        (integer labels are preferred to real values), and the output data is a polygonal mesh separating
+        labeled regions / objects.
+        (Note that on output each region [corresponding to a different segmented object] will share
+        points/edges on a common boundary, i.e., two neighboring objects will share the boundary that separates them).
+
+        Arguments:
+            value : (float, list)
+                single value or list of values to draw the isosurface(s).
+            nsmooth : (int)
+                number of iterations of smoothing (0 means no smoothing).
+
+        Examples:
+            - [isosurfaces2.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/isosurfaces2.py)
+        """
+        if not utils.is_sequence(value):
+            value = [value]
+
+        snets = vtki.new("SurfaceNets3D")
+        snets.SetInputData(self.dataset)
+
+        if nsmooth:
+            snets.SmoothingOn()
+            snets.AutomaticSmoothingConstraintsOn()
+            snets.GetSmoother().SetNumberOfIterations(nsmooth)
+            # snets.GetSmoother().SetRelaxationFactor(relaxation_factor)
+            # snets.GetSmoother().SetConstraintDistance(constraint_distance)
+        else:
+            snets.SmoothingOff()
+
+        for i, val in enumerate(value):
+            snets.SetValue(i, val)
+        snets.Update()
+        snets.SetOutputMeshTypeToTriangles()
+        snets.SetOutputStyleToBoundary()
+        snets.Update()
+
+        out = vedo.mesh.Mesh(snets.GetOutput())
+        out.pipeline = utils.OperationNode(
+            "isosurface_discrete",
+            parents=[self],
+            comment=f"#pts {out.dataset.GetNumberOfPoints()}",
+            c="#4cc9f0:#e9c46a",
+        )
+        return out
+
 
     def legosurface(
         self,

vedo/grids.py

@@ -608,7 +608,7 @@ def isosurface(self, value=None, flying_edges=False):
         Use flying_edges for faster results (but sometimes can interfere with `smooth()`).
 
         Examples:
-            - [isosurfaces.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/isosurfaces.py)
+            - [isosurfaces1.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/isosurfaces1.py)
 
                 ![](https://vedo.embl.es/images/volumetric/isosurfaces.png)
         """

vedo/vtkclasses.py

@@ -270,6 +270,7 @@
     "vtkResampleWithDataSet",
     "vtkReverseSense",
     "vtkStripper",
+    "vtkSurfaceNets3D",
     "vtkTensorGlyph",
     "vtkThreshold",
     "vtkTriangleFilter",