[Insight-users] compose an affine matrix and a deformation field

[Jeremy Lecoeur]

Hi,

I'm trying to use the ComposeTransform in order to apply an affine 
matric then a deformation to an image.

Here is my Code:

int main( int argc, char * argv[] )
{
     bool ok = true;
     MatrixType                   matrix3;
     VectorType                   vector3;
     typedef itk::AffineTransform<double,3> Affine3DType;

     const char* input_image = NULL;
     const char* reference_image = NULL;
     const char* output_image = NULL;
     const char* deformation_field = NULL;
     const char* affine_matrix = NULL;


     if (argc != 6) {
         std::cerr << argv[0]
     << "  input_image  reference_image deformation_fiel Affine_matrix 
Output_image";
         return EXIT_FAILURE;
     }

     input_image = argv[1];
     reference_image = argv[2];
     deformation_field = argv[3];
     affine_matrix=argv[4];
     output_image = argv[5];

     if (input_image == NULL) {
         std::cerr << "Input image not specified." << std::endl;
         ok = false;
     }
     if (reference_image == NULL) {
         std::cerr << "Reference image not specified." << std::endl;
         ok = false;
     }
     if (output_image == NULL) {
         std::cerr << "Output image not specified." << std::endl;
         ok = false;
     }
     if (deformation_field == NULL) {
         std::cerr << "Deformation field not specified." << std::endl;
         ok = false;
     }
     if (affine_matrix == NULL) {
         std::cerr << "Deformation field not specified." << std::endl;
         ok = false;
     }
     if (!ok) {
         std::cerr << argv[0]
     << "  input_image  reference_image deformation_field affine_matrix 
output_image";
         return EXIT_FAILURE;
     }


typedef itk::Transform<double, 3, 3> TransformType;
typename TransformType::Pointer transform;

typedef itk::DisplacementFieldTransform<double, 3> 
DisplacementFieldTransformType;


typedef typename DisplacementFieldTransformType::DisplacementFieldType
           DeformationFieldType;
typedef itk::ImageFileReader<DeformationFieldType> 
DeformationFieldReaderType;
typename DeformationFieldReaderType::Pointer field_reader = 
DeformationFieldReaderType::New();
     field_reader->SetFileName(deformation_field);

     try {
         field_reader->Update();
     } catch (itk::ExceptionObject& e) {
         std::cerr << e << std::endl;
         exit(EXIT_FAILURE);
     }

typename DisplacementFieldTransformType::Pointer 
deformationFieldTransform = DisplacementFieldTransformType::New();
deformationFieldTransform->SetDisplacementField( 
field_reader->GetOutput() );
transform = dynamic_cast<TransformType *>( 
deformationFieldTransform.GetPointer() );


loadAffinetransform(affine_matrix,matrix3,vector3);
Affine3DType::Pointer aff3 = Affine3DType::New();
aff3->SetMatrix(matrix3);
aff3->SetOffset(vector3);

typedef itk::CompositeTransform<double, 3> CompositeTransformType;
typename CompositeTransformType::Pointer compotransfo = 
CompositeTransformType::New();
compotransfo->AddTransform(aff3);
compotransfo->AddTransform(transform);

     InputImageReaderType::Pointer input_image_reader = 
InputImageReaderType::New();
     input_image_reader->SetFileName(input_image);

     try {
         input_image_reader->Update();
     } catch (itk::ExceptionObject& e) {
         std::cerr << e << std::endl;
         exit(EXIT_FAILURE);
     }

     InputImageReaderType::Pointer reference_image_reader = 
InputImageReaderType::New();
     reference_image_reader->SetFileName(reference_image);

     try {
         reference_image_reader->Update();
     } catch (itk::ExceptionObject& e) {
         std::cerr << e << std::endl;
         exit(EXIT_FAILURE);
     }

     C2DCastFilterType::Pointer c2d_caster = C2DCastFilterType::New();
     c2d_caster->SetInput(input_image_reader->GetOutput());

     try {
         c2d_caster->Update();
     } catch (itk::ExceptionObject& e) {
         std::cerr << e << std::endl;
         exit(EXIT_FAILURE);
     }

c2d_caster->GetOutput()->SetSpacing(reference_image_reader->GetOutput()->GetSpacing());
c2d_caster->GetOutput()->SetOrigin(reference_image_reader->GetOutput()->GetOrigin());
c2d_caster->GetOutput()->SetDirection(reference_image_reader->GetOutput()->GetDirection());



typedef itk::ResampleImageFilter<InternalImageType, InternalImageType, 
InternalPixelType> ResamplerType;
typename ResamplerType::Pointer resampleFilter = ResamplerType::New();

InterpolatorType::Pointer interpolator = InterpolatorType::New();

resampleFilter->SetInterpolator(interpolator);
resampleFilter->SetInput( c2d_caster->GetOutput() );
resampleFilter->SetOutputParametersFromImage( 
reference_image_reader->GetOutput() );
resampleFilter->SetTransform(compotransfo);

D2CCastFilterType::Pointer d2c_caster = D2CCastFilterType::New();
d2c_caster->SetInput(resampleFilter->GetOutput());

     try {
         d2c_caster->Update();
     } catch (itk::ExceptionObject& e) {
         std::cerr << e << std::endl;
         exit(EXIT_FAILURE);
     }

OutputImageWriterType::Pointer output_image_writer = 
OutputImageWriterType::New();
itk::NiftiImageIO::Pointer imageIO = itk::NiftiImageIO::New();
output_image_writer->SetImageIO(imageIO);
output_image_writer->SetFileName(output_image);
output_image_writer->SetInput(d2c_caster->GetOutput());

     try {
         output_image_writer->Update();
     } catch (itk::ExceptionObject& e) {
         std::cerr << e << std::endl;
         exit(EXIT_FAILURE);
     }

     exit(EXIT_SUCCESS);
}

The affine matrix is well applied (i played with different matrices to 
test that) but the deformation is never applied. The deformation field 
itself is well defined and applied correctly by another executable that 
only applies the deformation field. I am not sure what I am doing wrong.

Jeremy