Skip to content


# PineRootConnectivity

VTKEx/Cxx/VisualizationAlgorithms/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 (Python)

Question

If you have a simple question about this example contact us at VTKExProject If your question is more complex and may require extended discussion, please use the VTK Discourse Forum

### Code

PineRootConnectivity.cxx

    #include <vtkActor.h>
    #include <vtkCamera.h>
    #include <vtkCellArray.h>
    #ifdef VTK_CELL_ARRAY_V2
    #include <vtkCellArrayIterator.h>
    #endif // VTK_CELL_ARRAY_V2
    #include <vtkDataSetMapper.h>
    #include <vtkMCubesReader.h>
    #include <vtkNamedColors.h>
    #include <vtkOutlineFilter.h>
    #include <vtkPolyData.h>
    #include <vtkPolyDataConnectivityFilter.h>
    #include <vtkPolyDataMapper.h>
    #include <vtkProperty.h>
    #include <vtkRenderWindow.h>
    #include <vtkRenderWindowInteractor.h>
    #include <vtkRenderer.h>
    #include <vtkSmartPointer.h>

    int main(int argc, char* argv[])
    {
      if (argc < 2)
      {
        std::cout << "Usage: " << argv[0] << " filename [noConnectivity]"
                  << std::endl;
        std::cout << "where: filename is pine_root.tri." << std::endl;
        std::cout << "       if noConnectivity is nonzero then the connectivity "
                     "filter is ignored."
                  << std::endl;
        return EXIT_FAILURE;
      }

      /**
       * Count the triangles in the polydata.
       * @param pd: vtkPolyData.
       * @return The number of triangles.
       */
      auto NumberofTriangles = [](vtkPolyData* pd) {
        vtkCellArray* cells = pd->GetPolys();
        vtkIdType npts = 0;
        auto numOfTriangles = 0;
    #ifdef VTK_CELL_ARRAY_V2

        // Newer versions of vtkCellArray prefer local iterators:
        std::cout << pd->GetNumberOfPolys() << std::endl;
        auto cellIter = vtk::TakeSmartPointer(cells->NewIterator());
        for (cellIter->GoToFirstCell(); !cellIter->IsDoneWithTraversal();
             cellIter->GoToNextCell())
        {
          auto cell = cellIter->GetCurrentCell();
          if (cell == nullptr)
          {
            break;
          }
          // If a cell has three points it is a triangle.
          if (cell->GetNumberOfIds() == 3)
          {
            numOfTriangles++;
          }
        }

    #else  // VTK_CELL_ARRAY_V2

        // Older implementations of vtkCellArray use internal iterator APIs (not
        // thread safe):
        vtkIdType* pts;
        for (auto i = 0; i < pd->GetNumberOfPolys(); ++i)
        {
          int c = cells->GetNextCell(npts, pts);
          if (c == 0)
          {
            break;
          }
          // If a cell has three points it is a triangle.
          if (npts == 3)
          {
            numOfTriangles++;
          }
        }
    #endif // VTK_CELL_ARRAY_V2

        return numOfTriangles;
      };

      std::string fileName = argv[1];
      auto noConnectivity = false;
      if (argc > 2)
      {
        noConnectivity = atoi(argv[2]) != 0;
      }

      vtkSmartPointer<vtkNamedColors> colors =
          vtkSmartPointer<vtkNamedColors>::New();

      // Create the pipeline.
      vtkSmartPointer<vtkMCubesReader> reader =
          vtkSmartPointer<vtkMCubesReader>::New();
      reader->SetFileName(fileName.c_str());
      reader->FlipNormalsOff();
      if (!noConnectivity)
      {
        reader->Update();
        std::cout << "Before Connectivity." << std::endl;
        std::cout << "There are: " << NumberofTriangles(reader->GetOutput())
                  << " triangles." << std::endl;
      }

      vtkSmartPointer<vtkPolyDataConnectivityFilter> connect =
          vtkSmartPointer<vtkPolyDataConnectivityFilter>::New();
      connect->SetInputConnection(reader->GetOutputPort());
      connect->SetExtractionModeToLargestRegion();
      if (!noConnectivity)
      {
        connect->Update();
        std::cout << "After Connectivity." << std::endl;
        // Note the use of dynamic_cast<vtkPolyData*> here.
        std::cout << "There are: "
                  << NumberofTriangles(
                         dynamic_cast<vtkPolyData*>(connect->GetOutput()))
                  << " triangles." << std::endl;
      }

      vtkSmartPointer<vtkDataSetMapper> isoMapper =
          vtkSmartPointer<vtkDataSetMapper>::New();
      if (noConnectivity)
      {
        isoMapper->SetInputConnection(reader->GetOutputPort());
      }
      else
      {
        isoMapper->SetInputConnection(connect->GetOutputPort());
      }
      isoMapper->ScalarVisibilityOff();
      vtkSmartPointer<vtkActor> isoActor = vtkSmartPointer<vtkActor>::New();
      isoActor->SetMapper(isoMapper);
      isoActor->GetProperty()->SetColor(colors->GetColor3d("raw_sienna").GetData());

      // Get an outline of the data set for context.
      vtkSmartPointer<vtkOutlineFilter> outline =
          vtkSmartPointer<vtkOutlineFilter>::New();
      outline->SetInputConnection(reader->GetOutputPort());
      vtkSmartPointer<vtkPolyDataMapper> outlineMapper =
          vtkSmartPointer<vtkPolyDataMapper>::New();
      outlineMapper->SetInputConnection(outline->GetOutputPort());
      vtkSmartPointer<vtkActor> outlineActor = vtkSmartPointer<vtkActor>::New();
      outlineActor->SetMapper(outlineMapper);
      outlineActor->GetProperty()->SetColor(colors->GetColor3d("Black").GetData());

      // Create the Renderer, RenderWindow and RenderWindowInteractor.
      vtkSmartPointer<vtkRenderer> ren = vtkSmartPointer<vtkRenderer>::New();
      vtkSmartPointer<vtkRenderWindow> renWin =
          vtkSmartPointer<vtkRenderWindow>::New();
      renWin->AddRenderer(ren);
      vtkSmartPointer<vtkRenderWindowInteractor> iren =
          vtkSmartPointer<vtkRenderWindowInteractor>::New();
      iren->SetRenderWindow(renWin);

      // Add the actors to the renderer, set the background and size.
      ren->AddActor(outlineActor);
      ren->AddActor(isoActor);
      // renWin->SetSize(750, 750);
      renWin->SetSize(512, 512);
      ren->SetBackground(colors->GetColor3d("SlateGray").GetData());

      vtkCamera* 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();

      return EXIT_SUCCESS;
    }

