# 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:** 12
**Page:** 1
### Author:  [@divya\_chou](https://discourse.itk.org/u/divya_chou)
#### Post date: [June 25, 2020, 8:17am UTC](https://discourse.itk.org/t/generate-rotational-maximum-intensity-projections/3226/1 "2020-06-25T08:17:16Z")
I am very new to using ITK, and am working using SITK on Python.
I have 3D PET scans. From each 3D matrix, I have generated 2 MIPs at 90 degrees - Sagittal and Coronal. This was easily achieved using the sitk.MaximumProjection function.
How can I generate Maximum Intensity Projection (MIP), at 10 degree rotations?
So that every PET scan gives me 36 MIPs.
---
### Author:  [@zivy](https://discourse.itk.org/u/zivy)
#### Post date: [June 25, 2020, 1:02pm UTC](https://discourse.itk.org/t/generate-rotational-maximum-intensity-projections/3226/2 "2020-06-25T13:02:21Z")
Hello @divya_chou,
ITK/SimpleITK do not have this functionality (general volume rendering). Possibly take a look at [VTK](https://vtk.org/) which is also available in Python. [This very old posting](http://vtk.1045678.n5.nabble.com/multiple-Maximum-Intensity-Projection-volume-rendering-td1227873.html) shows how to do it using VTK (C++). The toolkit has changed since then but this should be a reasonable starting point. For VTK questions see the [VTK discourse](https://discourse.vtk.org/).
---
### Author:  [@divya\_chou](https://discourse.itk.org/u/divya_chou)
#### Post date: [June 25, 2020, 3:16pm UTC](https://discourse.itk.org/t/generate-rotational-maximum-intensity-projections/3226/3 "2020-06-25T15:16:59Z")
Thanks for pointing in the direction of VTK.
However I’m wondering that since ITK offers functionality to generate Coronal and Sagittal MIP (using sitk.MaximumProjection function), and it offers 3D transform using Euler3DTransform, is it not feasible to generate rotational MIPs using ITK?
---
### Author:  [@dchen](https://discourse.itk.org/u/dchen)
#### Post date: [June 25, 2020, 3:26pm UTC](https://discourse.itk.org/t/generate-rotational-maximum-intensity-projections/3226/4 "2020-06-25T15:26:05Z")
You ought to be able to rotate the image, resample it to a new orientation and then generate MIP images of the rotated image to get what you want.
---
### Author:  [@zivy](https://discourse.itk.org/u/zivy)
#### Post date: [June 25, 2020, 3:45pm UTC](https://discourse.itk.org/t/generate-rotational-maximum-intensity-projections/3226/5 "2020-06-25T15:45:42Z")
@dchen’s solves the problem by moving the image, my original thinking involved moving the camera. If you only want to use SimpleITK, I’d go with @dchen’s solution.
---
### Author:  [@divya\_chou](https://discourse.itk.org/u/divya_chou)
#### Post date: [June 26, 2020, 4:49am UTC](https://discourse.itk.org/t/generate-rotational-maximum-intensity-projections/3226/6 "2020-06-26T04:49:18Z")
Could you point me to the exact functions to be used in this case?
---
### Author:  [@zivy](https://discourse.itk.org/u/zivy)
#### Post date: [June 26, 2020, 12:43pm UTC](https://discourse.itk.org/t/generate-rotational-maximum-intensity-projections/3226/7 "2020-06-26T12:43:52Z")
Hello @divya_chou,
So the two relevant elements are the rigid transformation class, and the resample filter.
In this case the most appropriate transformation is [VersorRigid3DTransform](https://simpleitk.org/doxygen/latest/html/classitk_1_1simple_1_1VersorRigid3DTransform.html). You should set the rotation center to the center of your volume, `versor_tx.SetCenter(center_of_volume)`, and set the rotation using the axis-angle representation, `versor_tx.SetRotation(axis, angle)`.
You then use the [ResampleImageFilter](https://simpleitk.org/doxygen/latest/html/classitk_1_1simple_1_1ResampleImageFilter.html) to resample the original volume with the transform. You will need to define the resampling grid.
Finally, you apply the MIP to the resampled volume.
---
### Author:  [@divya\_chou](https://discourse.itk.org/u/divya_chou)
#### Post date: [June 27, 2020, 7:20pm UTC](https://discourse.itk.org/t/generate-rotational-maximum-intensity-projections/3226/8 "2020-06-27T19:20:16Z")
Thanks @zivy.
I followed your answer and have come up with the following code so far. And I have a few followup questions:
```
np_vol = pet_volume #shape is 180 X 512 X 512
sitk_vol = sitk.GetImageFromArray(pet_volume)
center_of_volume = [np_vol.shape[0]/2 , np_vol.shape[1]/2 , np_vol.shape[2]/2] # centre =90,256,256
axis = [0,0,1]
degrees = 20
radians = np.pi * degrees / 180
versor_tx = sitk.VersorRigid3DTransform()
versor_tx.SetCenter(center_of_volume)
versor_tx.SetRotation(axis, radians)
```
1. Am I right in calculating the center\_of\_volume as mid point of each dimension’s number of pixels?
2. How can I define the resampling grid?
3. How to use the ResampleImageFilter?
Many thanks for your guidance.
---
### Author:  [@zivy](https://discourse.itk.org/u/zivy)
#### Post date: [June 29, 2020, 12:51pm UTC](https://discourse.itk.org/t/generate-rotational-maximum-intensity-projections/3226/9 "2020-06-29T12:51:03Z")
Hello @divya_chou,
There are several issues with your code, primarily because you are thinking of an image as an array of pixels and not as a spatial object. This is a [fundamental concept](https://simpleitk.readthedocs.io/en/master/fundamentalConcepts.html) in ITK/SimpleITK.
1. The `sitk_vol` is constructed from a numpy array and is missing all of the spatial information (origin, spacing, direction cosine matrix). Please see [this Jupyter notebook](https://github.com/InsightSoftwareConsortium/SimpleITK-Notebooks/blob/master/Python/03_Image_Details.ipynb) section titled “Conversion between numpy and SimpleITK”.
2. The center of the volume is the midpoint index wise, but you need to convert it to physical coordinates using the image’s [TransformContinuousIndexToPhysicalPoint](https://simpleitk.org/doxygen/latest/html/classitk_1_1simple_1_1Image.html#a639da24d8df6d41989447d2ec95cb51d) .
3. Defining the resampling grid and working with transforms is described in [this Jupyter notebook](https://github.com/InsightSoftwareConsortium/SimpleITK-Notebooks/blob/master/Python/21_Transforms_and_Resampling.ipynb).
---
### Author:  [@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")
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.
---
### Author:  [@zivy](https://discourse.itk.org/u/zivy)
#### Post date: [July 2, 2020, 1:11am UTC](https://discourse.itk.org/t/generate-rotational-maximum-intensity-projections/3226/11 "2020-07-02T01:11:08Z")
Hello @divya_chou,
Below is a short script that does what you want. I think this is as far as I go, you should be able to modify it for your needs:
```auto
import SimpleITK as sitk
import numpy as np
file_name = 'training_001_ct.mha'
image = sitk.ReadImage(file_name)
projection = {'sum': sitk.SumProjection,
'mean': sitk.MeanProjection,
'std': sitk.StandardDeviationProjection,
'min': sitk.MinimumProjection,
'max': sitk.MaximumProjection}
ptype = 'mean'
paxis = 0
rotation_axis = [0,0,1]
rotation_angles = np.linspace(0.0, 2*np.pi, int(360.0/20.0))
rotation_center = image.TransformContinuousIndexToPhysicalPoint([(index-1)/2.0 for index in image.GetSize()])
rotation_transform = sitk.VersorRigid3DTransform()
rotation_transform.SetCenter(rotation_center)
#Compute bounding box of rotating volume and the resampling grid structure
image_indexes = list(zip([0,0,0], [sz-1 for sz in image.GetSize()]))
image_bounds = []
for i in image_indexes[0]:
for j in image_indexes[1]:
for k in image_indexes[2]:
image_bounds.append(image.TransformIndexToPhysicalPoint([i,j,k]))
all_points = []
for angle in rotation_angles:
rotation_transform.SetRotation(rotation_axis, angle)
all_points.extend([rotation_transform.TransformPoint(pnt) for pnt in image_bounds])
all_points = np.array(all_points)
min_bounds = all_points.min(0)
max_bounds = all_points.max(0)
#resampling grid will be isotropic so no matter which direction we project to
#the images we save will always be isotropic (required for image formats that
#assume isotropy - jpg,png,tiff...)
new_spc = [np.min(image.GetSpacing())]*3
new_sz = [int(sz/spc + 0.5) for spc,sz in zip(new_spc, max_bounds-min_bounds)]
proj_images = []
for angle in rotation_angles:
rotation_transform.SetRotation(rotation_axis, angle)
resampled_image = sitk.Resample(image1=image,
size=new_sz,
transform=rotation_transform,
interpolator=sitk.sitkLinear,
outputOrigin=min_bounds,
outputSpacing=new_spc,
outputDirection = [1,0,0,0,1,0,0,0,1],
defaultPixelValue = -1000, #HU unit for air in CT, possibly set to 0 in other cases
outputPixelType = image.GetPixelID())
proj_image = projection[ptype](resampled_image, paxis)
extract_size = list(proj_image.GetSize())
extract_size[paxis]=0
proj_images.append(sitk.Extract(proj_image, extract_size))
# Stack all images into fuax-volume for display
sitk.Show(sitk.JoinSeries(proj_images))
```
---
### Author:  [@divya\_chou](https://discourse.itk.org/u/divya_chou)
#### Post date: [July 7, 2020, 6:38am UTC](https://discourse.itk.org/t/generate-rotational-maximum-intensity-projections/3226/12 "2020-07-07T06:38:20Z")
Hello @zivy
I was able to achieve the objective by using a different overload of the Resample function.
> resampled\_image = sitk.Resample(sitk\_vol, output\_size, versor\_tx, sitk.sitkLinear, output\_origin, output\_spacing, output\_direction)
In this case, I did not have to define the grid. Rather, I just defined the extreme points of the volume needed to be transformed. I found the details to perform this under the section " Defining the Resampling Grid" of this [Jupyter notebook](https://github.com/InsightSoftwareConsortium/SimpleITK-Notebooks/blob/master/Python/21_Transforms_and_Resampling.ipynb)
Many thanks for your guidance towards the right direction.