VTK(四)---VMTK血管中心线提取

测试模型获取：
lumen.stl
代码：
CMakeLists.txt

cmake_minimum_required(VERSION 3.5)

project(VesselCenterline LANGUAGES CXX)

set(CMAKE_CXX_STANDARD 11)
set(CMAKE_CXX_STANDARD_REQUIRED ON)

# VMTK
find_package( VMTK REQUIRED)
include(${VMTK_USE_FILE})
# ITK
find_package( ITK REQUIRED)
include(${ITK_USE_FILE})
# VTK
find_package( VTK 8.1.2 EXACT REQUIRED)
include(${VTK_USE_FILE})


add_executable(VesselCenterline main.cpp)

target_link_libraries(VesselCenterline ${VTK_LIBRARIES})

main.cpp:

#include <vtkSmartPointer.h>
#include <vtkRenderWindow.h>
#include <vtkRenderer.h>
#include <vtkPolyDataMapper.h>
#include <vtkActor.h>

#include <vtkJPEGReader.h>
#include <vtkCamera.h>
#include <vtkProperty.h>
#include <vtkAxesActor.h>
#include <vtkSTLReader.h>

#include <vtkDiscreteMarchingCubes.h>
#include <vtkWindowedSincPolyDataFilter.h>
#include <vtkCleanPolyData.h>
#include <vtkTriangleFilter.h>
#include <vtkPolyData.h>
#include <vtkPointData.h>
#include <vtkDataArray.h>
#include <vtkIdList.h>
#include <vtkDecimatePro.h>
#include <vtkDoubleArray.h>
#include <vtkParametricSpline.h>
#include <vtkParametricFunctionSource.h>
#include <vtkLODActor.h>
#include <vtkInteractorStyleTrackballCamera.h>
#include <vtkRenderWindowInteractor.h>

#include <vtkvmtkCapPolyData.h>
#include <vtkvmtkPolyDataCenterlines.h>

#include <iostream>

using namespace std;

vtkSmartPointer<vtkPolyData> ExtractCenterline(vtkSmartPointer<vtkPolyData> inputSurface);

int main(int argc,char**argv)
{
    if(argc<2)
    {
        cout<<"please input a model\n";
        return -1;
    }
    string stl = argv[1];

    vtkSmartPointer<vtkSTLReader> reader = vtkSmartPointer<vtkSTLReader>::New();
    reader->SetFileName(stl.c_str());
    reader->Update();
    // STL 文件
    vtkSmartPointer<vtkPolyData> data = vtkSmartPointer<vtkPolyData>::New();
        data->DeepCopy(reader->GetOutput());
    vtkSmartPointer<vtkPolyDataMapper> mapper = vtkSmartPointer<vtkPolyDataMapper>::New();
        mapper->SetInputConnection(reader->GetOutputPort());
        mapper->Update();
    vtkSmartPointer<vtkActor> actor = vtkSmartPointer<vtkActor>::New();
        actor->SetMapper(mapper);
    vtkSmartPointer<vtkProperty> property = vtkSmartPointer<vtkProperty>::New();
        property->SetDiffuseColor(1, 0.49, 0.25);
        property->SetDiffuse(0.7);
        property->SetSpecular(0.3);
        property->SetSpecularPower(20);
        property->SetOpacity(0.3);
    actor->SetProperty(property); //设置透明度

    //获取中心线 B样条曲线拟合
    vtkSmartPointer<vtkPoints> pathpoints = vtkSmartPointer<vtkPoints>::New();//ExtractCenterline(data);
        pathpoints->DeepCopy(ExtractCenterline(data)->GetPoints());

    //-------------------------------------
        cout<<"The number of points:"<<pathpoints->GetNumberOfPoints()<<endl;
    vtkSmartPointer<vtkParametricSpline> spline = vtkSmartPointer<vtkParametricSpline>::New();
        spline->SetPoints(pathpoints);
        spline->ClosedOff();

    vtkSmartPointer<vtkParametricFunctionSource> splineSource = vtkSmartPointer<vtkParametricFunctionSource>::New();
        splineSource->SetParametricFunction(spline);

    vtkSmartPointer<vtkPolyDataMapper> pointMapper = vtkSmartPointer<vtkPolyDataMapper>::New();
        pointMapper->SetInputConnection(splineSource->GetOutputPort());

    vtkSmartPointer<vtkLODActor> actorPoints = vtkSmartPointer<vtkLODActor>::New();
        actorPoints->SetMapper(pointMapper);
        actorPoints->GetProperty()->SetColor(0,1,0);
        actorPoints->GetProperty()->SetLineWidth(2);
    //-------------------------------------

    //---交互----
    vtkSmartPointer<vtkRenderer> renderer = vtkSmartPointer<vtkRenderer>::New();
        renderer->AddActor(actor);
        renderer->AddActor(actorPoints);
        renderer->SetBackground(0.1,0.2,0.4);; // Background color green

    vtkSmartPointer<vtkRenderWindow> renderWindow =vtkSmartPointer<vtkRenderWindow>::New();
        renderWindow->AddRenderer(renderer);
        renderWindow->SetSize(1920,1080);
    vtkSmartPointer<vtkCamera> aCamera = vtkSmartPointer<vtkCamera>::New();
        aCamera->SetPosition ( 0, 0, -300 );
        aCamera->SetFocalPoint( 0, 0, 0);
        aCamera->SetViewUp ( 0, -1, 0 );
    vtkSmartPointer<vtkRenderWindowInteractor> interactor=vtkSmartPointer<vtkRenderWindowInteractor>::New();
        interactor->SetRenderWindow(renderWindow);

    vtkSmartPointer<vtkInteractorStyleTrackballCamera> style=vtkSmartPointer<vtkInteractorStyleTrackballCamera>::New();
        interactor->SetInteractorStyle(style);

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


    return 0;
}

