VTK曲面重建技术(CPR)

目录

代码

运行：

结果：

参考：

代码

#include <vtkAutoInit.h>
VTK_MODULE_INIT(vtkInteractionStyle)
VTK_MODULE_INIT(vtkRenderingFreeType)
VTK_MODULE_INIT(vtkRenderingOpenGL2)
VTK_MODULE_INIT(vtkRenderingVolumeOpenGL2)#include <vtkActor.h>
#include <vtkCamera.h>
#include <vtkCellArray.h>
#include <vtkDataSetMapper.h>
#include <vtkImageData.h>
#include <vtkImageReader2.h>
#include <vtkImageReader2Factory.h>
#include <vtkNamedColors.h>
#include <vtkNew.h>
#include <vtkPoints.h>
#include <vtkPolyData.h>
#include <vtkPolyDataReader.h>
#include <vtkProbeFilter.h>
#include <vtkProperty.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkRenderer.h>
#include <vtkSmartPointer.h>
#include <vtkSplineFilter.h>
#include <vtkWindowLevelLookupTable.h>#include <sstream>namespace {vtkSmartPointer<vtkPolyData> SweepLine(vtkPolyData* line, double direction[3],double distance, unsigned int cols);
}int main(int argc, char* argv[])
{vtkNew<vtkNamedColors> colors;// Verify argumentsif (argc < 4){std::cout << "Usage: " << argv[0] << " InputVolume PolyDataInput"<< " Resolution" << std::endl;std::cout << "e.g. HeadMRVolume.mhd polyline.vtk 200" << std::endl;return EXIT_FAILURE;}// Parse argumentsstd::string volumeFileName = argv[1];std::string polyDataFileName = argv[2];std::stringstream ssResolution;ssResolution << argv[3];unsigned int resolution;ssResolution >> resolution;// Output argumentsstd::cout << "InputVolume: " << volumeFileName << std::endl<< "PolyDataInput: " << polyDataFileName << std::endl<< "Resolution: " << resolution << std::endl;// Read the volume datavtkNew<vtkImageReader2Factory> imageFactory;vtkSmartPointer<vtkImageReader2> imageReader;imageReader.TakeReference(imageFactory->CreateImageReader2(volumeFileName.c_str()));imageReader->SetFileName(volumeFileName.c_str());imageReader->Update();// Read the PolylinevtkNew<vtkPolyDataReader> polyLineReader;polyLineReader->SetFileName(polyDataFileName.c_str());polyLineReader->Update();vtkNew<vtkSplineFilter> spline;spline->SetInputConnection(polyLineReader->GetOutputPort());spline->SetSubdivideToSpecified();spline->SetNumberOfSubdivisions(resolution);// Sweep the line to form a surfacedouble direction[3];direction[0] = 0.0;direction[1] = 0.0;direction[2] = 1.0;double distance = 164;spline->Update();auto surface =SweepLine(spline->GetOutput(), direction, distance, atoi(argv[3]));// Probe the volume with the extruded surfacevtkNew<vtkProbeFilter> sampleVolume;sampleVolume->SetInputConnection(1, imageReader->GetOutputPort());sampleVolume->SetInputData(0, surface);// Compute a simple window/level based on scalar rangevtkNew<vtkWindowLevelLookupTable> wlLut;double range = imageReader->GetOutput()->GetScalarRange()[1] -imageReader->GetOutput()->GetScalarRange()[0];double level = (imageReader->GetOutput()->GetScalarRange()[1] +imageReader->GetOutput()->GetScalarRange()[0]) /2.0;wlLut->SetWindow(range);wlLut->SetLevel(level);// Create a mapper and actor.vtkNew<vtkDataSetMapper> mapper;mapper->SetInputConnection(sampleVolume->GetOutputPort());mapper->SetLookupTable(wlLut);mapper->SetScalarRange(0, 255);vtkNew<vtkActor> actor;actor->SetMapper(mapper);// Create a renderer, render window, and interactorvtkNew<vtkRenderer> renderer;vtkNew<vtkRenderWindow> renderWindow;renderWindow->AddRenderer(renderer);renderWindow->SetWindowName("CurvedReformation");vtkNew<vtkRenderWindowInteractor> renderWindowInteractor;renderWindowInteractor->SetRenderWindow(renderWindow);// Add the actors to the scenerenderer->AddActor(actor);renderer->SetBackground(colors->GetColor3d("DarkSlateGray").GetData());// Set the camera for viewing medical imagesrenderer->GetActiveCamera()->SetViewUp(0, 0, 1);renderer->GetActiveCamera()->SetPosition(0, 0, 0);renderer->GetActiveCamera()->SetFocalPoint(0, 1, 0);renderer->ResetCamera();// Render and interactrenderWindow->Render();renderWindowInteractor->Start();return EXIT_SUCCESS;
}namespace {vtkSmartPointer<vtkPolyData> SweepLine(vtkPolyData* line, double direction[3],double distance, unsigned int cols){unsigned int rows = line->GetNumberOfPoints();double spacing = distance / cols;vtkNew<vtkPolyData> surface;// Generate the pointscols++;unsigned int numberOfPoints = rows * cols;unsigned int numberOfPolys = (rows - 1) * (cols - 1);vtkNew<vtkPoints> points;points->Allocate(numberOfPoints);vtkNew<vtkCellArray> polys;polys->Allocate(numberOfPolys * 4);double x[3];unsigned int cnt = 0;for (unsigned int row = 0; row < rows; row++){for (unsigned int col = 0; col < cols; col++){double p[3];line->GetPoint(row, p);x[0] = p[0] + direction[0] * col * spacing;x[1] = p[1] + direction[1] * col * spacing;x[2] = p[2] + direction[2] * col * spacing;points->InsertPoint(cnt++, x);}}// Generate the quadsvtkIdType pts[4];for (unsigned int row = 0; row < rows - 1; row++){for (unsigned int col = 0; col < cols - 1; col++){pts[0] = col + row * (cols);pts[1] = pts[0] + 1;pts[2] = pts[0] + cols + 1;pts[3] = pts[0] + cols;polys->InsertNextCell(4, pts);}}surface->SetPoints(points);surface->SetPolys(polys);return surface;}
} // namespace

运行：

demo.exe HeadMRVolume.mhd polyline.vtk 200

The image was generated with this volume data: src/Testing/Data/HeadMRVolume.mhd and src/Testing/Data/HeadMRVolume.raw and this polyLine data: src/Testing/Data/polyline.vtk.

https://github.com/Kitware/vtk-examples/blob/master/src/Testing/Data/HeadMRVolume.mhd

https://github.com/Kitware/vtk-examples/blob/master/src/Testing/Data/polyline.vtk

结果：

参考：

https://kitware.github.io/vtk-examples/site/Cxx/Visualization/CurvedReformation/

 CT 三维重建

http://radiol.dxy.cn/article/510899