# Generate Rotational Maximum Intensity Projections

**URL:** https://discourse.itk.org/t/generate-rotational-maximum-intensity-projections/3226
**Category:** Algorithms
**Tags:** python, rotation, simpleitk
**Created:** [June 25, 2020, 8:17am UTC](https://discourse.itk.org/t/generate-rotational-maximum-intensity-projections/3226 "2020-06-25T08:17:16Z")
**Posts on this page:** 1
**Showing post:** 10

<div class="post-metadata">

### Author: ![divya\_chou](https://discourse.itk.org/user_avatar/discourse.itk.org/divya_chou/32/1593_2.png) [@divya\_chou](https://discourse.itk.org/u/divya_chou)
#### Post date: [June 30, 2020, 5:58pm UTC](https://discourse.itk.org/t/generate-rotational-maximum-intensity-projections/3226/10 "2020-06-30T17:58:42Z")

</div>

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](https://simpleitk.org/doxygen/latest/html/classitk_1_1simple_1_1Image.html#a639da24d8df6d41989447d2ec95cb51d) 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:

1. 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)

2. Could you explain what the parameters used while defining the Grid mean? The documentation is not very clear about this.

3. While performing the rotation which interpolator should one use?

4. 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.

---

_[View the full topic](https://discourse.itk.org/t/generate-rotational-maximum-intensity-projections/3226)._