Hello there,
I am trying to rotate a CT-scan and then locate the centroid of a segmented lesion. I have a SITK image named image and a SITK image containing a binary segmentation mask named ann.
The idea is as follows:
- Convert both image and annotation to isotropic spacing
- Calculate the centroid indices of the annotation and convert this centroid to physical coordinates
- Apply a rotation transformation (euler transformation) to the image volume.
- Apply same rotation transformation to the physical centroid
- Calculate the new indices of the centroid, by converting the physical coordinates back to indices
For convenience, assume that step 1 has already been completed. I created the following function for step 2 3 4 and 5.
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
"""
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)
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))
# middle_of_scan = [sz/2 for sz in image.GetSize()]
# compute the resampling grid for the transformed image
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)
rotated_physical_centroid = euler_transform.TransformPoint(physical_centroid)
# 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)
new_centroid = resampled_image.TransformPhysicalPointToIndex(rotated_physical_centroid)
return resampled_image, new_centroid
The centroid indices returned by the function seem to be incorrect. Does anyone have an idea on why?
For convenience, I use the following code that plots the centroid as a red cross on top of the axial slice:
def plot_axial_slice_with_centroid(image, centroid):
# Convert the SimpleITK image to a NumPy array for easier manipulation
image_array = sitk.GetArrayFromImage(image)
# Get the axial slice containing the centroid
z, x, y = centroid
axial_slice = image_array[z, :, :]
# Create a plot
plt.imshow(axial_slice, cmap='gray')
# Add a red cross to indicate the centroid
plt.plot(y, x, 'rx', markersize=10)
# Set appropriate axis labels and title
plt.xlabel('Y')
plt.ylabel('X')
plt.title('Axial Slice with Centroid')
# Display the plot
plt.show()
rotated_image, rotated_centroid = rotation3d(image, ann, 0, 0, 0, -1024)
plot_axial_slice_with_centroid(rotated_image, rotated_centroid)
Any help is appreciated, thanks in advance!