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:
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
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)
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()
)
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!
I tried modifying your long code example to automatically extract centerline curves using python. However, I found that the extracted centerline curves are not labelled in the 3D view, but only in the slice views:
import vtk
import slicer
import ExtractCenterline
end_pts = slicer.util.loadMarkups("endpoints.json")
vessel_seg = slicer.util.loadSegmentation("segmentation.nii.gz")
segmentation = vessel_seg.GetSegmentation()
segmentation.SetConversionParameter("Smoothing factor", "0.1")
vessel_seg.CreateClosedSurfaceRepresentation() # This creates the 3D view
segmentID = segmentation.GetSegmentIdBySegmentName('Segment_1')
segmentVtkPolyData = vtk.vtkPolyData() # create empty vtkPolyData object
vessel_seg.GetClosedSurfaceRepresentation(segmentID, segmentVtkPolyData) # Get vtkPolyData representation of segment surface into segmentVtkPolyData variable
### Extract Centerline ###
extractLogic = ExtractCenterline.ExtractCenterlineLogic()
## Preprocess ##
inputSurfacePolyData = extractLogic.polyDataFromNode(vessel_seg, segmentID)
preprocessEnabled = True # (self._parameterNode.GetParameter("PreprocessInputSurface") == "true")
targetNumberOfPoints = 15000.0 # float(self._parameterNode.GetParameter("TargetNumberOfPoints"))
decimationAggressiveness = 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 = None
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, end_pts, curveSamplingDistance=0.001)
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."
)
Do you have any insights on what might be causing this display error? Moreover, how can I also generate the ‘Centerline Model’ (the model including all the centerline branches and usually displayed in green in the output from the ExtractCenterline extension)?
This section of ExtractCenterline.py shows how the centerline model is generated from the centerlinePolyData:
For the display of the names in the 3D view, I am guessing that something else is obscuring the markups curves (perhaps the voronoi diagram model, or a centerline curve tree model?). It is possible to get the markups curve “Properties label” to show only on the slice views and not in 3D, but it requires setting up multiple display nodes and tuning which views each display node applies to. Since you haven’t done that, the most likely explanation is that something is hiding the markups curves in the 3D view. Try hiding all other elements in the scene except for one of the markups curves (especially all models) and see if you can see the labels then.
Thanks! I tried hidding all the elements, but there is still no label of the branches in the 3D view. However, after running multiple experiments directly from the ExtractCenterline module on the 3D Slicer software, I found that whether a branch gets labelled in the 3D view seems to depend on whether it is visible in the 3D view at the moment when you are running the centerline extraction process (after you click “Apply”). If it was hidden, then there will be no label for it even after extracting the centerline and branches. Similarly, when I ran centerline extraction purely from my python script, technically no branch was visible in the 3D view, so in the end none of them was labelled.
It seems that the display setting is just not updated after extracting the branches.
Hi! I also found out that if I save the centerline curves in a mrml scene file, or just .mrk.json file, when I reload the files, the labels for the curves are not showing. I am struggling to figure out how I can let it show the label in 3D Slicer, and how I can perhaps change the saved .mrk.json file so that when it’s loaded the labels can show. Do you have any idea how to solve these problems?
Hmm, that does seem like odd behavior. If you want to chase this down further, I would try examining
whether the curves end up with multiple display nodes associated with them (myCurveNode.GetNumberOfDisplayNodes()),
whether they are independently assigned to be visible only in certain views (check myCurveNode.GetNthDisplayNode(0).GetViewNodeIDs(), where an empty list means "show in all views), and
what the current visibility setting is for the curve label (myDisplayNode.GetPropertiesLabelVisibility()) . This setting should also be being saved in the .mrk.json file, so you can see if that seems like it is the case, and you can also try turning it off and back on and then save to see if that triggers an update which is then recorded in the saved file.
It seems like all the outputs are as expected, but the centerline curve labels are still only showing in the slice views and not in the 3D view. Do you have any thoughts on this?
I agree, those are the expected outputs. It looks to me like you may have found a bug. If you can generate and share a minimal example which reproduces the problematic behavior that would be very helpful to anyone who is trying to find the source of the problem. As it is, we are now deeper in than I will likely be able to troubleshoot. I would suggest starting a new topic (this one is marked solved so probably has less visibility) on this forum with a clean description and example of how to reproduce the issue, and also file a bug report, though it is not clear to me whether the bug is present in VMTK, Slicer, or possibly vtk. You can link back to this discussion for details on things you’ve checked already.
