Thanks for your reply!
Why not use LabelGeometryImageFilter’s centroid functionality?
LabelGeometryImageFilter is not included when installing itk using pip (AttributeError: module 'itk' has no attribute 'LabelGeometryImageFilter'). I read somewhere that I have to compile ITK myself if I want that functionality, but I rather use the default ITK and SimpleITK packages.
With linear transforms, it is both better and faster to use TransformGeometryImageFilter than resampling.
Perhaps TransformGeometryImageFilter is better and faster, but so far I have not been able to achieve the desired rotation and translation using that class (perhaps you could give some example code that performs the same transform using TransformGeometryImageFilter as the resampling operation that I have provided?)
Regardless, I did find a mistake in my previously posted function, namely that the centroid coordinates were switched around. I have moved the centroid calculation code to the bottom of the function for clarity. The function is now working for all angles of theta_z!
However, the location of the transformed centroid is wrong when theta_x and theta_y are non-zero, and I do not understand why. Any insight would be appreciated!
def rotation3d(image, ann, theta_x, theta_y, theta_z, background_value=0.0):
"""
This function rotates an image across each of the x, y, z axes by theta_x, theta_y, and theta_z degrees
respectively (euler ZXY orientation) and resamples it to be isotropic.
:param image: An sitk 3D 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
:param background_value: The value that should be used to pad the rotated volume
:return: The rotated image
"""
euler_transform = sitk.Euler3DTransform(
image.TransformContinuousIndexToPhysicalPoint([(sz) / 2.0 for sz in image.GetSize()]),
np.deg2rad(theta_x),
np.deg2rad(theta_y),
np.deg2rad(theta_z))
max_indexes = [sz 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)
# isotropic ouput spacing
output_spacing = min(image.GetSpacing())
output_spacing = [output_spacing] * image.GetDimension()
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()
resampled_image = sitk.Resample(image,
output_size,
euler_transform.GetInverse(),
sitk.sitkLinear,
output_origin,
output_spacing,
output_direction,
background_value,
output_pixeltype)
#Calculate new centroid position
k = sitk.GetArrayFromImage(ann)
spot_indices = np.argwhere(k)
centroid = [int(q) for q in np.round(np.mean(spot_indices.astype(float), axis=0))] # Z,X,Y #
centroid = [centroid[1], centroid[2], centroid[0]] # ZXY to XYZ
physical_centroid = image.TransformIndexToPhysicalPoint(centroid)
rotated_physical_centroid = euler_transform.GetInverse().TransformPoint(physical_centroid)
new_centroid = resampled_image.TransformPhysicalPointToIndex(rotated_physical_centroid)
new_centroid = [new_centroid[2], new_centroid[0], new_centroid[1]] # XYZ to ZXY
return resampled_image, new_centroid
Thank you in advance!