Hello everyone,
I am working with a 3d curved tubular structure. I have extracted the centerline at various points and now I am trying to calculate the diameter of the tube at those points. At the curved regions, I am trying to use a 3D Euler rotation to align the tube with the z-axis so I can compute the diameter.
I am struggling to rotate the tube around the point so that I can extract the cross-section and the transformed point in the rotated space.
My current approach is to calculate the displacement between the two closest points in the centerline. From there, I compute the Euler rotation angles and apply the rotation to the data and a transformation to the center line point.
It appears that my rotation of the data is correct, however I am unable to verify this because the point transformation does not appear correct.
Here is the code I am using (please note the point is in numpy format so I have to change the axis)
def rotation3d(image, theta_x, theta_y, theta_z, point):
"""
This function rotates an image across each of the x, y, z axes by theta_x, theta_y, and theta_z degrees
respectively
:param image: An sitk MRI image
:param theta_x: The amount of degrees the user wants the image rotated around the x axis
:param theta_y: The amount of degrees the user wants the image rotated around the y axis
:param theta_z: The amount of degrees the user wants the image rotated around the z axis
:return: The rotated image
"""
point = (int(point[2]), int(point[1]), int(point[0]))
euler_transform = sitk.Euler3DTransform(image.TransformContinuousIndexToPhysicalPoint([(sz-1)/2.0 for sz in image.GetSize()]), theta_x, theta_y, theta_z, (0, 0, 0))
#euler_transform.SetCenter(point)
euler_transform.SetRotation(theta_x, theta_y, theta_z)
euler_transform.SetComputeZYX(True)
max_indexes = [sz-1 for sz in image.GetSize()]
extreme_indexes = list(itertools.product(*(list(zip([0]*image.GetDimension(),max_indexes)))))
extreme_points_transformed = [euler_transform.TransformPoint(image.TransformContinuousIndexToPhysicalPoint(p)) for p in extreme_indexes]
output_min_coordinates = np.min(extreme_points_transformed, axis=0)
output_max_coordinates = np.max(extreme_points_transformed, axis=0)
output_spacing = image.GetSpacing()
output_origin = output_min_coordinates
output_size = [int(((omx-omn)/ospc)+0.5) for ospc, omn, omx in zip(output_spacing, output_min_coordinates, output_max_coordinates)]
output_direction = [1,0,0,0,1,0,0,0,1]
output_pixeltype = image.GetPixelIDValue()
output = sitk.Resample(image,
output_size,
euler_transform,
sitk.sitkLinear,
output_origin,
output_spacing,
output_direction,
0,
output_pixeltype)
#resampled_image = resample(image, euler_transform)
point = euler_transform.TransformPoint(point)
point = (int(point[2]), int(point[1]), int(point[0]))
return output, point
This code is adapted from SimpleITK Euler3D transform, problem with output size/resampling - #2 by dzenanz
Thank you and any help is appreciated!