[Insight-users] Rotating a 3D image without losing information

[Charlotte Curtis]

Hi all,

I'm trying to rotate a 3D image about a single axis (the y axis in
this case).  I'm using a ResampleImageFilter with a
CenteredEuler3DTransform, and I've got a problem where I can't seem to
get the entire image within the output.  Basically, what I'm doing is:
1) Calculating the image centre and setting it using transform->SetCenter()
2) Calculating the largest possible output size (i.e. a rotation of 45
or 135 degrees, resulting in the hypoteneus of the xz plane becoming
the new height/width) and assigning this value to the output size
(resampler->SetSize())
3) Calculating the location of the new image origin based on the new
output size and centre location

All other output parameters are set using SetOutputParametersFromImage.

Everything works okay except for setting the origin.  I've tried so
many different things, including the thing that makes sense as well as
numerous stabs in the dark, and I can't seem to find anything that
works for an arbitrary rotation.  Even for a rotation of 0, I get a
blank image (although messing around with the "origin" code I can get
the image to appear, but I'd prefer it to make sense mathematically).
I'm guessing I just don't understand something about the way itk
images are laid out, but I'm a little stumped at this point.

Here's the relevant section of my code (theta is defined elsewhere, in radians):

const ImageType3D::PointType &origin = m_inputimage->GetOrigin();
const ImageType3D::SpacingType &spacing = m_inputimage->GetSpacing();
const ImageType3D::SizeType size =
m_inputimage->GetLargestPossibleRegion().GetSize();

itk::ResampleImageFilter< ImageType3D, ImageType3D >::Pointer resampler =
	itk::ResampleImageFilter< ImageType3D, ImageType3D >::New();

itk::CenteredEuler3DTransform< double >::Pointer transform =
	itk::CenteredEuler3DTransform< double >::New();
transform->SetIdentity();
resampler->SetOutputParametersFromImage(m_inputimage);

// Find the centre of the image
ImageType3D::PointType centre;
centre[0] = origin[2] + spacing[0]*(double)(size[0]-1)/2.0;
centre[1] = origin[1] + spacing[1]*(double)(size[1]-1)/2.0;
centre[2] = origin[0] + spacing[2]*(double)(size[2]-1)/2.0;
transform->SetCenter(centre);

// Make a bounding box large enough for the largest possible rotated volume
ImageType3D::SizeType outputsize;
outputsize[0] =	ceil(sqrt(pow((double)size[0]*spacing[0], 2.0) +
	pow((double)size[2]*spacing[2], 2.0))/spacing[0]);
outputsize[1] = size[1];
outputsize[2] = outputsize[0]*spacing[0]/spacing[2];
resampler->SetSize(outputsize);

// Need to shift the origin as well
ImageType3D::PointType outputorigin;
outputorigin[0] = centre[0] - spacing[0]*(double)(outputsize[0]-1)/2.0;
outputorigin[1] = centre[1] - spacing[1]*(double)(outputsize[1]-1)/2.0;
outputorigin[2] = centre[2] - spacing[2]*(double)(outputsize[2]-1)/2.0;
resampler->SetOutputOrigin(outputorigin);

transform->SetRotation(0.0, theta, 0.0);
resampler->SetTransform(transform);