vtkSmartPointer<vtkPolyData> ExtractCenterline(vtkSmartPointer<vtkPolyData> inputSurface) {

    //1.清理 并 三角化
    vtkSmartPointer<vtkPolyData> cappedSurface = vtkSmartPointer<vtkPolyData>::New();

    vtkSmartPointer<vtkCleanPolyData>surfaceCleaner = vtkSmartPointer<vtkCleanPolyData>::New();
        surfaceCleaner->SetInputData(inputSurface);
        surfaceCleaner->Update();

    vtkSmartPointer < vtkTriangleFilter> surfaceTriangulator = vtkSmartPointer<vtkTriangleFilter>::New();
        surfaceTriangulator->SetInputConnection(surfaceCleaner->GetOutputPort());
        surfaceTriangulator->PassLinesOff();
        surfaceTriangulator->PassVertsOff();
        surfaceTriangulator->Update();

    cappedSurface->DeepCopy(surfaceTriangulator->GetOutput());

    //2.填补管腔 cappedsurface
    vtkSmartPointer<vtkvmtkCapPolyData> capper = vtkSmartPointer<vtkvmtkCapPolyData>::New();
        capper->SetInputData(cappedSurface);
        capper->Update();
    vtkSmartPointer<vtkIdList> CapCenterIds = vtkSmartPointer<vtkIdList>::New();

    cappedSurface->DeepCopy(capper->GetOutput());
    CapCenterIds->DeepCopy(capper->GetCapCenterIds());

    cout<<"The number of the lumen is:" <<capper->GetCapCenterIds()->GetNumberOfIds()<<endl;

    //3.设置起始位置
        //起点
    vtkSmartPointer<vtkIdList> sourceIds = vtkSmartPointer<vtkIdList>::New();
        sourceIds->SetNumberOfIds(1);
        sourceIds->SetId(0, CapCenterIds->GetId(0));

    vtkSmartPointer<vtkIdList> targetIds = vtkSmartPointer<vtkIdList>::New();

        //多个目标
//        targetIds->SetNumberOfIds(CapCenterIds->GetNumberOfIds() - 1);
//        for (int i = 1; i < CapCenterIds->GetNumberOfIds(); i++)
//        {
//            targetIds->SetId(i - 1, CapCenterIds->GetId(i));
//        }
        //单个目标
        targetIds->SetNumberOfIds(1);
        targetIds->SetId(0, CapCenterIds->GetId(2));

    //4.中心线计算 calculate
    vtkSmartPointer<vtkvmtkPolyDataCenterlines> centerlinesFilter = vtkSmartPointer<vtkvmtkPolyDataCenterlines>::New();
        centerlinesFilter->SetInputData(cappedSurface);
        centerlinesFilter->SetSourceSeedIds(sourceIds);
        centerlinesFilter->SetTargetSeedIds(targetIds);
        centerlinesFilter->SetAppendEndPointsToCenterlines(true);
        centerlinesFilter->SetCenterlineResampling(true);
        centerlinesFilter->SetResamplingStepLength(1);
        centerlinesFilter->SetRadiusArrayName("Radius");
        centerlinesFilter->SetEdgeArrayName("Edge");
        centerlinesFilter->SetEdgePCoordArrayName("PCoord");
        centerlinesFilter->Update();

    //获取结果
    vtkSmartPointer<vtkPolyData> output = centerlinesFilter->GetOutput();
    vtkDoubleArray* centerlinesRadiusArray = static_cast<vtkDoubleArray*>(output->GetPointData()->GetArray("Radius"));

    //半径
    for (int i = 0; i < output->GetNumberOfPoints(); i++) {
        double radius = centerlinesRadiusArray->GetValue(i);
    }

    return output;
}

运行结果：

最后
起始sourceIds和终点targetIds可以自行设置。

参考代码：
https://github.com/jackyko1991/Vessel-Centerline-Extraction