VMTK Extract Centerline with python console

Hello, I am looking to utilize the vmtk module for the ‘Extract Centerline’ and ‘Cross-section Analysis’ features with a purely scripted method, similar to this issue which had a lot of great information.

A brief summary of my workflow:

  • I use other software to create an aorta segmentation, determine ideal endpoints for the centerline, convert those endpoint coordinates to RAS, and save them to a file ‘Endpoints.mrk.json’
  • I am running “Slicer --python-script extractCenter.py”, where the contents of extractCenter.py are:

CT_vol = slicer.util.loadVolume(‘CT.nii’)
end_pts = slicer.util.loadMarkups(‘Endpoints.mrk.json’)
aorta_seg = slicer.util.loadSegmentation(‘aorta.nii’)
aorta_seg_1 = aorta_seg.GetAttribute(“Segment_1”)
centerline = slicer.util.getModuleLogic(‘ExtractCenterline’).extractCenterline(aorta_seg_1, end_pts)

My issue:

  • I believe aorta_seg_1 is not what I think it is, how can I properly access the Segment_1 vtkDataObject, which is a child of the ‘aorta’ node? Currently if I print out type(aorta_seg_1) I get ‘NoneType’
  • ExtractCenterline has a lot of output parameters which I would like to set, such as a network model or centerline curve. How may I access these purely through python?

Thanks in advance for any advice. Below you can find the specific error I’m getting alongside a screenshot of my current scene.

centerline = slicer.util.getModuleLogic(‘ExtractCenterline’).extractCenterline(aorta_seg_1, end_pts)
Traceback (most recent call last):
File “”, line 1, in
File “C:/Users/me/AppData/Local/slicer.org/Slicer 5.4.0/slicer.org/Extensions-31938/SlicerVMTK/lib/Slicer-5.4/qt-scripted-modules/ExtractCenterline.py”, line 724, in extractCenterline
raise ValueError(“Failed to compute centerline (no Voronoi diagram was generated)”)
ValueError: Failed to compute centerline (no Voronoi diagram was generated)
[VTK] Input port 0 of algorithm vtkvmtkCapPolyData (000001ABAF8E1500) has 0 connections but is not optional.
[VTK] No points to subdivide
[VTK] No points to subdivide
[VTK] No points to subdivide
[VTK] Warning: In vtkDelaunay3D.cxx, line 455
[VTK] vtkDelaunay3D (000001ABA8D6B480): Cannot triangulate; no input points
[VTK] Centerline extraction failed: could not compute surface normals

Here is how to get the segment you want:

segmentation = aorta_seg.GetSegmentation()
segmentID = segmentation.GetSegmentIDBySegmentName('Segment_1')
segment = segmentation.GetSegment(segmentID)

The VMTK inputs/outputs will take a bit more research and I just have time for this quick answer at the moment.

Unfortunately the segment is not the right type:

TypeError: SetInputData argument 1: method requires a vtkDataObject, a vtkSegment was provided.

OK, I looked a little deeper into the ExtractCenterline.py code, and I think what you need is the vtkPolyData for the closed surface representation of the segment.

segmentation = aorta_seg.GetSegmentation()
segmentID = segmentation.GetSegmentIDBySegmentName('Segment_1')
aorta_seg.CreateClosedSurfaceRepresentation() # just to ensure this representation exists
segmentVtkPolyData = vtk.vtkPolyData() # create empty vtkPolyData object
# Get vtkPolyData representation of segment surface into segmentVtkPolyData variable
aorta_seg.GetClosedSurfaceRepresentation(segmentID, segmentVtkPolyData)

For the rest of the details of how to use the VMTK centerline extraction from python inside Slicer, the best place to look is at the code from ExtractCenterline.py, available here: https://github.com/vmtk/SlicerExtension-VMTK/blob/3787ea4a300da28ec5f0824f0715f2713b631155/ExtractCenterline/ExtractCenterline.py

Some longer example code which might be helpful

