[Insight-users] Rotation of extracted regions of images

[Charlotte Curtis]

On Thu, Aug 20, 2009 at 12:33 PM, Neal R. Harvey <harve at lanl.gov> wrote:
> Any assistance that anyone can provide with my problem would be much appreciated.
> If you can point me to any documents (beyond the textbook) that might enlighten me, that
> would be great. If you have dealt with this problem yourself and found a solution, if you
> could share your wisdom it would also be very much appreciated.
>
> Cheers
>
> Harve

Hi Harve,

I feel your pain and struggled with a similar problem for a similar
length of time.  In the end, I could never get it to work out as it
should on paper, so I came up with the following (somewhat
computationally expensive) workaround:

1) Create a new blank image volume the size of the largest possible
rotated image
2) Paste the extracted region into the centre
3) Apply the rotation
4) Crop off the blank regions

I only rotate around a single axis (x) so this code probably won't be
identical to what you need, but it works for me so hopefully it'll be
of some help:

itk::PasteImageFilter< ImageType3D, ImageType3D, ImageType3D >::Pointer paster =
	itk::PasteImageFilter< ImageType3D, ImageType3D, ImageType3D >::New();

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

itk::CenteredEuler3DTransform< double >::Pointer transform =
itk::CenteredEuler3DTransform< double >::New();

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

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

// Calculate the centre of the volume
ImageType3D::PointType centre;
centre[0] = origin[0] + spacing[0]*(double)(outputsize[0]-1)/2.0;
centre[1] = origin[1] + spacing[1]*(double)(outputsize[1]-1)/2.0;
centre[2] = origin[2] + spacing[2]*(double)(outputsize[2]-1)/2.0;

// Shift the origin over to get the whole image in
origin[0] -= spacing[0]*static_cast<double>(outputsize[0]-size[0])/2.0;
origin[1] -= spacing[1]*static_cast<double>(outputsize[1]-size[1])/2.0;
origin[2] -= spacing[2]*static_cast<double>(outputsize[2]-size[2])/2.0;

// Create a blank image volume to extend the bounds
ImageType3D::Pointer boundingimg = ImageType3D::New();
boundingimg->SetSpacing(spacing);
boundingimg->SetOrigin(origin);
ImageType3D::IndexType index;
index[0] = 0;
index[1] = 0;
index[2] = 0;
ImageType3D::RegionType region;
region.SetSize(outputsize);
region.SetIndex(index);
boundingimg->SetRegions( region );
boundingimg->Allocate();
boundingimg->FillBuffer(0);
boundingimg->Update();

// Paste the image in the middle of the new blank image
index[0] = (outputsize[0] - size[0])/2;
index[1] = (outputsize[1] - size[1])/2;
index[2] = (outputsize[2] - size[2])/2;
paster->SetDestinationIndex(index);
paster->SetDestinationImage(boundingimg);
paster->SetSourceImage(m_inputimage);
paster->SetSourceRegion(m_inputimage->GetLargestPossibleRegion());
paster->Update();

// Set the values of the transformer
transform->SetCenter(centre);
transformer->SetOutputParametersFromImage(paster->GetOutput());
transform->SetRotation(angle);
transformer->SetTransform(transform);
transformer->SetInput(paster->GetOutput());

// Crop off the blank edges
ImageType3D::RegionType cropregion =
ComputeAutoCropRegion(transformer->GetOutput());
itk::RegionOfInterestImageFilter< ImageType2D, ImageType2D >::Pointer =
	itk::RegionOfInterestImageFilter< ImageType2D, ImageType2D >::New();
cropper->SetInput(transformer->GetOutput());
cropper->SetRegionOfInterest(cropregion);
cropper->update

The "ComputeAutoCropRegion" function just iterates line by line and
determines where the blank image ends and the rotated image begins.
In my case the 3D rotated image is actually projected to a 2D image,
so it's pretty simple, but I imagine you could do the same thing in
3D.

In any case, I hope this helps you a little.  It's not a nice clean
solution but it seems to work...

Charlotte