PolyDataToImageDataVisualizationSTL

VTKExamples/Cxx/PolyData/PolyDataToImageDataVisualizationSTL

Code

PolyDataToImageDataVisualizationSTL.cxx

#include <iostream>
#include <string>
#include <cstdlib>
#include <chrono>

#include <vtkVersion.h>
#include <vtkSmartPointer.h>
#include <vtkCleanPolyData.h>
#include <vtkPolyData.h>
#include <vtkPolyDataMapper.h>
#include <vtkPolyDataToImageStencil.h>
#include <vtkPointData.h>
#include <vtkImageConstantPad.h>
#include <vtkTransformFilter.h>
#include <vtkTransform.h>
#include <vtkRenderer.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkActor.h>
#include <vtkProperty.h>
#include <vtkSTLReader.h>
#include <vtkPoints.h>
#include <vtkGlyph3DMapper.h>
#include <vtkCellArray.h>
#include <vtkTransformPolyDataFilter.h>
#include <vtkPolyDataNormals.h>
#include <vtkTriangleFilter.h>
#include <vtkCubeSource.h>
#include <vtkImageData.h>
#include <vtkImageStencilToImage.h>


/**
 * This program reads a STL mesh and converts it into a volume 
 * representation (vtkImageData) where the foreground voxels are 1
 * and the background voxels are 0. The resulting image is visualized via 
 * vtk.
 **/


// function taken from 3d slicer and modified. original version:
// https://github.com/Slicer/Slicer/blob/master/Libs/vtkSegmentationCore/vtkClosedSurfaceToBinaryLabelmapConversionRule.cxx#L118

bool ConvertPolyDataToBinaryLabelMap(
 vtkSmartPointer<vtkPolyData> closedSurfacePolyData, 
 vtkSmartPointer<vtkImageData> binaryLabelMap)
{
  // Check for consistency
  if (closedSurfacePolyData->GetNumberOfPoints() < 2 || closedSurfacePolyData->GetNumberOfCells() < 2)
  {
    std::cout 
      << "Convert: Cannot create binary labelmap from surface with number of points: "
      << closedSurfacePolyData->GetNumberOfPoints() << " and number of cells: " 
      << closedSurfacePolyData->GetNumberOfCells() << std::endl;

      return false;
  }

  // Compute polydata normals
  vtkNew<vtkPolyDataNormals> normalFilter;
  normalFilter->SetInputData(closedSurfacePolyData);
  normalFilter->ConsistencyOn();
//  normalFilter->AutoOrientNormalsOn();
  normalFilter->SplittingOff();

  // Make sure that we have a clean triangle polydata
  vtkNew<vtkTriangleFilter> triangle;
  triangle->SetInputConnection(normalFilter->GetOutputPort());

  // Convert polydata to stencil
  vtkNew<vtkPolyDataToImageStencil> polyDataToImageStencil;
  polyDataToImageStencil->SetInputConnection(triangle->GetOutputPort());
  polyDataToImageStencil->SetOutputSpacing(binaryLabelMap->GetSpacing());
  polyDataToImageStencil->SetOutputOrigin(binaryLabelMap->GetOrigin());
  polyDataToImageStencil->SetOutputWholeExtent(binaryLabelMap->GetExtent());
  polyDataToImageStencil->Update();

  // Convert stencil to image
  vtkNew<vtkImageStencilToImage> imageStencilToImage;
  imageStencilToImage->SetInputConnection(polyDataToImageStencil->GetOutputPort());
  imageStencilToImage->SetOutsideValue(0);
  imageStencilToImage->SetInsideValue(1);
  imageStencilToImage->SetOutput(binaryLabelMap);
  imageStencilToImage->Update();

  return true;
}

