Merge pull request #104 from vtkiorg/split-bodies

Add new split_bodies filter

docs/examples/quick/compute_volume.rst

@@ -86,3 +86,48 @@ cell sizes, then extract the volumes of each body:
     Low grade volume: 518.0
     High grade volume: 35.0
     Original volume: 729.0
+
+
+
+Splitting Volumes
+=================
+
+What if instead, we wanted to split all the different connected bodies/volumes
+in a dataset like the one above? We could use the
+:func:`vtki.DataSetFilters.split_bodies` filter to extract all the different
+connected volumes in a dataset into blocks in a :class:`vtki.MultiBlock`
+dataset. For example, lets split the thresholded volume in the example above:
+
+
+.. testcode:: python
+
+    import numpy as np
+    import vtki
+    from vtki import examples
+    vtki.set_plot_theme('document')
+
+    # Load a simple example mesh
+    dataset = examples.load_uniform()
+    dataset.set_active_scalar('Spatial Cell Data')
+    threshed = dataset.threshold_percent([0.15, 0.50], invert=True)
+
+    bodies = threshed.split_bodies()
+
+    for i, body in enumerate(bodies):
+        print('Body {} volume: {:.3f}'.format(i, body.volume))
+
+
+.. testoutput:: python
+   :hide:
+   :options: -ELLIPSIS, +NORMALIZE_WHITESPACE
+
+    Body 0 volume: 518.000
+    Body 1 volume: 35.000
+
+
+.. code-block:: python
+
+    bodies.plot(show_bounds=True, multi_colors=True)
+
+
+.. image:: ../../images/split-bodies.png

tests/test_filters.py

@@ -247,3 +247,17 @@ def test_glyph():
     sphere.point_arrays['arr'] = np.ones(sphere.n_points)
     result = sphere.glyph(scale='arr')
     result = sphere.glyph(scale='arr', orient='Normals', factor=0.1)
+
+
+def test_split_and_connectivity():
+    # Load a simple example mesh
+    dataset = examples.load_uniform()
+    dataset.set_active_scalar('Spatial Cell Data')
+    threshed = dataset.threshold_percent([0.15, 0.50], invert=True)
+
+    bodies = threshed.split_bodies()
+
+    volumes = [518.0, 35.0]
+    assert len(volumes) == bodies.n_blocks
+    for i, body in enumerate(bodies):
+        assert np.allclose(body.volume, volumes[i], rtol=0.1)

vtki/filters.py

@@ -666,7 +666,7 @@ def compute_cell_sizes(dataset, length=False, area=True, volume=True):
         alg.Update()
         return _get_output(alg)
 
-    def cell_centers(self, vertex=True):
+    def cell_centers(dataset, vertex=True):
         """Generate points at the center of the cells in this dataset.
         These points can be used for placing glyphs / vectors.
 
@@ -676,14 +676,14 @@ def cell_centers(self, vertex=True):
             Enable/disable the generation of vertex cells.
         """
         alg = vtk.vtkCellCenters()
-        alg.SetInputDataObject(self)
+        alg.SetInputDataObject(dataset)
         alg.SetVertexCells(vertex)
         alg.Update()
         output = _get_output(alg)
         return output
 
 
-    def glyph(self, orient=True, scale=True, factor=1.0, geom=None):
+    def glyph(dataset, orient=True, scale=True, factor=1.0, geom=None):
         """
         Copies a geometric representation (called a glyph) to every
         point in the input dataset.  The glyph may be oriented along
@@ -711,16 +711,57 @@ def glyph(self, orient=True, scale=True, factor=1.0, geom=None):
         alg = vtk.vtkGlyph3D()
         alg.SetSourceData(geom)
         if isinstance(scale, str):
-            self.active_scalar_name = scale
+            dataset.active_scalar_name = scale
             scale = True
         if scale:
             alg.SetScaleModeToScaleByScalar()
         if isinstance(orient, str):
-            self.active_vectors_name = orient
+            dataset.active_vectors_name = orient
             orient = True
         alg.SetOrient(orient)
-        alg.SetInputData(self)
+        alg.SetInputData(dataset)
         alg.SetVectorModeToUseVector()
         alg.SetScaleFactor(factor)
         alg.Update()
         return _get_output(alg)
+
+
+    def connectivity(dataset):
+        """Find and label connected bodies/volumes. This adds an ID array to
+        the point and cell data to distinguish seperate connected bodies.
+        This applies a ``vtkConnectivityFilter`` filter which extracts cells
+        that share common points and/or meet other connectivity criterion.
+        (Cells that share vertices and meet other connectivity criterion such
+        as scalar range are known as a region.)
+        """
+        alg = vtk.vtkConnectivityFilter()
+        alg.SetInputData(dataset)
+        alg.SetExtractionModeToAllRegions()
+        alg.SetColorRegions(True)
+        alg.Update()
+        return _get_output(alg)
+
+
+    def split_bodies(dataset, label=False):
+        """Find, label, and split connected bodies/volumes. This splits
+        different connected bodies into blocks in a MultiBlock dataset.
+
+        Parameters
+        ----------
+        label : bool
+            A flag on whether to keep the ID arrays given by the
+            ``connectivity`` filter.
+        """
+        # Get the connectivity and label different bodies
+        labeled = dataset.connectivity()
+        classifier = labeled.cell_arrays['RegionId']
+        bodies = vtki.MultiBlock()
+        for vid in np.unique(classifier):
+            # Now extract it:
+            b = labeled.threshold([vid-0.5, vid+0.5], scalars='RegionId')
+            if not label:
+                del b.cell_arrays['RegionId']
+                del b.point_arrays['RegionId']
+            bodies.append(b)
+
+        return bodies