3D血管表面重建

Question

我有一个包含一根血管的 3D 血管 free-hand 超声体积，我正在尝试重建血管的表面。 3D 体积由 2D images/B-scans 的堆栈构成，每个 B-scan 中的血管轮廓已被分割；也就是说，我有一个椭圆代表体积中每个 B-scan 中血管的轮廓。我试图通过遵循 'GenerateModelsFromLabels.cxx' (http://www.vtk.org/Wiki/VTK/Examples/Cxx/Medical/GenerateModelsFromLabels) 的 VTK 示例来重建血管的轮廓。然而，结果并不是像我希望的那样从一帧到另一帧是光滑的表面。它是不连续和不规则的，如果椭圆之间的位移很大，则该表面不会连接体积中两个相邻框架之间的血管轮廓。在我的方法中，我基本上使用了 DiscreteMarchingCubes -> WindowedSincPolyDataFilter -> GeometryFilter。

我尝试了 passband、smoothingIterations 和 featureAngle 参数，并获得了以下最佳结果：

如你所见，它不是一个光滑的连续表面，相邻帧之间有很多未插值的"holes"，但没关系。可以做得更好吗？我还尝试使用 3D Delaunay 三角剖分，但它只给了我凸包，这不是我预期的输出。我想知道是否有更好的方法来重建一个表面，该表面紧密跟随从一个 B-scan 到下一个容器的轮廓？

下面显示了一个最小的工作示例：

    vtkSmartPointer<vtkImageData> vesselVolume =
    vtkSmartPointer<vtkImageData>::New();

int totalImages = 210;

for (int z = 0; z < totalImages; z++)
{
    std::string strFile = "E:/datasets/vasc/rendering/contour/" + std::to_string(z + 1) + ".png";

    cv::Mat im = cv::imread(strFile, CV_LOAD_IMAGE_GRAYSCALE);

    if (z == 0)
    {
        vesselVolume->SetExtent(0, im.cols, 0, im.rows, 0, totalImages - 1);
        vesselVolume->SetSpacing(1, 1, 1);
        vesselVolume->SetOrigin(0, 0, 0);
        vesselVolume->AllocateScalars(VTK_UNSIGNED_CHAR, 0);
    }   

    std::vector<cv::Point2i> locations;   // output, locations of non-zero pixels 
    cv::findNonZero(im, locations);

    for (int nzi = 0; nzi < locations.size(); nzi++)
    {
        unsigned char* pixel = static_cast<unsigned char*>(vesselVolume->GetScalarPointer(locations[nzi].x, locations[nzi].y, z));
        pixel[0] = 255;
    }
}

vtkSmartPointer<vtkDiscreteMarchingCubes> discreteCubes =
    vtkSmartPointer<vtkDiscreteMarchingCubes>::New();

discreteCubes->SetInputData(vesselVolume);
discreteCubes->GenerateValues(1, 255, 255);
discreteCubes->ComputeNormalsOn();

vtkSmartPointer<vtkWindowedSincPolyDataFilter> smoother =
    vtkSmartPointer<vtkWindowedSincPolyDataFilter>::New();

unsigned int smoothingIterations = 10;
double passBand = 2;
double featureAngle = 360.0;

smoother->SetInputConnection(discreteCubes->GetOutputPort());
smoother->SetNumberOfIterations(smoothingIterations);
smoother->BoundarySmoothingOff();
//smoother->FeatureEdgeSmoothingOff();
smoother->FeatureEdgeSmoothingOn();
smoother->SetFeatureAngle(featureAngle);
smoother->SetPassBand(passBand);
smoother->NonManifoldSmoothingOn();
smoother->BoundarySmoothingOn();
smoother->NormalizeCoordinatesOn();
smoother->Update();

vtkSmartPointer<vtkThreshold> selector =
    vtkSmartPointer<vtkThreshold>::New();

selector->SetInputConnection(smoother->GetOutputPort());
selector->SetInputArrayToProcess(0, 0, 0,
    vtkDataObject::FIELD_ASSOCIATION_CELLS,
    vtkDataSetAttributes::SCALARS);

vtkSmartPointer<vtkMaskFields> scalarsOff =
    vtkSmartPointer<vtkMaskFields>::New();

// Strip the scalars from the output
scalarsOff->SetInputConnection(selector->GetOutputPort());
scalarsOff->CopyAttributeOff(vtkMaskFields::POINT_DATA,
    vtkDataSetAttributes::SCALARS);
scalarsOff->CopyAttributeOff(vtkMaskFields::CELL_DATA,
    vtkDataSetAttributes::SCALARS);

vtkSmartPointer<vtkGeometryFilter> geometry =
    vtkSmartPointer<vtkGeometryFilter>::New();

geometry->SetInputConnection(scalarsOff->GetOutputPort());
geometry->Update();

vtkSmartPointer<vtkPolyDataMapper> mapper =
    vtkSmartPointer<vtkPolyDataMapper>::New();  
mapper->SetInputConnection(geometry->GetOutputPort());  
mapper->ScalarVisibilityOff();
mapper->Update();

vtkSmartPointer<vtkRenderWindow> renderWindow =
    vtkSmartPointer<vtkRenderWindow>::New();

vtkSmartPointer<vtkRenderWindowInteractor> renderWindowInteractor =
    vtkSmartPointer<vtkRenderWindowInteractor>::New();
renderWindowInteractor->SetRenderWindow(renderWindow);

    vtkSmartPointer<vtkRenderer> renderer =
    vtkSmartPointer<vtkRenderer>::New();

renderWindow->AddRenderer(renderer);
renderer->SetBackground(.2, .3, .4);

vtkSmartPointer<vtkActor> actor =
    vtkSmartPointer<vtkActor>::New();
actor->SetMapper(mapper);
renderer->AddActor(actor);

renderer->ResetCamera();

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

Answer 1

假设您的问题是切片之间的手抖动，改善结果的一种可能方法是应用切片到切片配准。尝试使用 ImageJ 应该很容易。使用切片之间的转换也可以转换标记的图像。然后运行通过当前管道转换后的标签图像。

3D血管表面重建

3D vessel surface reconstruction

vtk

itk