tools/clitkAffineTransformGenericFilter.txx

   1 /*=========================================================================
   2   Program:   vv                     http://www.creatis.insa-lyon.fr/rio/vv
   3
   4   Authors belong to:
   5   - University of LYON              http://www.universite-lyon.fr/
   6   - Léon Bérard cancer center       http://www.centreleonberard.fr
   7   - CREATIS CNRS laboratory         http://www.creatis.insa-lyon.fr
   8
   9   This software is distributed WITHOUT ANY WARRANTY; without even
  10   the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  11   PURPOSE.  See the copyright notices for more information.
  12
  13   It is distributed under dual licence
  14
  15   - BSD        See included LICENSE.txt file
  16   - CeCILL-B   http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
  17   ===========================================================================**/
  18 #ifndef clitkAffineTransformGenericFilter_txx
  19 #define clitkAffineTransformGenericFilter_txx
  20
  21 #include <sstream>
  22 #include <istream>
  23 #include <iterator>
  24 #include <itkCenteredEuler3DTransform.h>
  25 #include "clitkElastix.h"
  26
  27 namespace clitk
  28 {
  29
  30   //-----------------------------------------------------------
  31   // Constructor
  32   //-----------------------------------------------------------
  33   template<class args_info_type>
  34   AffineTransformGenericFilter<args_info_type>::AffineTransformGenericFilter()
  35   {
  36     m_Verbose=false;
  37     m_InputFileName="";
  38   }
  39   //-------------------------------------------------------------------
  40
  41
  42   //-----------------------------------------------------------
  43   // Update
  44   //-----------------------------------------------------------
  45   template<class args_info_type>
  46   void AffineTransformGenericFilter<args_info_type>::Update()
  47   {
  48     // Read the Dimension and PixelType
  49     int Dimension, Components;
  50     std::string PixelType;
  51     ReadImageDimensionAndPixelType(m_InputFileName, Dimension, PixelType, Components);
  52
  53     // Call UpdateWithDim
  54     if(Dimension==2) UpdateWithDim<2>(PixelType, Components);
  55     else
  56       if(Dimension==3) UpdateWithDim<3>(PixelType, Components);
  57       else if (Dimension==4)UpdateWithDim<4>(PixelType, Components);
  58       else {
  59         std::cout<<"Error, Only for 2, 3 or 4  Dimensions!!!"<<std::endl ;
  60         return;
  61       }
  62   }
  63   //-------------------------------------------------------------------
  64
  65
  66   //-------------------------------------------------------------------
  67   // Update with the number of dimensions
  68   //-------------------------------------------------------------------
  69   template<class args_info_type>
  70   template<unsigned int Dimension>
  71   void
  72   AffineTransformGenericFilter<args_info_type>::UpdateWithDim(std::string PixelType, int Components)
  73   {
  74     if (m_Verbose) std::cout << "Image was detected to be "<<Dimension<<"D and "<<Components<<" component(s) of "<<  PixelType<<"..."<<std::endl;
  75
  76     if (Components==1) {
  77       if(PixelType == "short") {
  78         if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and signed short..." << std::endl;
  79         UpdateWithDimAndPixelType<Dimension, signed short>();
  80       }
  81       else if(PixelType == "unsigned_short"){
  82         if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and unsigned_short..." << std::endl;
  83         UpdateWithDimAndPixelType<Dimension, unsigned short>();
  84       }
  85
  86       else if (PixelType == "unsigned_char") {
  87         if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and unsigned_char..." << std::endl;
  88         UpdateWithDimAndPixelType<Dimension, unsigned char>();
  89       }
  90
  91       //     else if (PixelType == "char"){
  92       //       if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and signed_char..." << std::endl;
  93       //       UpdateWithDimAndPixelType<Dimension, signed char>();
  94       //     }
  95       else {
  96         if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and float..." << std::endl;
  97         UpdateWithDimAndPixelType<Dimension, float>();
  98       }
  99     }
 100
 101     else if (Components==3) {
 102       if (m_Verbose) std::cout  << "Launching transform in "<< Dimension <<"D and 3D float (DVF)" << std::endl;
 103       UpdateWithDimAndVectorType<Dimension, itk::Vector<float, Dimension> >();
 104     }
 105
 106     else std::cerr<<"Number of components is "<<Components<<", not supported!"<<std::endl;
 107
 108   }
 109   //-------------------------------------------------------------------
 110
 111
 112   //-------------------------------------------------------------------
 113   // Update with the number of dimensions and the pixeltype
 114   //-------------------------------------------------------------------
 115   template<class args_info_type>
 116   template <unsigned int Dimension, class  PixelType>
 117   void
 118   AffineTransformGenericFilter<args_info_type>::UpdateWithDimAndPixelType()
 119   {
 120
 121     // ImageTypes
 122     typedef itk::Image<PixelType, Dimension> InputImageType;
 123     typedef itk::Image<PixelType, Dimension> OutputImageType;
 124
 125     // Read the input
 126     typedef itk::ImageFileReader<InputImageType> InputReaderType;
 127     typename InputReaderType::Pointer reader = InputReaderType::New();
 128     reader->SetFileName( m_InputFileName);
 129     reader->Update();
 130     typename InputImageType::Pointer input= reader->GetOutput();
 131
 132     //Filter
 133     typedef  itk::ResampleImageFilter< InputImageType,OutputImageType >  ResampleFilterType;
 134     typename ResampleFilterType::Pointer resampler = ResampleFilterType::New();
 135
 136     // Matrix
 137     typename itk::Matrix<double, Dimension+1, Dimension+1> matrix;
 138     if (m_ArgsInfo.rotate_given || m_ArgsInfo.translate_given)
 139       {
 140         if (m_ArgsInfo.matrix_given)
 141           {
 142             std::cerr << "You must use either rotate/translate or matrix options" << std::cout;
 143             return;
 144           }
 145         itk::Array<double> transformParameters(2 * Dimension);
 146         transformParameters.Fill(0.0);
 147         if (m_ArgsInfo.rotate_given)
 148           {
 149             if (Dimension == 2)
 150               transformParameters[0] = m_ArgsInfo.rotate_arg[0];
 151             else
 152               for (unsigned int i = 0; i < 3; i++)
 153                 transformParameters[i] = m_ArgsInfo.rotate_arg[i];
 154           }
 155         if (m_ArgsInfo.translate_given)
 156           {
 157             int pos = 3;
 158             if (Dimension == 2)
 159               pos = 1;
 160             for (unsigned int i = 0; i < Dimension && i < 3; i++)
 161               transformParameters[pos++] = m_ArgsInfo.translate_arg[i];
 162           }
 163         if (Dimension == 4)
 164           {
 165             matrix.SetIdentity();
 166             itk::Matrix<double, 4, 4> tmp = GetForwardAffineMatrix3D(transformParameters);
 167             for (unsigned int i = 0; i < 3; ++i)
 168               for (unsigned int j = 0; j < 3; ++j)
 169                 matrix[i][j] = tmp[i][j];
 170             for (unsigned int i = 0; i < 3; ++i)
 171               matrix[i][4] = tmp[i][3];
 172           }
 173         else
 174           matrix = GetForwardAffineMatrix<Dimension>(transformParameters);
 175       }
 176     else
 177       {
 178         if (m_ArgsInfo.matrix_given)
 179           {
 180             matrix= clitk::ReadMatrix<Dimension>(m_ArgsInfo.matrix_arg);
 181             if (m_Verbose) std::cout << "Reading the matrix..." << std::endl;
 182           }
 183         else {
 184           if (m_ArgsInfo.elastix_given) {
 185             std::vector<std::string> s;
 186             for(uint i=0; i<m_ArgsInfo.elastix_given; i++) s.push_back(m_ArgsInfo.elastix_arg[i]);
 187             matrix = createMatrixFromElastixFile<Dimension>(s, m_Verbose);
 188           }
 189           else
 190             matrix.SetIdentity();
 191         }
 192       }
 193     if (m_Verbose)
 194       std::cout << "Using the following matrix:" << std::endl
 195                 << matrix << std::endl;
 196     typename itk::Matrix<double, Dimension, Dimension> rotationMatrix = clitk::GetRotationalPartMatrix(matrix);
 197     typename itk::Vector<double, Dimension> translationPart = clitk::GetTranslationPartMatrix(matrix);
 198
 199     // Transform
 200     typedef itk::AffineTransform<double, Dimension> AffineTransformType;
 201     typename AffineTransformType::Pointer affineTransform=AffineTransformType::New();
 202     affineTransform->SetMatrix(rotationMatrix);
 203     affineTransform->SetTranslation(translationPart);
 204
 205     // Interp
 206     typedef clitk::GenericInterpolator<args_info_type, InputImageType, double> GenericInterpolatorType;
 207     typename GenericInterpolatorType::Pointer genericInterpolator=GenericInterpolatorType::New();
 208     genericInterpolator->SetArgsInfo(m_ArgsInfo);
 209
 210     // Properties
 211     if (m_ArgsInfo.like_given) {
 212       typename InputReaderType::Pointer likeReader=InputReaderType::New();
 213       likeReader->SetFileName(m_ArgsInfo.like_arg);
 214       likeReader->Update();
 215       resampler->SetOutputParametersFromImage(likeReader->GetOutput());
 216     } else if(m_ArgsInfo.transform_grid_flag) {
 217       typename itk::Matrix<double, Dimension+1, Dimension+1> invMatrix( matrix.GetInverse() );
 218       typename itk::Matrix<double, Dimension, Dimension> invRotMatrix( clitk::GetRotationalPartMatrix(invMatrix) );
 219       typename itk::Vector<double,Dimension> invTrans =  clitk::GetTranslationPartMatrix(invMatrix);
 220
 221       // Display warning
 222       if (m_ArgsInfo.spacing_given)
 223         std::cout << "Warning --spacing ignored (because --transform_grid_flag)" << std::endl;
 224       if (m_ArgsInfo.origin_given)
 225         std::cout << "Warning --origin ignored (because --transform_grid_flag)" << std::endl;
 226
 227       // Spacing is influenced by affine transform matrix and input direction
 228       typename InputImageType::SpacingType outputSpacing;
 229       outputSpacing = invRotMatrix *
 230         input->GetDirection() *
 231         input->GetSpacing();
 232
 233       // Origin is influenced by translation but not by input direction
 234       typename InputImageType::PointType outputOrigin;
 235       outputOrigin = invRotMatrix *
 236         input->GetOrigin() +
 237         invTrans;
 238
 239       // Size is influenced by affine transform matrix and input direction
 240       // Size is converted to double, transformed and converted back to size type.
 241       vnl_vector<double> vnlOutputSize(Dimension);
 242       for(unsigned int i=0; i< Dimension; i++) {
 243         vnlOutputSize[i] = input->GetLargestPossibleRegion().GetSize()[i];
 244       }
 245       vnlOutputSize = invRotMatrix *
 246         input->GetDirection().GetVnlMatrix() *
 247         vnlOutputSize;
 248       typename OutputImageType::SizeType outputSize;
 249       for(unsigned int i=0; i< Dimension; i++) {
 250         // If the size is negative, we have a flip and we must modify
 251         // the origin and the spacing accordingly.
 252         if(vnlOutputSize[i]<0.) {
 253           vnlOutputSize[i] *= -1.;
 254           outputOrigin[i] = outputOrigin[i] + outputSpacing[i] * (vnlOutputSize[i]-1);
 255           outputSpacing[i] *= -1.;
 256         }
 257         outputSize[i] = lrint(vnlOutputSize[i]);
 258       }
 259       resampler->SetSize( outputSize );
 260       resampler->SetOutputSpacing( outputSpacing );
 261       resampler->SetOutputOrigin( outputOrigin );
 262     } else {
 263       //Size
 264       typename OutputImageType::SizeType outputSize;
 265       if (m_ArgsInfo.size_given) {
 266         for(unsigned int i=0; i< Dimension; i++)
 267           outputSize[i]=m_ArgsInfo.size_arg[i];
 268       } else outputSize=input->GetLargestPossibleRegion().GetSize();
 269
 270       //Spacing
 271       typename OutputImageType::SpacingType outputSpacing;
 272       if (m_ArgsInfo.spacing_given) {
 273         for(unsigned int i=0; i< Dimension; i++)
 274           outputSpacing[i]=m_ArgsInfo.spacing_arg[i];
 275       } else outputSpacing=input->GetSpacing();
 276
 277       //Origin
 278       typename OutputImageType::PointType outputOrigin;
 279       if (m_ArgsInfo.origin_given) {
 280         for(unsigned int i=0; i< Dimension; i++)
 281           outputOrigin[i]=m_ArgsInfo.origin_arg[i];
 282       } else outputOrigin=input->GetOrigin();
 283
 284       // Set
 285       resampler->SetSize( outputSize );
 286       resampler->SetOutputSpacing( outputSpacing );
 287       resampler->SetOutputOrigin(  outputOrigin );
 288
 289     }
 290
 291     if (m_ArgsInfo.spacinglike_given) {
 292       typename InputReaderType::Pointer likeReader=InputReaderType::New();
 293       likeReader->SetFileName(m_ArgsInfo.spacinglike_arg);
 294       likeReader->Update();
 295
 296       // set the support like the image
 297       if (m_ArgsInfo.like_given) {
 298         typename OutputImageType::SizeType outputSize;
 299         outputSize[0] = ceil(resampler->GetSize()[0]*resampler->GetOutputSpacing()[0]
 300                              /likeReader->GetOutput()->GetSpacing()[0]);
 301         outputSize[1] = ceil(resampler->GetSize()[1]*resampler->GetOutputSpacing()[1]
 302                              /likeReader->GetOutput()->GetSpacing()[1]);
 303         outputSize[2] = ceil(resampler->GetSize()[2]*resampler->GetOutputSpacing()[2]
 304                              /likeReader->GetOutput()->GetSpacing()[2]);
 305         if (m_ArgsInfo.verbose_flag) {
 306           std::cout << "Compute the number of pixels such as the support is like " << m_ArgsInfo.like_arg << std::endl;
 307         }
 308         resampler->SetSize( outputSize );
 309       }
 310
 311       resampler->SetOutputSpacing( likeReader->GetOutput()->GetSpacing() );
 312     }
 313
 314     if (m_ArgsInfo.verbose_flag) {
 315       std::cout << "Setting the output size to " << resampler->GetSize() << "..." << std::endl;
 316       std::cout << "Setting the output spacing to " << resampler->GetOutputSpacing() << "..." << std::endl;
 317       std::cout << "Setting the output origin to " << resampler->GetOutputOrigin() << "..." << std::endl;
 318     }
 319
 320     resampler->SetInput( input );
 321     resampler->SetTransform( affineTransform );
 322     resampler->SetInterpolator( genericInterpolator->GetInterpolatorPointer());
 323     resampler->SetDefaultPixelValue( static_cast<PixelType>(m_ArgsInfo.pad_arg) );
 324
 325     try {
 326       resampler->Update();
 327     } catch(itk::ExceptionObject) {
 328       std::cerr<<"Error resampling the image"<<std::endl;
 329     }
 330
 331     typename OutputImageType::Pointer output = resampler->GetOutput();
 332
 333     // Output
 334     typedef itk::ImageFileWriter<OutputImageType> WriterType;
 335     typename WriterType::Pointer writer = WriterType::New();
 336     writer->SetFileName(m_ArgsInfo.output_arg);
 337     writer->SetInput(output);
 338     writer->Update();
 339
 340   }
 341   //-------------------------------------------------------------------
 342
 343
 344   //-------------------------------------------------------------------
 345   // Update with the number of dimensions and the pixeltype (components)
 346   //-------------------------------------------------------------------
 347   template<class args_info_type>
 348   template<unsigned int Dimension, class PixelType>
 349   void AffineTransformGenericFilter<args_info_type>::UpdateWithDimAndVectorType()
 350   {
 351     // ImageTypes
 352     typedef itk::Image<PixelType, Dimension> InputImageType;
 353     typedef itk::Image<PixelType, Dimension> OutputImageType;
 354
 355     // Read the input
 356     typedef itk::ImageFileReader<InputImageType> InputReaderType;
 357     typename InputReaderType::Pointer reader = InputReaderType::New();
 358     reader->SetFileName( m_InputFileName);
 359     reader->Update();
 360     typename InputImageType::Pointer input= reader->GetOutput();
 361
 362     //Filter
 363     typedef  itk::VectorResampleImageFilter< InputImageType,OutputImageType, double >  ResampleFilterType;
 364     typename ResampleFilterType::Pointer resampler = ResampleFilterType::New();
 365
 366     // Matrix
 367     typename itk::Matrix<double, Dimension+1, Dimension+1> matrix;
 368     if (m_ArgsInfo.rotate_given || m_ArgsInfo.translate_given)
 369       {
 370         if (m_ArgsInfo.matrix_given)
 371           {
 372             std::cerr << "You must use either rotate/translate or matrix options" << std::cout;
 373             return;
 374           }
 375         itk::Array<double> transformParameters(2 * Dimension);
 376         transformParameters.Fill(0.0);
 377         if (m_ArgsInfo.rotate_given)
 378           {
 379             if (Dimension == 2)
 380               transformParameters[0] = m_ArgsInfo.rotate_arg[0];
 381             else
 382               for (unsigned int i = 0; i < 3; i++)
 383                 transformParameters[i] = m_ArgsInfo.rotate_arg[i];
 384           }
 385         if (m_ArgsInfo.translate_given)
 386           {
 387             int pos = 3;
 388             if (Dimension == 2)
 389               pos = 1;
 390             for (unsigned int i = 0; i < Dimension && i < 3; i++)
 391               transformParameters[pos++] = m_ArgsInfo.translate_arg[i];
 392           }
 393         if (Dimension == 4)
 394           {
 395             matrix.SetIdentity();
 396             itk::Matrix<double, 4, 4> tmp = GetForwardAffineMatrix3D(transformParameters);
 397             for (unsigned int i = 0; i < 3; ++i)
 398               for (unsigned int j = 0; j < 3; ++j)
 399                 matrix[i][j] = tmp[i][j];
 400             for (unsigned int i = 0; i < 3; ++i)
 401               matrix[i][4] = tmp[i][3];
 402           }
 403         else
 404           matrix = GetForwardAffineMatrix<Dimension>(transformParameters);
 405       }
 406     else
 407       {
 408         if (m_ArgsInfo.matrix_given)
 409           {
 410             matrix= clitk::ReadMatrix<Dimension>(m_ArgsInfo.matrix_arg);
 411             if (m_Verbose) std::cout << "Reading the matrix..." << std::endl;
 412           }
 413         else
 414           matrix.SetIdentity();
 415       }
 416     if (m_Verbose)
 417       std::cout << "Using the following matrix:" << std::endl
 418                 << matrix << std::endl;
 419     typename itk::Matrix<double, Dimension, Dimension> rotationMatrix = clitk::GetRotationalPartMatrix(matrix);
 420     typename itk::Vector<double, Dimension> translationPart = clitk::GetTranslationPartMatrix(matrix);
 421
 422     // Transform
 423     typedef itk::AffineTransform<double, Dimension> AffineTransformType;
 424     typename AffineTransformType::Pointer affineTransform=AffineTransformType::New();
 425     affineTransform->SetMatrix(rotationMatrix);
 426     affineTransform->SetTranslation(translationPart);
 427
 428     // Interp
 429     typedef clitk::GenericVectorInterpolator<args_info_type, InputImageType, double> GenericInterpolatorType;
 430     typename GenericInterpolatorType::Pointer genericInterpolator=GenericInterpolatorType::New();
 431     genericInterpolator->SetArgsInfo(m_ArgsInfo);
 432
 433     // Properties
 434     if (m_ArgsInfo.like_given) {
 435       typename InputReaderType::Pointer likeReader=InputReaderType::New();
 436       likeReader->SetFileName(m_ArgsInfo.like_arg);
 437       likeReader->Update();
 438       resampler->SetSize( likeReader->GetOutput()->GetLargestPossibleRegion().GetSize() );
 439       resampler->SetOutputSpacing( likeReader->GetOutput()->GetSpacing() );
 440       resampler->SetOutputOrigin(  likeReader->GetOutput()->GetOrigin() );
 441     } else {
 442       //Size
 443       typename OutputImageType::SizeType outputSize;
 444       if (m_ArgsInfo.size_given) {
 445         for(unsigned int i=0; i< Dimension; i++)
 446           outputSize[i]=m_ArgsInfo.size_arg[i];
 447       } else outputSize=input->GetLargestPossibleRegion().GetSize();
 448       std::cout<<"Setting the size to "<<outputSize<<"..."<<std::endl;
 449
 450       //Spacing
 451       typename OutputImageType::SpacingType outputSpacing;
 452       if (m_ArgsInfo.spacing_given) {
 453         for(unsigned int i=0; i< Dimension; i++)
 454           outputSpacing[i]=m_ArgsInfo.spacing_arg[i];
 455       } else outputSpacing=input->GetSpacing();
 456       std::cout<<"Setting the spacing to "<<outputSpacing<<"..."<<std::endl;
 457
 458       //Origin
 459       typename OutputImageType::PointType outputOrigin;
 460       if (m_ArgsInfo.origin_given) {
 461         for(unsigned int i=0; i< Dimension; i++)
 462           outputOrigin[i]=m_ArgsInfo.origin_arg[i];
 463       } else outputOrigin=input->GetOrigin();
 464       std::cout<<"Setting the origin to "<<outputOrigin<<"..."<<std::endl;
 465
 466       // Set
 467       resampler->SetSize( outputSize );
 468       resampler->SetOutputSpacing( outputSpacing );
 469       resampler->SetOutputOrigin(  outputOrigin );
 470
 471     }
 472
 473     resampler->SetInput( input );
 474     resampler->SetTransform( affineTransform );
 475     resampler->SetInterpolator( genericInterpolator->GetInterpolatorPointer());
 476     resampler->SetDefaultPixelValue( static_cast<PixelType>(m_ArgsInfo.pad_arg) );
 477
 478     try {
 479       resampler->Update();
 480     } catch(itk::ExceptionObject) {
 481       std::cerr<<"Error resampling the image"<<std::endl;
 482     }
 483
 484     typename OutputImageType::Pointer output = resampler->GetOutput();
 485
 486     // Output
 487     typedef itk::ImageFileWriter<OutputImageType> WriterType;
 488     typename WriterType::Pointer writer = WriterType::New();
 489     writer->SetFileName(m_ArgsInfo.output_arg);
 490     writer->SetInput(output);
 491     writer->Update();
 492
 493   }
 494   //-------------------------------------------------------------------
 495
 496 } //end clitk
 497
 498 #endif //#define clitkAffineTransformGenericFilter_txx