[Insight-users] CenteredRigid2DTransform/ResampleImageFilter

[Simon Warfield]

You are creating a centered transform which has 5 parameters, and 
initializing it to its default settings.
Then you are setting an angle, which creates a matrix rotating about the 
default center.
Then you set a new center. The documentation says:
> WARNING: When using the Center, we strongly recommend only changing 
> the matrix and translation to define a transform. Changing a 
> transform's center, changes the mapping between spaces - specifically, 
> translation is not changed with respect to that new center, and so the 
> offset is updated to * maintain the consistency with translation. If a 
> center is not used, or is set before the matrix and the offset, then 
> it is safe to change the offset directly. As a rule of thumb, if you 
> wish to set the center explicitly, set before Offset 
> <http://www.itk.org/Doxygen/html/classitk_1_1Offset.html> computations 
> are done

Then you set a new translation.
The transform you are ending up with is not what you are expecting.

The bottom line is that all of the Centered* family of transforms 
operate in this way that frequently surprises new users,
and has a behavior that is hard for optimizers to handle because of the 
interaction between the different parameters.
I recommend not to use any of the Centered* transforms.

Instead, just use the Rigid2DTransform.  Note that the Rigid2DTransform 
has a center that operates as you are expecting.
It is called the 'fixed parameters'.

--
Simon