int main(int argc, char* argv[])
{
  // constants
  int imageResolutionX = 200;
  int imageResolutionY = 200;
  int imageResolutionZ = 200;

  static constexpr unsigned char inval = 1;
  static constexpr unsigned char outval = 0;


  // get commandline argument
  if(argc < 2)
  {
    std::string commandName(argv[0]);
    std::cout << "Usage: " << commandName << "<stl-filename>" << std::endl;

    return EXIT_SUCCESS;
  }

  // load STL file
  vtkNew<vtkSTLReader> reader;
  reader->SetFileName(argv[1]);

  vtkNew<vtkCleanPolyData>clean;
  clean->SetInputConnection(reader->GetOutputPort());
  clean->Update();

  vtkNew<vtkPolyData> pd;
  pd->DeepCopy(clean->GetOutput());

  // compute bounds for stl mesh polydata
  double bounds[6];
  pd->GetBounds(bounds);


  // vtkImageData for voxel representation storage  
  vtkNew<vtkImageData> voxelImage;

  // Specify the size of the image data
  voxelImage->SetDimensions( imageResolutionX, 
                             imageResolutionY, 
                             imageResolutionZ);

  double spacing[3];                      // desired volume spacing
  spacing[0] = 1.0 / imageResolutionX;
  spacing[1] = 1.0 / imageResolutionY;
  spacing[2] = 1.0 / imageResolutionZ;

  std::cout << "spacing: "
            << spacing[0] << ", "
            << spacing[1] << ", "
            << spacing[2] << std::endl;
  voxelImage->SetSpacing(spacing);
  voxelImage->SetExtent(0, imageResolutionX, 
                        0, imageResolutionY, 
                        0, imageResolutionZ);

  double origin[3];
  origin[0] = bounds[0] + spacing[0] / 2;
  origin[1] = bounds[2] + spacing[1] / 2;
  origin[2] = bounds[4] + spacing[2] / 2;
  voxelImage->SetOrigin(origin);
  voxelImage->AllocateScalars(VTK_UNSIGNED_CHAR, 1);


  // fill the image with background voxels:
  voxelImage->GetPointData()->GetScalars()->Fill(outval);


  // start measure time
  std::chrono::high_resolution_clock::time_point t1 = std::chrono::high_resolution_clock::now();


  // convert to voxelimage
  ConvertPolyDataToBinaryLabelMap(pd, voxelImage);


  // end measure time
  std::chrono::high_resolution_clock::time_point t2 = std::chrono::high_resolution_clock::now();
  auto duration = std::chrono::duration_cast<std::chrono::microseconds>( t2 - t1 ).count();
  cout << "duration: " << static_cast<double>(duration) / 1000000.0 << std::endl;

  /////////////////////////////////////////////////////////////////
  // Visualization
  /////////////////////////////////////////////////////////////////


  // create point data for visualization via vtkGlyph3DMapper
  // based on the example code from 
  // https://www.vtk.org/Wiki/VTK/Examples/Cxx/Visualization/Glyph3DMapper 
  vtkNew<vtkPoints> points;

  unsigned char* voxelImageData = static_cast<unsigned char*>(voxelImage->GetScalarPointer());

  for (int x = 0; x < imageResolutionX; x++)
  {
    for (int y = 0; y < imageResolutionY; y++)
    {
      for (int z = 0; z < imageResolutionZ; z++)
      {
        unsigned char pixelValue = voxelImageData[x + imageResolutionX * (y + imageResolutionY * z)];

        if(pixelValue != outval)
        {
          points->InsertNextPoint(x, y, z);
        }
      }
    }
  }
  points->Print(std::cout);

  vtkNew<vtkPolyData> polydata;
  polydata->SetPoints(points);

  vtkNew<vtkPolyData> glyph;

  vtkNew<vtkCubeSource> cubeSource;

  vtkNew<vtkGlyph3DMapper> glyph3Dmapper;
  glyph3Dmapper->SetSourceConnection(cubeSource->GetOutputPort());
  glyph3Dmapper->SetInputData(polydata);
  glyph3Dmapper->Update();


  vtkNew<vtkActor> actor;
  actor->SetMapper(glyph3Dmapper);
  actor->GetProperty()->SetColor(0.0, 1.0, 0.0);
  actor->GetProperty()->SetOpacity(0.1); // uncomment for transparency

  vtkNew<vtkRenderer> renderer;
  renderer->UseFXAAOn();

  vtkNew<vtkRenderWindow> renderWindow;
  renderWindow->AddRenderer(renderer);

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

  renderer->AddActor(actor);
  renderer->SetBackground(1.0, 1.0, 1.0);

  renderWindow->Render();
  renderWindowInteractor->Start();

  return EXIT_SUCCESS;
}

CMakeLists.txt

cmake_minimum_required(VERSION 2.8)

PROJECT(PolyDataToImageDataVisualizationSTL)

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

add_executable(PolyDataToImageDataVisualizationSTL MACOSX_BUNDLE PolyDataToImageDataVisualizationSTL.cxx )

target_link_libraries(PolyDataToImageDataVisualizationSTL ${VTK_LIBRARIES})

Download and Build PolyDataToImageDataVisualizationSTL

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

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

./PolyDataToImageDataVisualizationSTL

WINDOWS USERS

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