C++ VTK VMTK 提取血管中心线

VMTK是一个基于VTK和ITK的工具包，主要用于血管的3D重建、几何分析、网格生成、血管分割。可以直接官网下载下来按照它的PypeS规则，结合Python，命令行直接使用；也可以下载源码自己编译，在代码中使用。

算法

VMTK提供了一个准确的血管或管状物体的中线生成算法。这个算法是由 Luca Antiga 在他的博士论文中提出，算法的输入是血管的表面数据和中线的起止点。主要思路是用Delaunay三角剖分算法算出血管Voronoi图，图上的点是血管最大内接球的球心，再由提供的起止点，在这些球心点中根据半径找到最短路径。查找最短路径的算法是Fast Marching算法。算法的最后输出可以得到中线上点的坐标和半径。

参考

• 源码：https://github.com/vmtk/vmtk
• Fast Marching算法：https://www.cnblogs.com/shushen/p/5381753.html
• vmtkcenterlines的命令行用法可以看官网：http://www.vmtk.org/tutorials/Centerlines.html
• 这个博主的介绍：https://blog.csdn.net/a15005784320/article/details/103834840

基本使用方法

// include所有需要的class

void getCenterline(vtkSmartPointer<vtkPolyData> data) {
	
	// 1. 得到血管表面的数据
	vtkSmartPointer<vtkDiscreteMarchingCubes> marchingCubes = vtkSmartPointer<vtkDiscreteMarchingCubes>::New();
	marchingCubes->SetInputData(data);
	marchingCubes->GenerateValues(1, 1, 1); // 如果是多标签可以改成多个值
	marchingCubes->Update();
	
	// 根据需要，做一些网格平滑处理
	vtkSmartPointer<vtkWindowedSincPolyDataFilter> smoother = vtkSmartPointer<vtkWindowedSincPolyDataFilter>::New();
	smoother->SetInputConnection(marchingCubes->GetOutputPort());
	smoother->SetNumberOfIterations(8);
	smoother->BoundarySmoothingOn();
	smoother->Update();
	
	// 根据需要，清理一下数据
	vtkSmartPointer<vtkCleanPolyData> cleaner = vtkSmartPointer<vtkCleanPolyData>::New();
	cleaner->SetInputConnection(smoother->GetOutputPort());
	cleaner->Update();
	
	// 2. 生成血管表面的三角形网格数据
	vtkSmartPointer<vtkTriangleFilter> triangleFilter = vtkSmartPointer<vtkTriangleFilter>::New();
    triangleFilter->SetInputConnection(cleaner->GetOutputPort());
    triangleFilter->PassLinesOff();
    triangleFilter->PassVertsOff();
    triangleFilter->Update();
	
	vtkSmartPointer<vtkPolyData> surface = triangleFilter->GetOutput();
	
	// 根据需要，如果数据量过大，可以缩减一些三角形
	// vtkSmartPointer<vtkDecimatePro> decimate = vtkSmartPointer<vtkDecimatePro>::New();
	// decimate->SetInputConnection(triangleFilter->GetOutputPort()); 
	// decimate->SetTargetReduction(0.2); 
	// decimate->PreserveTopologyOn(); 
	// decimate->Update();
	// surface = decimate->GetOutput();
	
	// 根据需要，如果数据量比较小，可以增加一些三角形
	// vtkSmartPointer<vtkLinearSubdivisionFilter> linearSubdivisionFilter = vtkSmartPointer<vtkLinearSubdivisionFilter>::New();
	// linearSubdivisionFilter->SetInputConnection(triangleFilter->GetOutputPort());
	// linearSubdivisionFilter->SetNumberOfSubdivisions(4);
	// linearSubdivisionFilter->Update();
	// surface = linearSubdivisionFilter->GetOutput();
	
	// 根据需要，填补表面网格上的孔洞
	vtkSmartPointer<vtkvmtkCapPolyData> capper = vtkSmartPointer<vtkvmtkCapPolyData>::New();
    capper->SetInputConnection(surface->GetOutputPort());
    capper->SetDisplacement(0);
    capper->SetInPlaneDisplacement(0);
    capper->Update();
	
	// 3. 创建起止点
	vtkSmartPointer<vtkIdList> sourceSeeds, targetSeeds;
	sourceSeeds->InsertNextId(0);
	targetSeeds->InsertNextId(1);
	
	// 4. 计算中心线
	vtkvmtkPolyDataCenterlines* centerlineFilter = vtkvmtkPolyDataCenterlines::New();
    centerlineFilter->SetInputData(capper->GetOutput());
    centerlineFilter->SetSourceSeedIds(sourceSeeds);
    centerlineFilter->SetTargetSeedIds(targetSeeds);
    centerlineFilter->SetRadiusArrayName("MaximumInscribedSphereRadius");
	centerlineFilter->SetCostFunction("1/R");
    centerlineFilter->SetFlipNormals(0);
    centerlineFilter->SetAppendEndPointsToCenterlines(1);
    centerlineFilter->SetSimplifyVoronoi(0);
    centerlineFilter->SetCenterlineResampling(1);
    centerlineFilter->SetResamplingStepLength(1);
    centerlineFilter->Update();
	
	// 5. 得到中线上的点和每个点对应的半径
	vtkSmartPointer<vtkPolyData> output = centerlineFilter->GetOutput();
	vtkDoubleArray* centerlinesRadiusArray = output->GetPointData()->GetArray("MaximumInscribedSphereRadius");
	for (int i = 0; i < output->GetNumberOfPoints(); i++) {
		double* point = new double[3];
		output->GetPoint(i, point);
		double radius = centerlinesRadiusArray->GetValue(i);
	}
	
}