3D vessel surface reconstruction

Question

I have a 3D vascular free-hand ultrasound volume containing one vessel, and I am trying to reconstruct the surface of the vessel. The 3D volume is constructed from a stack of 2D images/B-scans, and the contour of the vessel in each B-scan has been segmented; that is, I have an ellipse representing the contour of the vessel in each B-scan in the volume. I have tried to reconstruct the contour of the vessel by following the VTK example of 'GenerateModelsFromLabels.cxx' (http://www.vtk.org/Wiki/VTK/Examples/Cxx/Medical/GenerateModelsFromLabels). However, the result is not a smooth surface from one frame to another as I would have hoped for it to be. It is discontinuous and irregular, and the surface doesn't connect the vessel contours between two adjacent frames in the volume if the displacement between the ellipses is large. In my approach, I basically used DiscreteMarchingCubes -> WindowedSincPolyDataFilter -> GeometryFilter.

I played around with the passband, smoothingIterations and featureAngle parameters, and I was able to obtain the best following result:

As you can see, it is not a smooth continuous surface with a lot of uninterpolated "holes" between adjacent frames, but it is all right. Can it be made better? I also tried using a 3D Delaunay triangulation, but it only gave me the convex hull, which is not the output I expected. I would like to know if there is a better approach towards reconstructing a surface that closely follows the contour of the vessel from one B-scan to the next in a volume?

A minimal working example is shown below:

    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

Assuming that your problem is hand shaking between slices, one possible way to improve your result is to apply slice to slice registration. It should be easy to try using ImageJ. Use the transforms between slices to also transform your labeled images. Then run your transformed label images through your current pipeline.