I’ve been trying to implement a python method for automatically measuring the aortic diameter from CT segmentations. I was previously trying to do this all in numpy/scipy, but I think I need the efficiency and thought put into the SimpleITK package to improve the output. I’m trying to use SimpleITK’s ResampleImageFilter for the transformation. In my case, the aorta is saved in a binary 3d numpy array, and I have calculated the centerline, it’s tangent vectors, and the vectors orthogonal to the tangents to define the plane. Currently my script looks something like this:
import numpy as np
import SimpleITK as sitk
# Smoothed centerline is the generated centerline from the 3d array, smoothed with Savitzky-Golay filtering using the kimimaro package
centerline_z = smoothed_centerline[:, 0] # Is this the proper way to realign the centerline to itk's orientation?
centerline_y = smoothed_centerline[:, 1]
centerline_x = smoothed_centerline[:, 2]
sitk_centerline = np.column_stack((centerline_x, centerline_y, centerline_z))
sitk_tangents = compute_tangent_vectors(sitk_centerline) # using scipy.interpolate.CubicSpline here
sitk_normals, sitk_binormals = compute_orthogonal_vectors(sitk_tangents)
# tangents, normals, and binormals are all normalized vectors
# Transform the aorta segmentation to sitk image
aorta_itk = sitk.GetImageFromArray(aorta)
# This is my attempt at the function to extract the orthogonal slice using the centerline
def extract_orthogonal_slice_sitk(sitk_image, center, tangent, normal, binormal, slice_resolution):
# Calculate orientation
rotation_matrix = np.vstack([normal, binormal, tangent])
# Create affine transform with rotation and translation
transform = sitk.AffineTransform(3)
transform.SetMatrix(rotation_matrix.ravel())
transform.SetTranslation(np.array(center) - np.array(sitk_image.GetSize()) / 2)
# Define resampler
resampler = sitk.ResampleImageFilter()
resampler.SetReferenceImage(sitk_image)
resampler.SetSize((slice_resolution, slice_resolution, 1))
resampler.SetTransform(transform)
resampler.SetInterpolator(sitk.sitkNearestNeighbor)
# Execute resampling
resampled_image = resampler.Execute(sitk_image)
# Convert back to numpy array
slice_2d = sitk.GetArrayFromImage(resampled_image).squeeze()
return slice_2d
idx = 50
slice_2d = extract_oblique_slice_sitk(aorta_itk, sitk_centerline[idx], sitk_tangents[idx], sitk_normals[idx], sitk_binormals[idx], 100) # this returns an empty array
Any help would be greatly appreciated as I’ve been working on this for some time now without much success. Thanks!
Thank you for your response! I’m not sure if this is exactly what I’m after. I implemented the GetEquivalentEllipsoidDiameter to my aorta_itk image, and I’m confused by the output. Is the first index in the returned tuple the measured diameter?
If so, I know this is not correct as the spacing is in (1mm, 1mm, 1mm) (I resampled as a part of preprocessing), and I had a professional verify that in this image the diameter should be 59mm. What would be best is if I could reconstruct the entire aorta straightened, to verify that my transformations worked correctly. Thanks again for your response!
Perhaps it is easiest to use the danielssondistancemap function to create a distancemap for the inside of the aorta. Then the diameter is twice the value of the distancemap at the centerline points.
import itk
ct_img = itk.imread("ct.nii") # Read in your original CT image so that you know spacing/info
# assuming your aorta array has 1 on the inside and zeros on the outside, we need to flip so
# that it is 0 inside and 1 outside.
aorta = 1 - aorta
# Convert aorta array to itk image and copy info (spacing) from original ct.
# spacing is important so that distances are in physical units, not as voxels.
aorta_img = itk.GetImageFromArray(aorta)
aorta_img.CopyInformation(ct_img)
ddFilter = itk.DanielssonDistanceMapImageFilter(aorta_img)
ddFilter.SetInputIsBinary(True)
ddFilter.SetSquaredDistance(False)
ddFilter.SetUseImageSpacing(True) #important for real-world data
ddFilter.Update()
dist_map_img = ddFilter.GetDistanceMap()
dist_map = itk.GetArrayFromImage(dMap_img)
radius = [ dist_map[centerline_z[i], centerline_y[i], centerline_x[i]] for i in range(len(centerline_x)) ]
diameter = 2 * radius
This will give you the radius of the maximally inscribed sphere at that centerline point. This is NOT the same as fitting an ellipse to the cross section at sharp curves, at sections with rapidly changing radii, at non-circular cross-sections, and at the ends. So, the “correctness” of this approach is application dependent.
Also, you may need to check the ordering of indexing into dMap/aorta. Not certain if it is aorta[z,y,x] or aorta[x,y,z].
And…I haven’t tested this, so there might be bugs, but it should be close.
