Hello!

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!