[Insight-users] 3D Affine transformation

[Luis Ibanez]

Hi Petter,

You seem to be doing a number of manual
computations for making sure that the
rotation is made around the center of the
image.

You should use the CenteredAffineTransform
http://www . itk . org/Insight/Doxygen/html/classitk_1_1CenteredAffineTransform . html

instead of the AffineTransform.
http://www . itk . org/Insight/Doxygen/html/classitk_1_1AffineTransform . html

The CenteredAffine transform will do the
rotation center corrections for you.

You will find a description of this class
in the software guide

   http://www . itk . org/ItkSoftwareGuide . pdf

Section 8.5.3, pdf-page 278.

where it is used in the context of image
registration.




   Regards,


     Luis


----------------------
Petter Risholm wrote:
> hi,
> 
> I'm trying to do a 3D affine transformation on a 3D volume but am having
> huge difficulties getting a resonable result.
> 
> When I run my program without doing any transformations (no rotation and
> no translation), the first slice of the output volume looks exactly like
> the input volume. However, the rest of the slices have no information
> except for the DefaultPixelValue.
> 
> Is there something wrong with the resampling I'm doing?
> 
> This program is just an extension to 3D of the ResampleImageFilter4 found
> in the examples.
> 
> Next follows the metaheader for the volume, next is the code for my
> program.
> 
> regards Petter Risholm
> 
>  NDims = 3 DimSize = 128 128 20
> ElementSpacing = 1 1 1
> Position = 0 0 0
> ElementByteOrderMSB = False
> ElementType = MET_USHORT
> HeaderSize = -1
> ElementDataFile = LIST 2
> Image_ref.raw
> 
> int main( int argc, char * argv[] )
> {
>   if( argc < 5 )
>     {
>     std::cerr << "Usage: " << std::endl;
>     std::cerr << argv[0] << "  inputImageFile  outputImageFile  degrees
> x-trans y-trans z-trans" << std::endl;
>     return 1;
>     }
> 
>   const     unsigned int   Dimension = 3;
>   typedef   unsigned short  InputPixelType;
>   typedef   unsigned short  OutputPixelType;
> 
>   typedef itk::Image< InputPixelType,  Dimension >   InputImageType;
>   typedef itk::Image< OutputPixelType, Dimension >   OutputImageType;
> 
>   typedef itk::ImageFileReader< InputImageType  >  ReaderType;
>   typedef itk::ImageFileWriter< OutputImageType >  WriterType;
>  reader->SetFileName( argv[1] );
>   writer->SetFileName( argv[2] );
> 
>   const double angleInDegrees = atof( argv[3] );
>   const int xtrans = atoi( argv[4] );
>   const int ytrans = atoi( argv[5] );
>   const int ztrans = atoi( argv[6] );
> 
>   typedef itk::ResampleImageFilter<
>                   InputImageType, OutputImageType >  FilterType;
> 
>   FilterType::Pointer filter = FilterType::New();
> 
>   typedef itk::AffineTransform< double, Dimension >  TransformType;
>   TransformType::Pointer transform = TransformType::New();
> 
>   typedef itk::BSplineInterpolateImageFunction<
>                        InputImageType, double >  InterpolatorType;
>   InterpolatorType::Pointer interpolator = InterpolatorType::New();
> 
>   filter->SetInterpolator( interpolator );
> 
>   filter->SetDefaultPixelValue( 13 );
> 
>   reader->Update();
>   const double * spacing = reader->GetOutput()->GetSpacing();
>   const double * origin  = reader->GetOutput()->GetOrigin();
>   InputImageType::SizeType size =
>       reader->GetOutput()->GetLargestPossibleRegion().GetSize();
>   std::cout << "size: " <<  size << "\n";
>   filter->SetOutputOrigin( origin );
>   filter->SetOutputSpacing( spacing );
>   filter->SetSize( size );
> 
> 
>   filter->SetInput( reader->GetOutput() );
>   writer->SetInput( filter->GetOutput() );
> 
>   TransformType::OutputVectorType translation1;
> 
>   const double imageCenterX = origin[0] + spacing[0] * size[0] / 2.0;
>   const double imageCenterY = origin[1] + spacing[1] * size[1] / 2.0;
>   const double imageCenterZ = origin[2] + spacing[2] * size[2] / 2.0;
> 
>   translation1[0] =   -imageCenterX;
>   translation1[1] =   -imageCenterY;
>   translation1[2] =   -imageCenterZ;
> 
>    transform->Translate( translation1 );
> 
>   std::cout << "imageCenterX = " << imageCenterX << std::endl;
>   std::cout << "imageCenterY = " << imageCenterY << std::endl;
>   std::cout << "imageCenterZ = " << imageCenterZ << std::endl;
> TransformType::OutputVectorType axis;
>   axis[0] = 0;
>   axis[1] = 0;
>   axis[2] = 1;
> 
>   const double degreesToRadians = atan(1.0) / 45.0;
>   const double angle = angleInDegrees * degreesToRadians;
>   transform->Rotate3D(axis, angle, false );
> 
>   TransformType::OutputVectorType translation2;
>   translation2[0] =   imageCenterX + xtrans;
>   translation2[1] =   imageCenterY + ytrans;
>   translation2[2] =   imageCenterZ + ztrans;
>   transform->Translate( translation2, false );
> 
>   filter->SetTransform( transform );
> 
>   try
>     {
>     writer->Update();
>     }
>   catch( itk::ExceptionObject & excep )
>     {
>     std::cerr << "Exception catched !" << std::endl;
>     std::cerr << excep << std::endl;
>     }
>   // Software Guide : EndCodeSnippet
> return 0;
> }
> 
> 
> _______________________________________________
> Insight-users mailing list
> Insight-users at itk . org
> http://www . itk . org/mailman/listinfo/insight-users
>