[Insight-users] A working code for rotating a 3D array with VersorRigid3DTransform

[rama]

Hi All,

After a very long discussion on the ITK forum regarding rotation of a 3D 
cube, I am able to get it work properly. Thanks to Luis, Karthik, 
Marius, Tamburo and others for their suggestions in this case. Actually 
the problem of memory leak that I raised before is my mistake. I think I 
am overshadowed by the ITK stuff and could not point that memory leak 
before hand. But I was able to solve it. Please refer to messages with 
the title "[Insight-users] confusion with specifying origin for 3D 
affine transformations - itkAffineTransform" for the actual discussion.

Here I am giving the working code snippet, with which one can rotate an 
array cube (a 3D array) about a given center point. Someone might find 
it useful in future. Also if anyone can suggest any optimizations to the 
code below, I will be thankful. One thing which is making me skeptical 
from the beginning is, why should I use ResampleImageFilter when I don't 
need it. Other than using a filter, I did not find any way of feeding in 
data directly to the VersorRigid3DTransfom class. so, any suggestions in 
this case will be helpful.

Here is the working code using the class VersorRigid3DTransform, which 
rotates an array about a given center point around X-Axis.

const unsigned int Dimension=3;
typedef unsigned char InputPixelType;
typedef unsigned char OutputPixelType;
typedef itk::Image<InputPixelType, Dimension> InputImage;
typedef itk::Image<OutputPixelType, Dimension> OutputImage;

typedef itk::ResampleImageFilter<InputImage, OutputImage> Resampler;
Resampler::Pointer resampler=Resampler::New();

typedef itk::ImportImageFilter<InputPixelType, Dimension> 
ImportImgFromArray;   //filter to get data from a C++ array
ImportImgFromArray::Pointer importFilter = ImportImgFromArray::New();

ImportImgFromArray::SizeType size;
size[0]=xPixels;
size[1]=yPixels;
size[2]=zPixels;

ImportImgFromArray::IndexType start;
start.Fill(0);

ImportImgFromArray::RegionType region;
region.SetIndex(start);
region.SetSize(size);
importFilter->SetRegion(region);

double origin[Dimension];
origin[0]=0.0;      //values indicating that origin is set to the 
(0,0,0) index in the array
origin[1]=0.0;
origin[2]=0.0;
importFilter->SetOrigin(origin);

double spacing[Dimension];
spacing[0]=PixelXWidth;
spacing[1]=PixelYWidth;
spacing[2]=PixelZWidth;
importFilter->SetSpacing(spacing);

unsigned char *tempBuffer;
tempBuffer = new unsigned char[size[0]*size[1]*size[2]];

//here comes your data. store your data into this buffer or directly 
specify a pointer of your data in the function below

importFilter->SetImportPointer(tempBuffer, size[0]*size[1]*size[2], false);

typedef itk::NearestNeighborInterpolateImageFunction<InputImage, double> 
Interpolator;
Interpolator::Pointer interpolator=Interpolator::New();
resampler->SetInterpolator(interpolator);

resampler->SetDefaultPixelValue(0);
resampler->SetSize(size);
resampler->SetOutputOrigin(origin);
resampler->SetOutputSpacing(spacing);

typedef itk::VersorRigid3DTransform<double> TransformType;
TransformType::Pointer transform=TransformType::New();

TransformType::InputPointType centerPoint;      //you specify a center 
point here.
centerPoint[0]=(xPixels * PixelXWidth)/2.0;      // here the trick is, 
you are actually shifting
centerPoint[1]=(yPixels * PixelYWidth)/2.0;      //the center point from 
the origin to the point
centerPoint[2]=(zPixels * PixelZWidth)/2.0;      //you want, here it is 
the center point of the cube.
transform->SetCenter(centerPoint);                  //Consider that 
origin initially is at the top-back-left of the cube
                                                                         
//However, this will vary from case to case. This just provides a 
helpful hint.


typedef TransformType::VersorType VersorType;
typedef VersorType::VectorType VectorType;

VersorType rotation;
VectorType axis;

axis[0]=1.0;         //rotate around X-Axis
axis[1]=0.0;
axis[2]=0.0;

rotation.Set(axis, angleInDeg * DegToRad);   //angleInDeg is a variable 
which carries angle of rotation in degrees

transform->SetRotation(rotation);

resampler->SetTransform(transform);

resampler->SetInput(importFilter->GetOutput());
resampler->Update();

InputImage::PixelContainer *container;
container=resampler->GetOutput()->GetPixelContainer();
//container->SetContainerManageMemory(false);   //it is disabled here. 
If you don't want your buffer to be deallocated,
                                                                         
      //uncomment this line here.
unsigned char *tempBuffer2=container->GetImportPointer();

// at this point store the resultant data from tempBuffer2 to the place 
you want. This will be destroyed automatically
//when the 'resampler' filter goes out of scope.

delete [] tempBuffer;

I hope this will help someone else in future in cases of urgency.

thanks and regards,
Rama.