### CMakeLists.txt

    cmake_minimum_required(VERSION 3.3 FATAL_ERROR)

    project(PineRootConnectivity)

    find_package(VTK COMPONENTS
      vtkCommonColor
      vtkCommonCore
      vtkCommonDataModel
      vtkFiltersCore
      vtkFiltersModeling
      vtkIOGeometry
      vtkInteractionStyle
      vtkRenderingContextOpenGL2
      vtkRenderingCore
      vtkRenderingFreeType
      vtkRenderingGL2PSOpenGL2
      vtkRenderingOpenGL2 QUIET)
    if (NOT VTK_FOUND)
      message("Skipping PineRootConnectivity: ${VTK_NOT_FOUND_MESSAGE}")
      return ()
    endif()
    message (STATUS "VTK_VERSION: ${VTK_VERSION}")
    if (VTK_VERSION VERSION_LESS "8.90.0")
      # old system
      include(${VTK_USE_FILE})
      add_executable(PineRootConnectivity MACOSX_BUNDLE PineRootConnectivity.cxx )
      target_link_libraries(PineRootConnectivity PRIVATE ${VTK_LIBRARIES})
    else ()
      # include all components
      add_executable(PineRootConnectivity MACOSX_BUNDLE PineRootConnectivity.cxx )
      target_link_libraries(PineRootConnectivity PRIVATE ${VTK_LIBRARIES})
      # vtk_module_autoinit is needed
      vtk_module_autoinit(
        TARGETS PineRootConnectivity
        MODULES ${VTK_LIBRARIES}
        )
    endif ()

### Download and Build PineRootConnectivity

Click here to download PineRootConnectivity and its CMakeLists.txt file. Once the tarball PineRootConnectivity.tar has been downloaded and extracted,

    cd PineRootConnectivity/build

If VTK is installed:

    cmake ..

If VTK is not installed but compiled on your system, you will need to specify the path to your VTK build:

    cmake -DVTK_DIR:PATH=/home/me/vtk_build ..

Build the project:

    make

and run it:

    ./PineRootConnectivity

WINDOWS USERS

Be sure to add the VTK bin directory to your path. This will resolve the VTK dll's at run time.