[vtkusers] 3D vessel surface reconstruction

[Andras Lasso]

Tejas,

Contours are not aligned to each other well - there is a very visible periodic displacement. You need to do something with this serious problem before trying to reconstruct the surface. If you are reconstructing the volume from a single sweep then the displacements are probably due to cardiac or respiratory motion, so it's not a tracking or calibration issue that you could solve by more careful data acquisition.  I see two main approaches to address this.


A.      Compensate for the displacements. Compute the center-of-gravity of each contour, fit a curve (e.g., polynomial with an order of 2-3) to these center points, then translate each contour so that its center is on the curve. You can then apply VTK's generic-purpose surface reconstruction filters. For a simple tubular shape like a vessel you may simply use vtkDelaunay3D filter to construct a surface and optionally smooth it with vtkWindowedSincPolyDataFilter.

B.      Average displacements. Probably the easiest is to do this in the image domain. Reconstruct a volume from the image slices (you can use the open-source Plus toolkit for this, www.plustoolkit.org<http://www.plustoolkit.org>). Then, segment the vessel on the 3D volume using standard volumetric image segmentation tools (for example, in Segment Editor module in 3D Slicer, www.slicer.org<http://www.slicer.org>).

Andras

From: vtkusers [mailto:vtkusers-bounces at vtk.org] On Behalf Of Paul Korir
Sent: January 23, 2017 4:28
To: vtkusers at vtk.org
Subject: Re: [vtkusers] 3D vessel surface reconstruction


Hi Tejas,

You can check out vtkRuledSurfaceFilter for that. Here are some starters:

http://www.paraview.org/Wiki/VTK/Examples/Python/PolyData/RuledSurfaceFilter

Best,

Paul

On 20/01/2017 21:05, Tejas Sudharshan Mathai wrote:
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






_______________________________________________

Powered by www.kitware.com<http://www.kitware.com>



Visit other Kitware open-source projects at http://www.kitware.com/opensource/opensource.html



Please keep messages on-topic and check the VTK FAQ at: http://www.vtk.org/Wiki/VTK_FAQ



Search the list archives at: http://markmail.org/search/?q=vtkusers



Follow this link to subscribe/unsubscribe:

http://public.kitware.com/mailman/listinfo/vtkusers

--
Paul x4422
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://public.kitware.com/pipermail/vtkusers/attachments/20170123/e49ab0f4/attachment.html>