[Insight-users] CenteredRigid2DTransform/ResampleImageFilter inconsistent behavior of transformPoint() and image transform(bug?)

[Matthias Seise]

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,
-------------- next part --------------
A non-text attachment was scrubbed...
Name: out.png
Type: image/png
Size: 737 bytes
Desc: not available
URL: <http://www.itk.org/pipermail/insight-users/attachments/20090831/52b39d55/attachment-0002.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: white.png
Type: image/png
Size: 163 bytes
Desc: not available
URL: <http://www.itk.org/pipermail/insight-users/attachments/20090831/52b39d55/attachment-0003.png>