Hello @zivy .
I’ve incorporated your inputs and my code looks like this. Right now I’ve hard coded many values as I’m trying to get the code to work for 1 PET Scan. The source of the values are mentioned as a comment.
np_vol = pet_resampled[0]
sitk_vol = sitk.GetImageFromArray(pet_resampled[1])
type(sitk_vol)
#Adding components of spatial objects
sitk_vol.SetDirection([1,0,0,0,1,0,0,0,1]) #Identity Matrix
sitk_vol.SetSpacing([5,5,5]) # Since pixel spacing is 5 x 5 and the slice thickness is 5mm
sitk_vol.SetOrigin([-337.88238150935,-510.94804353575,-768]) # Values from Image Patient Position (IPP) from DICOM
#Function to convert the centre index to origin, from index to millimeter. Created this as was unable to get the TransformContinuousIndexToPhysicalPoint to work
def idx_mm(orig,idx,spacing):
mm_x= idx[0]*spacing[0] + orig[0]
mm_y= idx[1]*spacing[1] + orig[1]
mm_z= idx[2]*spacing[2] + orig[2]
cent = [mm_x,mm_y,mm_z]
return cent
center_of_volume_idx = np.array([np_vol.shape[0]/2,np_vol.shape[1]/2,np_vol.shape[2]/2])
centre_of_vol_mm = idx_mm([-337.88238150935,-510.94804353575,-768],center_of_volume_idx,[5,5,5]) # Output = [47.11761849064999, -170.94804353575, -428.0]
axis = [0,1,0]
degrees = 20
radians = np.pi * degrees / 180
versor_tx = sitk.VersorRigid3DTransform()
versor_tx.SetCenter(centre_of_vol_mm)
versor_tx.SetRotation(axis, radians)
#Defining the grid
grid = sitk.GridImageSource()
grid.SetOutputPixelType(sitk.sitkUInt16)
grid.SetSize([180 X 512 X 512]) #same as np_vol.shape
grid.SetSigma([0.5, 0.5,0.5]) # Used value from code example. What does this mean?
grid.SetGridSpacing([5.0, 5.0,5.0]) #same as Spacing of sitk_vol
grid.SetOrigin([-337.88238150935,-510.94804353575,-768]) #Value of IPP
grid.SetSpacing([5.0, 5.0,5.0]) # Used value from code example. Unclear about difference from GridSpacing
#Performing rotation
rotated_sitk = sitk.Resample(grid, sitk_vol, versor_tx ,sitk.sitkCosineWindowedSinc, 0,sitk.sitkUInt16)
rotated_ndarr= sitk.GetArrayFromImage(rotated_sitk)
plt.imshow(rotated_ndarr[:,:,101]) # Plotting a random slice from rotated object shows Blank
plt.imshow(np_vol[:,:,101]) # Whereas Plotting a random slice from original nd array of the volume showed values of the PET scan
I have the following questions:
-
Are the identified sources of information for the hardcoded values in line with the correct interpretation?
Eg: while defining the sitk volumne as a spatial object,
Origin = IPP
Spacing = Pixel_spacing_X, Pixel_spacing_Y, Slice_Thickness
Direction Cosine matrix = 100(for x), 010 (for y) and 001 (for z) -
Could you explain what the parameters used while defining the Grid mean? The documentation is not very clear about this.
-
While performing the rotation which interpolator should one use?
-
Finally on plotting the rotated object, I get a blank output - This is documented as one of the common errors, and is attributed to setting the wrong grid - But the details of the parameters expected by the Grid and what they mean are fuzzy to me.
Thanks for humoring the plethora of questions. Coming from a Python background, I’m finding it difficult to understand the Doxygen documentation, and am relying on the Jupyter examples which contain very few examples of 3D rotation.