[vtkusers] 3D vessel surface reconstruction

[Tejas Sudharshan Mathai]

Hi,

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:

https://snag.gy/yqK3Fd.jpg
https://snag.gy/txu6Ur.jpg
https://snag.gy/ksJfUa.jpg


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();


Thanks,
Tejas
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://public.kitware.com/pipermail/vtkusers/attachments/20170120/acd0b0bb/attachment.html>