VoxelsOnBoundary

Description¶

This example uses vtkMultiThreshold to extract voxels that are inside an isosurface and on the boundary of the isosurface. The result is a vtkUnstructuredGrid for each set of voxels. Before processing, vtkImageShrink3D reduces the resolution by a factor of 4.

Compare these results with MedicalDemo1 that extracts the surface using vtkFlyingEdges3D or vtkMarchingCubes to extract an interpolated isosurface.

Info

The example uses src/Testing/Data/FullHead.mhd which references src/Testing/Data/FullHead.raw.gz.

Other languages

See (PythonicAPI)

Question

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

Code¶

VoxelsOnBoundary.cxx

#include <vtkActor.h>
#include <vtkCamera.h>
#include <vtkDataSetMapper.h>
#include <vtkImageDataToPointSet.h>
#include <vtkImageShrink3D.h>
#include <vtkMetaImageReader.h>
#include <vtkMultiBlockDataSet.h>
#include <vtkMultiThreshold.h>
#include <vtkNamedColors.h>
#include <vtkNew.h>
#include <vtkPlane.h>
#include <vtkProperty.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkRenderer.h>
#include <vtkUnstructuredGrid.h>

#include <iostream>
#include <string>
#include <vector>

int main(int argc, char* argv[])
{
  if (argc < 2)
  {
    cout << "Usage: " << argv[0] << " file.mhd e.g. FullHead.mhd" << endl;
    return EXIT_FAILURE;
  }

  vtkNew<vtkMetaImageReader> reader;
  reader->SetFileName(argv[1]);
  reader->Update();

  vtkNew<vtkImageShrink3D> shrink;
  shrink->SetShrinkFactors(4, 4, 4);
  shrink->SetInputConnection(reader->GetOutputPort());
  shrink->Update();

  // Create a PointSet from the ImageData
  vtkNew<vtkImageDataToPointSet> imageDataToPointSet;
  imageDataToPointSet->SetInputConnection(reader->GetOutputPort());
  imageDataToPointSet->SetInputConnection(shrink->GetOutputPort());
  imageDataToPointSet->Update();

  // Extract voxels on the border between the inside and outside.
  vtkNew<vtkMultiThreshold> threshold;
  // Inside points have one or more points above the isosurface
  int insideId = threshold->AddIntervalSet(
      501, 20000, vtkMultiThreshold::CLOSED, vtkMultiThreshold::CLOSED,
      vtkDataObject::FIELD_ASSOCIATION_POINTS, "ImageScalars", 0, 0);
  // Border points have points that straddle the boundary
  int borderId = threshold->AddIntervalSet(
      499.9999, 501.0000, vtkMultiThreshold::OPEN, vtkMultiThreshold::OPEN,
      vtkDataObject::FIELD_ASSOCIATION_POINTS, "ImageScalars", 0, 0);

  threshold->SetInputData(imageDataToPointSet->GetOutput());

  // Select the intervals to be output
  threshold->OutputSet(insideId);
  threshold->OutputSet(borderId);
  threshold->Update();

  // Visualize
  vtkNew<vtkNamedColors> colors;
  // vtkColor3d outsideColor    = colors->GetColor3d("Crimson");
  vtkColor3d insideColor = colors->GetColor3d("Banana");
  vtkColor3d borderColor = colors->GetColor3d("Mint");
  // vtkColor3d surfaceColor    = colors->GetColor3d("Peacock");
  vtkColor3d backgroundColor = colors->GetColor3d("Silver");

  vtkNew<vtkPlane> plane;
  plane->SetOrigin(dynamic_cast<vtkUnstructuredGrid*>(
                       vtkMultiBlockDataSet::SafeDownCast(
                           threshold->GetOutput()->GetBlock(insideId))
                           ->GetBlock(0))
                       ->GetCenter());
  plane->SetNormal(1.0, -1.0, -1.0);

  // Inside
  vtkNew<vtkDataSetMapper> insideMapper;
  insideMapper->SetInputData(dynamic_cast<vtkUnstructuredGrid*>(
      vtkMultiBlockDataSet::SafeDownCast(
          threshold->GetOutput()->GetBlock(insideId))
          ->GetBlock(0)));
  insideMapper->ScalarVisibilityOff();

  vtkNew<vtkActor> insideActor;
  insideActor->SetMapper(insideMapper);
  insideActor->GetProperty()->SetDiffuseColor(insideColor.GetData());
  insideActor->GetProperty()->SetSpecular(.6);
  insideActor->GetProperty()->SetSpecularPower(30);
  insideActor->GetProperty()->EdgeVisibilityOn();

  // Border
  vtkNew<vtkDataSetMapper> borderMapper;
  borderMapper->SetInputData(dynamic_cast<vtkUnstructuredGrid*>(
      vtkMultiBlockDataSet::SafeDownCast(
          threshold->GetOutput()->GetBlock(borderId))
          ->GetBlock(0)));
  borderMapper->ScalarVisibilityOff();

  vtkNew<vtkActor> borderActor;
  borderActor->SetMapper(borderMapper);
  borderActor->GetProperty()->SetDiffuseColor(borderColor.GetData());
  borderActor->GetProperty()->SetSpecular(.6);
  borderActor->GetProperty()->SetSpecularPower(30);
  borderActor->GetProperty()->EdgeVisibilityOn();

  vtkNew<vtkRenderer> renderer;
  vtkNew<vtkRenderWindow> renderWindow;
  renderWindow->AddRenderer(renderer);
  renderWindow->SetSize(640, 480);

  vtkNew<vtkRenderWindowInteractor> renderWindowInteractor;
  renderWindowInteractor->SetRenderWindow(renderWindow);

  renderer->SetBackground(backgroundColor.GetData());
  renderer->UseHiddenLineRemovalOn();

  renderer->AddActor(insideActor);
  renderer->AddActor(borderActor);

  renderWindow->SetWindowName("CellsOnBoundary");
  renderWindow->Render();

  // Setup a good view
  renderer->GetActiveCamera()->SetViewUp(0, 0, -1);
  renderer->GetActiveCamera()->SetPosition(0, -1, 0);
  renderer->GetActiveCamera()->SetFocalPoint(0, 0, 0);
  renderer->GetActiveCamera()->Azimuth(30.0);
  renderer->GetActiveCamera()->Elevation(30.0);
  renderer->ResetCamera();
  renderer->GetActiveCamera()->Dolly(1.5);

  renderWindowInteractor->Start();

  return EXIT_SUCCESS;
}

CMakeLists.txt¶

cmake_minimum_required(VERSION 3.12 FATAL_ERROR)

project(VoxelsOnBoundary)

find_package(VTK COMPONENTS 
  CommonColor
  CommonCore
  CommonDataModel
  FiltersGeneral
  IOImage
  ImagingCore
  InteractionStyle
  RenderingContextOpenGL2
  RenderingCore
  RenderingFreeType
  RenderingGL2PSOpenGL2
  RenderingOpenGL2
)

if (NOT VTK_FOUND)
  message(FATAL_ERROR "VoxelsOnBoundary: 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(VoxelsOnBoundary MACOSX_BUNDLE VoxelsOnBoundary.cxx )
  target_link_libraries(VoxelsOnBoundary PRIVATE ${VTK_LIBRARIES}
)
# vtk_module_autoinit is needed
vtk_module_autoinit(
  TARGETS VoxelsOnBoundary
  MODULES ${VTK_LIBRARIES}
)

Download and Build VoxelsOnBoundary¶

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

cd VoxelsOnBoundary/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:

./VoxelsOnBoundary

WINDOWS USERS

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