Skip to content

PineRootConnectivity

Repository source: 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), (PythonicAPI)

Question

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

Code

PineRootConnectivity.cxx

#include <vtkActor.h>
#include <vtkCamera.h>
#include <vtkCellArray.h>
#include <vtkDataSetMapper.h>
#include <vtkMCubesReader.h>
#include <vtkNamedColors.h>
#include <vtkNew.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>

#ifdef VTK_CELL_ARRAY_V2
#include <vtkCellArrayIterator.h>
#endif // VTK_CELL_ARRAY_V2

#include <iostream>
#include <string>

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:
    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;
  }

  vtkNew<vtkNamedColors> colors;

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

  vtkNew<vtkPolyDataConnectivityFilter> connect;
  connect->SetInputConnection(reader->GetOutputPort());
  connect->SetExtractionModeToLargestRegion();
  if (!noConnectivity)
  {
    connect->Update();
    std::cout << "After Connectivity there are:  "
              << NumberofTriangles(
                     dynamic_cast<vtkPolyData*>(connect->GetOutput()))
              << " triangles." << std::endl;
  }

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

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

  // Create the Renderer, RenderWindow and RenderWindowInteractor.
  vtkNew<vtkRenderer> ren;
  vtkNew<vtkRenderWindow> renWin;
  renWin->AddRenderer(ren);
  vtkNew<vtkRenderWindowInteractor> iren;
  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);
  renWin->SetWindowName("PineRootConnectivity");

  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.12 FATAL_ERROR)

project(PineRootConnectivity)

find_package(VTK COMPONENTS 
  CommonColor
  CommonCore
  CommonDataModel
  FiltersCore
  FiltersModeling
  IOGeometry
  InteractionStyle
  RenderingContextOpenGL2
  RenderingCore
  RenderingFreeType
  RenderingGL2PSOpenGL2
  RenderingOpenGL2
)

if (NOT VTK_FOUND)
  message(FATAL_ERROR "PineRootConnectivity: Unable to find the VTK build folder.")
endif()

# Prevent a "command line is too long" failure in Windows.
set(CMAKE_NINJA_FORCE_RESPONSE_FILE "ON" CACHE BOOL "Force Ninja to use response files.")
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}
)

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.