> Date: Mon, 31 Aug 2009 14:42:13 +0200
> From: Matthias Seise <Matthias.Seise at fit.fraunhofer.de>
> Subject: [Insight-users] CenteredRigid2DTransform/ResampleImageFilter
>         inconsistent behavior of transformPoint() and image transform(bug?)
> To: insight-users at itk.org
> Message-ID: <4A9BC525.5000706 at fit.fraunhofer.de>
> Content-Type: text/plain; charset="iso-8859-15"
>
> Hey all,
>
> I'm struggling with the usage of itkResampleImageFilter together with
> CenteredRigid2DTransform. I just want to transform an image and write
> the full image to a PNG-file.  (I'm using ITK 3.14, "Review ON")
>
> Basic "algorithm" (sample code is attached as well)
> I calculate the bounding box of the transformed image (Basically I take
> the corners of the original image, use
> transform->TransformPoints(corners) to get the new coordinates and take
> min/max coordinates as bounding box)
> Since PNG ignores origin settings I change the translation in the
> transform so that the bounding box has origin (0,0). I set the correct
> size of the output image (=size of bounding box) and the origin (0,0).
>
> The resulting image should show the full rotated image (all corners are
> within the image region which is written to the file) but the file shows
> only a cropped image (wrong translation? see attached image "out.png" )
>
> Here is the code and its output, small sample image are attached for
> testing. Any help is appreciated!
>
> Thanks for the help!
>
> CODE:
>
> #include "itkImage.h"
> #include "itkImageFileReader.h"
> #include "itkImageFileWriter.h"
> #include "itkResampleImageFilter.h"
> #include "itkIdentityTransform.h"
> #include "itkLinearInterpolateImageFunction.h"
> #include "itkRGBPixel.h"
>
> #include "itkCenteredRigid2DTransform.h"
> #include<vector>
>
> const     unsigned int   Dimension = 2;
> typedef   unsigned char  PixelType;
>
> typedef itk::Image< PixelType,  Dimension >   ImageType;
>
>
> typedef itk::ImageFileReader< ImageType >  ReaderType;
> typedef itk::ImageFileWriter< ImageType >  WriterType;
>
> typedef itk::CenteredRigid2DTransform < double > TransformType;
> //  typedef itk::IdentityTransform< double, Dimension >  TransformType;
> typedef itk::ResampleImageFilter<ImageType, ImageType >  FilterType;
> typedef itk::LinearInterpolateImageFunction<ImageType, double >
> InterpolatorType;
>
> int main( int argc, char * argv[] )
> {
>   ReaderType::Pointer reader = ReaderType::New();
>   WriterType::Pointer writer = WriterType::New();
>
>   reader->SetFileName( "white.png" );
>   writer->SetFileName( "out.png" );
>   reader->UpdateOutputInformation();
>
>   TransformType::Pointer transform = TransformType::New();
>   transform->SetAngle(0.3);
>   double centre[2]={25,35};
>   transform->SetCenter(centre);
>   double transl[2]={10,20};
>   transform->SetTranslation(transl);
>
>   ImageType::SpacingType spacing;
>   spacing = reader->GetOutput()->GetSpacing();
>   ImageType::PointType origin;
>   origin = reader->GetOutput()->GetOrigin();
>   ImageType::DirectionType direction;
>   direction=reader->GetOutput()->GetDirection();
>   ImageType::SizeType size;
>   size = reader->GetOutput()->GetLargestPossibleRegion().GetSize();
>
>   std::vector<ImageType::PointType> pIn(4);
>   std::vector<ImageType::PointType> pOut(4);
>   ImageType::PointType orgN;
>
> //calculate corners
>   for (int k=0;k<4;k++){
>     pIn[k].Fill(0);
>   }
>   pIn[1][0]=origin[0];
>   pIn[1][1]=origin[1]+size[1];
>   pIn[2][0]=origin[0]+size[0];
>   pIn[2][1]=origin[1];
>   pIn[3][0]=origin[0]+size[0];
>   pIn[3][1]=origin[1]+size[1];
>
>   std::cout<<"Original origin: "<<origin[0]<<","<<origin[1]<<","<<std::endl;
>   for (int k=0;k<4;k++)
>     std::cout<<"Original corner p"<<k<<":
> "<<pIn[k][0]<<","<<pIn[k][1]<<","<<std::endl;
>
> //Transform all Points
>   orgN=transform->TransformPoint(origin);
>   for (int k=0;k<4;k++)
>     pOut[k]=transform->TransformPoint(pIn[k]);
>
>   TransformType::OutputVectorType trans=transform->GetTranslation();
>
> //calculate bounding box
>   double minX,maxX;
>   for (int z=0;z<2;z++){
>     minX =1000000.0;
>     maxX =-1000000.0;
>     for (int k=0;k<4;k++){
>       if (minX>pOut[k][z])
>         minX=pOut[k][z];
>       if (maxX<pOut[k][z])
>         maxX=pOut[k][z];
>     }
>     size[z] = ((int) ceil(maxX)) - ((int) floor(minX));
>     trans[z]  -= minX;//translate bounding box to origin (0,0)
>   }
>
> // set transform so that bounding box has origin(0,0)
>   transform->SetTranslation(trans);
>   orgN=transform->TransformPoint(origin);
>   for (int k=0;k<4;k++)
>     pOut[k]=transform->TransformPoint(pIn[k]);
>
>   for (int k=0;k<4;k++)
>     std::cout<<"Transformed corner p"<<k<<":
> "<<pOut[k][0]<<","<<pOut[k][1]<<","<<std::endl;
>
>
>   std::cout<<"Size: "<<size[0]<<","<<size[1]<<","<<std::endl;
>   std::cout<<"Origin: "<<origin[0]<<","<<origin[1]<<","<<std::endl;
>
>   FilterType::Pointer filter = FilterType::New();
>   InterpolatorType::Pointer interpolator = InterpolatorType::New();
>
>   filter->SetInterpolator( interpolator );
>   filter->SetDefaultPixelValue( 0  );
>
>   filter->SetTransform( transform );
>
>   filter->SetOutputSpacing( spacing);
>   filter->SetOutputOrigin( origin );
>   filter->SetOutputDirection( direction );
>   filter->SetSize(size);
>
>   filter->SetInput( reader->GetOutput() );
>   writer->SetInput( filter->GetOutput() );
>   writer->Update();
> }
>
> OUTPUT:
> Original origin: 0,0,
> Original corner p0: 0,0,
> Original corner p1: 0,70,
> Original corner p2: 50,0,
> Original corner p3: 50,70,
> Transformed corner p0: 20.6864,0,
> Transformed corner p1: -3.55271e-15,66.8736,
> Transformed corner p2: 68.4532,14.776,
> Transformed corner p3: 47.7668,81.6496,
> Size: 70,82,
> Origin: 0,0,
>   


-- 
Simon