PineRootConnectivity

Description¶

Demonstrates how to apply the connectivity filter to remove noisy isosurfaces.

To illustrate the application of connectivity analysis, we will use an MRI dataset generated by Janet MacFall at the Center for In Vivo Microscopy at Duke University. The dataset is a volume of dimensions 256^3. The data is of the root system of a small pine tree. Using the class vtkSliceCubes , an implementation of marching cubes for large volumes, we generate an initial isosurface represented by 351,536 triangles. This is a faster way of manipulating this data. If you have a large enough computer you can process the volume directly with a vtk image reader and vtkMarchingCubes. The example on this other page shows the initial dataset. Notice that there are many small, disconnected isosurfaces due to noise and isolated moisture in the data. We use vtkConnectivityFilter to remove these small, disconnected surfaces. The example on this page shows the result of applying the filter. Over 50,000 triangles were removed, leaving 299,380 triangles. The vtkConnectivityFilter is a general filter taking datasets as input, and generating an unstructured grid as output. It functions by extracting cells that are connected at points (i.e., share common points). In this example the single largest surface is extracted. It is also possible to specify cell ids and point ids and extract surfaces connected to these.

Info

To count the number of triangles in the polydata we do the following:

C++¶

We use a lambda function c++ auto NumberofTriangles = [](auto *pd)

Python¶

We just implement: python def NumberOfTriangles(pd): within the scope of the calling function.

Cite

Comparative Water Uptake by Roots of Different Ages in Seedlings of Loblolly Pine (Pinus taeda L.) December 1991. New Phytologist 119(4):551 - 560.

Info

See Figure 9-51 in Chapter 9 in the VTK Textbook.

Other languages

See (Cxx), (PythonicAPI)

Question

If you have a question about this example, please use the VTK Discourse Forum

Code¶

PineRootConnectivity.py

#!/usr/bin/env python

# noinspection PyUnresolvedReferences
import vtkmodules.vtkInteractionStyle
# noinspection PyUnresolvedReferences
import vtkmodules.vtkRenderingOpenGL2
from vtkmodules.vtkCommonColor import vtkNamedColors
from vtkmodules.vtkCommonCore import vtkIdList
from vtkmodules.vtkFiltersCore import (
    vtkPolyDataConnectivityFilter
)
from vtkmodules.vtkFiltersModeling import vtkOutlineFilter
from vtkmodules.vtkIOGeometry import vtkMCubesReader
from vtkmodules.vtkRenderingCore import (
    vtkActor,
    vtkPolyDataMapper,
    vtkRenderWindow,
    vtkRenderWindowInteractor,
    vtkRenderer
)


def pine_root_connectivity(fileName, noConnectivity):
    def NumberOfTriangles(pd):
        """
        Count the number of triangles.
        :param pd: vtkPolyData.
        :return: The number of triangles.
        """
        cells = pd.GetPolys()
        numOfTriangles = 0
        idList = vtkIdList()
        for i in range(0, cells.GetNumberOfCells()):
            cells.GetNextCell(idList)
            # If a cell has three points it is a triangle.
            if idList.GetNumberOfIds() == 3:
                numOfTriangles += 1
        return numOfTriangles

    colors = vtkNamedColors()

    # Create the pipeline.
    reader = vtkMCubesReader()
    reader.SetFileName(fileName)
    if not noConnectivity:
        reader.Update()
        print('Before Connectivity.')
        print('There are: ', NumberOfTriangles(reader.GetOutput()), 'triangles')

    connect = vtkPolyDataConnectivityFilter()
    connect.SetInputConnection(reader.GetOutputPort())
    connect.SetExtractionModeToLargestRegion()
    if not noConnectivity:
        connect.Update()
        print('After Connectivity.')
        print('There are: ', NumberOfTriangles(connect.GetOutput()), 'triangles')

    isoMapper = vtkPolyDataMapper()
    if noConnectivity:
        isoMapper.SetInputConnection(reader.GetOutputPort())
    else:
        isoMapper.SetInputConnection(connect.GetOutputPort())
    isoMapper.ScalarVisibilityOff()
    isoActor = vtkActor()
    isoActor.SetMapper(isoMapper)
    isoActor.GetProperty().SetColor(colors.GetColor3d('raw_sienna'))

    #  Get an outline of the data set for context.
    outline = vtkOutlineFilter()
    outline.SetInputConnection(reader.GetOutputPort())
    outlineMapper = vtkPolyDataMapper()
    outlineMapper.SetInputConnection(outline.GetOutputPort())
    outlineActor = vtkActor()
    outlineActor.SetMapper(outlineMapper)
    outlineActor.GetProperty().SetColor(colors.GetColor3d('Black'))

    # Create the Renderer, RenderWindow and RenderWindowInteractor.
    ren = vtkRenderer()
    renWin = vtkRenderWindow()
    renWin.AddRenderer(ren)
    iren = vtkRenderWindowInteractor()
    iren.SetRenderWindow(renWin)

    # Add the actors to the renderer, set the background and size.
    ren.AddActor(outlineActor)
    ren.AddActor(isoActor)
    renWin.SetSize(512, 512)
    renWin.SetWindowName('PineRootConnectivity')

    ren.SetBackground(colors.GetColor3d('SlateGray'))

    # render the image
    #
    # iren AddObserver UserEvent {wm deiconify .vtkInteract}
    cam = ren.GetActiveCamera()
    cam.SetFocalPoint(40.6018, 37.2813, 50.1953)
    cam.SetPosition(40.6018, -280.533, 47.0172)
    cam.ComputeViewPlaneNormal()
    cam.SetClippingRange(26.1073, 1305.36)
    cam.SetViewAngle(20.9219)
    cam.SetViewUp(0.0, 0.0, 1.0)

    iren.Initialize()
    renWin.Render()
    iren.Start()


def main():
    fileName, noConnectivity = get_program_parameters()
    noConnectivity = noConnectivity != 0
    pine_root_connectivity(fileName, noConnectivity)


def get_program_parameters():
    import argparse
    description = 'Applying connectivity filter to remove noisy isosurfaces.'
    epilogue = '''
        Applying connectivity filter to remove noisy isosurfaces.

This example demonstrates how to use the vtkConnectivityFilter.
If the extra parameter 'noConnectivity' is non zero, the connectivity filter will not be used.

   '''
    parser = argparse.ArgumentParser(description=description, epilog=epilogue,
                                     formatter_class=argparse.RawDescriptionHelpFormatter)
    parser.add_argument('filename', help='pine_root.tri.')
    parser.add_argument('noConnectivity', default=0, type=int, nargs='?',
                        help='If non-zero do not use the connectivity filter.')
    args = parser.parse_args()
    return args.filename, args.noConnectivity


if __name__ == '__main__':
    main()