PineRootConnectivity

VTKExamples/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.

Code

PineRootConnectivity.cxx

#include <vtkActor.h>
#include <vtkCamera.h>
#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>

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;
    vtkIdType* pts;
    auto numOfTriangles = 0;
    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++;
      }
    }
    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 vtkPolyData::SafeDownCast here.
    std::cout << "There are: " << NumberofTriangles(vtkPolyData::SafeDownCast(
                                    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 2.8)

PROJECT(PineRootConnectivity)

find_package(VTK REQUIRED)
include(${VTK_USE_FILE})

add_executable(PineRootConnectivity MACOSX_BUNDLE PineRootConnectivity.cxx )

target_link_libraries(PineRootConnectivity ${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.