Don't lose data with rotation?

Beginner Questions
1

Hello. I’m rotating CTs with a Euler3DTransform. Currently any data that rotates out of the image array is lost. I’ve compensated for this by adding blank columns/rows/slices before rotating, and then taking away blank data after. That works, but slows down the transform considerably. Is there a way to configure the transform so that no data is lost? Thanks!

2

Hello @radiplab,

This is not about changing the transformation, it is about how we think of images in ITK/SimpleITK. They are not arrays of pixels, they are physical objects that occupy a region in physical space, a sampling of physical space on a regular grid. This is a fundamental tenet of ITK/SimpleITK. I highly recommend reading the linked page.

Now that we’ve started you on the path to insight, below is code that does what you requested:

import SimpleITK as sitk
from itertools import product
import numpy as np

image = sitk.ReadImage('training_001_ct.mha')

#Rotate around image center
tx = sitk.Euler3DTransform()
tx.SetCenter(image.TransformContinuousIndexToPhysicalPoint([(sz-1)/2 for sz in image.GetSize()]))
tx.SetRotation(angleX=0.0, angleY=0.0, angleZ=0.785398)

#Compute the bounding box of the rotated volume
start_index = [0,0,0]
end_index = [sz-1 for sz in image.GetSize()]
physical_corners = [image.TransformIndexToPhysicalPoint(corner) for corner in list(product(*zip(start_index, end_index)))]
transformed_corners = [tx.TransformPoint(corner) for corner in physical_corners]
min_bounds = np.min(transformed_corners,0)
max_bounds = np.max(transformed_corners,0)

# We resample onto an axis aligned grid, with origin defined by min_bounds, using
# the original spacing.
physical_size = max_bounds - min_bounds
spacing = image.GetSpacing()
origin = min_bounds
new_size = [int(np.ceil(phsz/spc)) for phsz, spc in zip(physical_size, spacing)] 
air_HU_value = -1000

new_image = sitk.Resample(image1=image,
                          size=new_size,
                          transform=tx,
                          interpolator=sitk.sitkLinear,
                          outputOrigin=origin,
                          outputSpacing=spacing,
                          outputDirection = [1,0,0,0,1,0,0,0,1],
                          defaultPixelValue=air_HU_value) 

sitk.Show(new_image)
4 Likes
3

Great - thanks for your response.