All–
I’ve been having a bit of trouble with the ITKMinimalPathExtraction remote module [1]. In the code below, I create a simple axis-aligned cylinder (foreground 1, background 0), and put my starting and ending points at either end. I tried using the RegularStepGradientDescentOptimizer, but found that the filter’s ability to calculate the path is dependent upon the order in which I provide the start and end points; if I comment out the line that swaps them, the filter fails (zero vertices returned in path). I thought that perhaps this might somehow be related to this issue [2], so I decided to try using the IterateNeighborhoodOptimizer (translating this c++ to python as closely as I could [3]), but I get a segfault on updating the path filter. (I tried enabling faulthandler, but it was unrevealing.) I’m just starting to get familiar with the MinimalPathExtraction module, and was hoping someone with more experience with it might have an idea of how to solve either the first problem or the second.
Best, and thanks,
–Davis
#!/usr/bin/env python
import numpy as np
import itk
import faulthandler
faulthandler.enable()
# Create a simple synthetic image of cylinder oriented along the final numpy axis.
# The inside of the cylinder is 1, the outside is 0.
x, y = np.mgrid[:256, :256]
c = (x - 128) ** 2 + (y - 128) ** 2
c = np.where(c < 1000, 1, 0)
c = np.repeat(np.reshape(c, (256,256,1)), 512, axis = 2)
im = itk.image_from_array(c.astype(np.float32))
IndexType = itk.Index[3]
PointType = itk.Point[itk.D,3]
PathInformationType = itk.SpeedFunctionPathInformation[PointType]
# Defaults suggested in the test suite.
step_length_factor = 1.0
step_length_relax = 0.999
number_of_iterations = 1000
termination_value = 2.0
# Define the start and end points.
s_index = IndexType()
s_index[0] = 0
s_index[1] = 128
s_index[2] = 128
e_index = IndexType()
e_index[0] = 511
e_index[1] = 128
e_index[2] = 128
## If I swap the starting and enxing points, it works as expected.
## If I comment this line, the filter fails (zero vertices are returned).
s_index, e_index = e_index, s_index
# Checking to make sure that my starting and ending points are within the cylinder.
if im.GetPixel(s_index) == 0:
print("WARNING: The starting point is outside the structure of interest.")
if im.GetPixel(e_index) == 0:
print("WARNING: The ending point is outside the structure of interest.")
# Converting from index to physical space.
s_point = im.TransformIndexToPhysicalPoint(s_index)
e_point = im.TransformIndexToPhysicalPoint(e_index)
# Set up and run the filter.
path_information = PathInformationType.New()
path_information.SetStartPoint(s_point)
path_information.SetEndPoint(e_point)
interpolator = itk.LinearInterpolateImageFunction.New(im)
cost_function = itk.SingleImageCostFunction[type(im)].New(interpolator=interpolator)
cost_function.SetInterpolator(interpolator)
spacing = list(im.GetSpacing())
min_spacing = min(spacing)
## Attempting to use the NeighborhoodOptimizer causes a segfault.
#optimizer = itk.IterateNeighborhoodOptimizer.New()
#optimizer.MinimizeOn()
#optimizer.FullyConnectedOn()
#size = optimizer.GetNeighborhoodSize()
#size.SetSize(3)
#for i in range(3):
# size.SetElement(i, im.GetSpacing()[i] * step_length_factor)
#optimizer.SetNeighborhoodSize(size)
# The RegularStepGradientDescentOptimizer works, but only if I swap the start/end points.
optimizer = itk.RegularStepGradientDescentOptimizer.New(
number_of_iterations=number_of_iterations,
maximum_step_length=1.0*step_length_factor*min_spacing,
minimum_step_length=0.5*step_length_factor*min_spacing,
relaxation_factor=step_length_relax)
path_filter = itk.SpeedFunctionToPathFilter.New(
im,
cost_function=cost_function,
optimizer=optimizer,
termination_value=termination_value,
)
path_filter.AddPathInformation(path_information)
path_filter.Update()
path = path_filter.GetOutput(0)
if path.GetVertexList().Size() == 0:
print("WARNING: Zero vertices returned.")
for i in range(path.GetVertexList().Size()):
print(path.GetVertexList().GetElement(i))
[1] GitHub - InsightSoftwareConsortium/ITKMinimalPathExtraction: Vessel and tube centerline extraction
[2] Gradient sometimes zero (which prevents a gradient decent optimizer from decenting) · Issue #61 · InsightSoftwareConsortium/ITKMinimalPathExtraction · GitHub
[3] ITKMinimalPathExtraction/MinimalPathTest.h at master · InsightSoftwareConsortium/ITKMinimalPathExtraction · GitHub