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

[Karthik Krishnan]

rama wrote:

> 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.

An itk::Transform is not a ProcessObject.  In other words it is not 
intended to process images (which are data objects). It is merely 
intended to represent a transformation from one co-ordinate system to 
another. It is upto a filter (a ProcessObject ) to use the transform to 
process images.

>
> 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.
>
>
> _______________________________________________
> Insight-users mailing list
> Insight-users at itk.org
> http://www.itk.org/mailman/listinfo/insight-users
>