Approximate anatomical axis (curve) of long bone

Development
1

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:

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 mean-centered 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((distance-parcelationLength/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?

2

For long curved structures I would use VMTK for simple geometry based or Voronoi diagram based centerline extraction. The resulting curve can be used to reslice and quantify the cross-sections using Cross-section analysis module.

Maybe you could also discuss applications with @jmhuie (developer of SegmentGeometry extension), who so far worked with straight axes but might be interested in curved axes as well.

3

This seems like a neat feature.

Ultimately, it could be cool to implement the ability to slice through a curve structure, while following the trajectory curvature. Although, I am not sure that this the way to do it. It could be used to measure the length of a curved structure for sure.

4

You can already do this using SlicerVMTK extension. You can even straighten a curved structure (preserving length and cross-section shape) using Sandbox extension’s Curved Planar Reformat module.

5

I was attempting to give this script a try, but I get an error that slicer.vtkFastMarchingGeodesicDistance() does not exist. Might I be missing a needed extension? If so, which one?

6

You just need latest Slicer Preview Release.
It’s the same filter used by DynamicModeler’s SelectByPoints tool

1 Like
7

OK, thanks, I have been using 2021-10-14, so I’ll update again.