In some old code, I found a method from inside a custom scripted Slicer module I had developed which uses the ExtractCenterline module’s logic.

    def onGenerateCenterlineButtonClick(self):
        """Trying to mimic the implementation in ExtractCenterline.py from VMTK module"""
        qt.QApplication.setOverrideCursor(qt.Qt.WaitCursor)
        try:
            ## Gather inputs ##
            linkedAirwaySegmentationNode = self._parameterNode.GetNodeReference(
                "LinkedAirwaySegmentationNode"
            )
            linkedAirwaySegmentID = (
                self.ui.linkedAirwaySegmentSelectorWidget.currentSegmentID()
            )
            endPointsMarkupsNode = self._parameterNode.GetNodeReference(
                "EndpointsFiducialNode"
            )
            import ExtractCenterline

            extractLogic = ExtractCenterline.ExtractCenterlineLogic()
            ## Preprocess ##
            inputSurfacePolyData = extractLogic.polyDataFromNode(
                linkedAirwaySegmentationNode, linkedAirwaySegmentID
            )
            if (
                not inputSurfacePolyData
                or inputSurfacePolyData.GetNumberOfPoints() == 0
            ):
                raise ValueError("Valid input surface is required")
            # These are the default values in the VMTK module, which I have never had to change (until I ran into an error)
            # OK, now I've had to change decimationAggressiveness down to 3.5 from 4.0.  Preprocessed surface had inverted segments after decimation
            # which lead to failed centerline curve.  It is possible that we could detect failed curve if it has only two control points? (That was the
            # case in the single failure so far). If so, we could re-run with lower decimation aggressiveness.
            preprocessEnabled = True  # (self._parameterNode.GetParameter("PreprocessInputSurface") == "true")
            targetNumberOfPoints = 5000.0  # float(self._parameterNode.GetParameter("TargetNumberOfPoints"))
            decimationAggressiveness = 3.5  # 4.0 #float(self._parameterNode.GetParameter("DecimationAggressiveness"))
            subdivideInputSurface = False  # (self._parameterNode.GetParameter("SubdivideInputSurface") == "true")
            slicer.util.showStatusMessage(
                "Preprocessing surface before centerline extraction..."
            )
            slicer.app.processEvents()
            preprocessedPolyData = extractLogic.preprocess(
                inputSurfacePolyData,
                targetNumberOfPoints,
                decimationAggressiveness,
                subdivideInputSurface,
            )
            ## Extract centerline ##
            centerlineCurveNode = self._parameterNode.GetNodeReference(
                "CenterlineCurveNode"
            )
            if centerlineCurveNode is None:
                centerlineCurveNode = slicer.mrmlScene.AddNewNodeByClass(
                    "vtkMRMLMarkupsCurveNode", "Centerline curve"
                )
            slicer.util.showStatusMessage("Extracting centerline...")
            slicer.app.processEvents()  # force update
            centerlinePolyData, voronoiDiagramPolyData = extractLogic.extractCenterline(
                preprocessedPolyData, endPointsMarkupsNode
            )
            slicer.util.showStatusMessage(
                "Generating curves and quantification results table..."
            )
            slicer.app.processEvents()  # force update
            centerlinePropertiesTableNode = None
            extractLogic.createCurveTreeFromCenterline(
                centerlinePolyData, centerlineCurveNode, centerlinePropertiesTableNode
            )
            if centerlineCurveNode.GetNumberOfControlPoints() == 2:
                # Extraction had an error, likely due to too high decimation aggressiveness
                # Try again
                slicer.util.errorDisplay(
                    "Centerline generation failed, possibly due to decimationAggressiveness being too high."
                )
                # TODO TODO TODO: implement automatically trying again, with message to user about what's going on
        except Exception as e:
            slicer.util.errorDisplay("Failed to compute results: " + str(e))
            import traceback

            traceback.print_exc()
        qt.QApplication.restoreOverrideCursor()
        slicer.util.showStatusMessage("Centerline analysis complete.", 3000)
        self._parameterNode.SetNodeReferenceID(
            "CenterlineCurveNode", centerlineCurveNode.GetID()
        )
2 Likes

Thank you Mike! Your input is very appreciated, with your help I was able to generate a centerline curve from command line. I was hoping that understanding how to use one part of the vtmk module would make it easier for the next part, but I am stuck again for the cross section analysis tool.
While referencing the CrossSectionAnalysis code I found the functions process() and createMarksupCurve(), but neither of these are callable from CrossSectionAnalysisLogic() because of being NoneType. I noticed the vital code for extractCenterlines was within an onApplyButton() function, but the CrossSectionAnalysis class does fancy stuff with widgets that has thrown me off. Could you offer any more guidance around using these modules through code? Thanks again!

my code:

logic = CrossSectionAnalysis.CrossSectionAnalysisLogic()
input_surface_lumen_data = logic.lumenSurfaceNode(aorta_seg_node, aorta_seg_ID)
slicer.app.processEvents()
crosssection_curve = slicer.mrmlScene.AddNewNodeByClass(“vtkMRMLMarkupsCurveNode”, “Cross-section curve”)
crosssection_analysis_table_node = slicer.mrmlScene.AddNewNodeByClass(“vtkMRMLTableNode”, “Centerline cross-sections”)

You are pointing to a former module that no longer exists as such. Please consider the latest single module ‘Cross-section analysis’ in SlicerVMTK.

1 Like

This is the current code as referenced by @chir.set : https://github.com/vmtk/SlicerExtension-VMTK/blob/master/CrossSectionAnalysis/CrossSectionAnalysis.py
Hopefully that will be helpful. I have not ended up using the CrossSectionAnalysis module at all in my work so I don’t have any experience or outside code that uses it. If you run into a specific problem, it’s fine to post back here and we’ll see if we can help.

1 Like