[Insight-users] 3D Affine transformation

[Petter Risholm]

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;
}