Hi devs.
I developed a script to parcelate a long bone. This parcells can then be used to calculate their centroids and with that points create the anatomical axis (curve).
Right now the script is missing somo thresholds and consideration of special cases to find the centroids of the parcelated pieces and create the curve from them.
Do you think I should finish it? or there is some limitations I’m not considering that would make this algorithm not work?
It searches the surface’s first principal component vector and uses it as normal of a plane with origin in the centroid of the bone to cut it and use that intersection as input for the geodesicDistanceFilter.
Here are some pictures:

Big parcelations

Small parcelations
Here is the script:
# find anatomical axis curve of bone
import numpy as np
model = getNode('mymodel')
parcelationLength = 9
numberOfSampledPointsOfModel = 2000
if model.GetPolyData().GetNumberOfPoints() > numberOfSampledPointsOfModel:
maskPointsFilter = vtk.vtkMaskPoints()
maskPointsFilter.SetInputData(model.GetPolyData())
#
ratio = round(model.GetPolyData().GetNumberOfPoints()/numberOfSampledPointsOfModel)
#
# This works but the sampling could be biased spatially I think
# If you have vtk9 you should use UNIFORM_SPATIAL_SURFACE option
maskPointsFilter.SetOnRatio(ratio)
maskPointsFilter.RandomModeOn()
maskPointsFilter.Update()
#
polydata = vtk.vtkPolyData()
polydata.ShallowCopy(maskPointsFilter.GetOutput())
#
# Calculate the SVD and mean
from vtk.util.numpy_support import vtk_to_numpy
modelPoints = vtk_to_numpy(polydata.GetPoints().GetData())
else:
from vtk.util.numpy_support import vtk_to_numpy
modelPoints = vtk_to_numpy(model.GetPolyData().GetPoints().GetData())
# Calculate the mean of the points, i.e. the 'center' of the cloud
modelPointsMean = modelPoints.mean(axis=0)
# Do an SVD on the meancentered data.
uu, dd, rightSingularVectors = np.linalg.svd(modelPoints  modelPointsMean)
# Create a frame for the model
modelZ = np.zeros(3)
modelX = rightSingularVectors[0]
# Cut the bone at the half and find the centroid
plane = vtk.vtkPlane()
plane.SetOrigin(modelPointsMean)
plane.SetNormal(modelX)
cutter = vtk.vtkCutter()
cutter.SetInputData(model.GetPolyData())
cutter.SetCutFunction(plane)
cutter.Update()
pointData = cutter.GetOutput().GetPoints().GetData()
from vtk.util.numpy_support import vtk_to_numpy
seedPoints = vtk_to_numpy(pointData)
meshLocator = vtk.vtkPointLocator()
meshLocator.SetDataSet(model.GetMesh())
seeds = vtk.vtkIdList()
for i in range(len(seedPoints)):
pointIDOfClosestPoint = meshLocator.FindClosestPoint(seedPoints[i]);
seeds.InsertNextId(pointIDOfClosestPoint);
geodesicDistance = slicer.vtkFastMarchingGeodesicDistance()
geodesicDistance.SetInputData(model.GetPolyData());
geodesicDistance.SetFieldDataName('GeodesicDistance')
geodesicDistance.SetSeeds(seeds)
geodesicDistance.Update()
parcelatedRegions = vtk.vtkPolyData()
parcelatedRegions.ShallowCopy(geodesicDistance.GetOutput())
pointScalars = parcelatedRegions.GetPointData()
distanceArray = pointScalars.GetArray("GeodesicDistance")
tempParcelationArray = vtk.vtkIntArray()
tempParcelationArray.SetName('TempParcelationArray')
numberOfPoints = parcelatedRegions.GetNumberOfPoints()
tempParcelationArray.SetNumberOfValues(numberOfPoints)
tempParcelationArray.Fill(0)
for i in range(parcelatedRegions.GetNumberOfPoints()):
distance = distanceArray.GetTuple1(i)
if distance < (parcelationLength/2):
parcellID = 0
else:
parcellID = int((distanceparcelationLength/2) // parcelationLength) +1
tempParcelationArray.SetTuple1(i,parcellID)
pointScalars.AddArray(tempParcelationArray)
model.SetAndObservePolyData(parcelatedRegions)
# Set up coloring by selection array
model.GetDisplayNode().SetActiveScalar("TempParcelationArray", vtk.vtkAssignAttribute.POINT_DATA)
model.GetDisplayNode().SetAndObserveColorNodeID("vtkMRMLProceduralColorNodeRandomIntegers")
scalarRange = model.GetDisplayNode().GetScalarRange()
model.GetDisplayNode().SetAutoScalarRange(False)
model.GetDisplayNode().SetScalarRange(1,scalarRange[1])
model.GetDisplayNode().SetScalarVisibility(True)
I wonder how it would work with deformed bones because if it works well it could be used by researchers.
I was thinking that maybe the geodesic distance should be weighted by the inverse of the curvature (there is an option for this in the filter) but I’m not sure it would give me a usable result. I think the propagation speed may vary too much. What do you think?