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 <itkRecursiveGaussianImageFilter.h>
  26 #include "clitkElastix.h"
  27
  28 namespace clitk
  29 {
  30
  31   //-----------------------------------------------------------
  32   // Constructor
  33   //-----------------------------------------------------------
  34   template<class args_info_type>
  35   AffineTransformGenericFilter<args_info_type>::AffineTransformGenericFilter()
  36   {
  37     m_Verbose=false;
  38     m_InputFileName="";
  39   }
  40   //-------------------------------------------------------------------
  41
  42
  43   //-----------------------------------------------------------
  44   // Update
  45   //-----------------------------------------------------------
  46   template<class args_info_type>
  47   void AffineTransformGenericFilter<args_info_type>::Update()
  48   {
  49     // Read the Dimension and PixelType
  50     int Dimension, Components;
  51     std::string PixelType;
  52     ReadImageDimensionAndPixelType(m_InputFileName, Dimension, PixelType, Components);
  53
  54     // Call UpdateWithDim
  55     if(Dimension==2) UpdateWithDim<2>(PixelType, Components);
  56     else
  57       if(Dimension==3) UpdateWithDim<3>(PixelType, Components);
  58       else if (Dimension==4)UpdateWithDim<4>(PixelType, Components);
  59       else {
  60         std::cout<<"Error, Only for 2, 3 or 4  Dimensions!!!"<<std::endl ;
  61         return;
  62       }
  63   }
  64   //-------------------------------------------------------------------
  65
  66
  67   //-------------------------------------------------------------------
  68   // Update with the number of dimensions
  69   //-------------------------------------------------------------------
  70   template<class args_info_type>
  71   template<unsigned int Dimension>
  72   void
  73   AffineTransformGenericFilter<args_info_type>::UpdateWithDim(std::string PixelType, int Components)
  74   {
  75     if (m_Verbose) std::cout << "Image was detected to be "<<Dimension<<"D and "<<Components<<" component(s) of "<<  PixelType<<"..."<<std::endl;
  76
  77     if (Components==1) {
  78       if(PixelType == "short") {
  79         if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and signed short..." << std::endl;
  80         UpdateWithDimAndPixelType<Dimension, signed short>();
  81       }
  82       else if(PixelType == "unsigned_short"){
  83         if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and unsigned_short..." << std::endl;
  84         UpdateWithDimAndPixelType<Dimension, unsigned short>();
  85       }
  86
  87       else if (PixelType == "unsigned_char") {
  88         if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and unsigned_char..." << std::endl;
  89         UpdateWithDimAndPixelType<Dimension, unsigned char>();
  90       }
  91
  92       //     else if (PixelType == "char"){
  93       //       if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and signed_char..." << std::endl;
  94       //       UpdateWithDimAndPixelType<Dimension, signed char>();
  95       //     }
  96       else {
  97         if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and float..." << std::endl;
  98         UpdateWithDimAndPixelType<Dimension, float>();
  99       }
 100     }
 101
 102     else if (Components==3) {
 103       if (m_Verbose) std::cout  << "Launching transform in "<< Dimension <<"D and 3D float (DVF)" << std::endl;
 104       UpdateWithDimAndVectorType<Dimension, itk::Vector<float, Dimension> >();
 105     }
 106
 107     else std::cerr<<"Number of components is "<<Components<<", not supported!"<<std::endl;
 108
 109   }
 110   //-------------------------------------------------------------------
 111
 112
 113   //-------------------------------------------------------------------
 114   // Update with the number of dimensions and the pixeltype
 115   //-------------------------------------------------------------------
 116   template<class args_info_type>
 117   template <unsigned int Dimension, class  PixelType>
 118   void
 119   AffineTransformGenericFilter<args_info_type>::UpdateWithDimAndPixelType()
 120   {
 121
 122     // ImageTypes
 123     typedef itk::Image<PixelType, Dimension> InputImageType;
 124     typedef itk::Image<PixelType, Dimension> OutputImageType;
 125
 126     // Read the input
 127     typedef itk::ImageFileReader<InputImageType> InputReaderType;
 128     typename InputReaderType::Pointer reader = InputReaderType::New();
 129     reader->SetFileName( m_InputFileName);
 130     reader->Update();
 131     typename InputImageType::Pointer input= reader->GetOutput();
 132
 133     //Gaussian pre-filtering
 134     typename itk::Vector<double, Dimension> gaussianSigma;
 135     gaussianSigma.Fill(0);
 136     bool gaussianFilteringEnabled(false);
 137     bool autoGaussEnabled(false);
 138     if (m_ArgsInfo.autogauss_given) { // Gaussian filter auto
 139       autoGaussEnabled = m_ArgsInfo.autogauss_flag;
 140     }
 141     if (m_ArgsInfo.gauss_given) { // Gaussian filter set by user
 142       gaussianFilteringEnabled = true;
 143       if (m_ArgsInfo.gauss_given == 1)
 144       {
 145         for (unsigned int i=0; i<Dimension; i++)
 146         {
 147           gaussianSigma[i] = m_ArgsInfo.gauss_arg[0];
 148         }
 149       }
 150       else if (m_ArgsInfo.gauss_given == Dimension)
 151       {
 152         for (unsigned int i=0; i<Dimension; i++)
 153         {
 154           gaussianSigma[i] = m_ArgsInfo.gauss_arg[i];
 155         }
 156       }
 157       else
 158       {
 159         std::cerr << "Gaussian sigma dimension is incorrect" << std::endl;
 160         return;
 161       }
 162     }
 163
 164     //Filter
 165     typedef  itk::ResampleImageFilter< InputImageType,OutputImageType >  ResampleFilterType;
 166     typename ResampleFilterType::Pointer resampler = ResampleFilterType::New();
 167
 168     // Matrix
 169     typename itk::Matrix<double, Dimension+1, Dimension+1> matrix;
 170     if (m_ArgsInfo.rotate_given || m_ArgsInfo.translate_given)
 171       {
 172         if (m_ArgsInfo.matrix_given)
 173           {
 174             std::cerr << "You must use either rotate/translate or matrix options" << std::endl;
 175             return;
 176           }
 177         itk::Array<double> transformParameters(2 * Dimension);
 178         transformParameters.Fill(0.0);
 179         if (m_ArgsInfo.rotate_given)
 180           {
 181             if (Dimension == 2)
 182               transformParameters[0] = m_ArgsInfo.rotate_arg[0];
 183             else
 184               for (unsigned int i = 0; i < 3; i++)
 185                 transformParameters[i] = m_ArgsInfo.rotate_arg[i];
 186           }
 187         if (m_ArgsInfo.translate_given)
 188           {
 189             int pos = 3;
 190             if (Dimension == 2)
 191               pos = 1;
 192             for (unsigned int i = 0; i < Dimension && i < 3; i++)
 193               transformParameters[pos++] = m_ArgsInfo.translate_arg[i];
 194           }
 195         if (Dimension == 4)
 196           {
 197             matrix.SetIdentity();
 198             itk::Matrix<double, 4, 4> tmp = GetForwardAffineMatrix3D(transformParameters);
 199             for (unsigned int i = 0; i < 3; ++i)
 200               for (unsigned int j = 0; j < 3; ++j)
 201                 matrix[i][j] = tmp[i][j];
 202             for (unsigned int i = 0; i < 3; ++i)
 203               matrix[i][4] = tmp[i][3];
 204           }
 205         else
 206           matrix = GetForwardAffineMatrix<Dimension>(transformParameters);
 207       }
 208     else
 209       {
 210         if (m_ArgsInfo.matrix_given)
 211           {
 212             matrix= clitk::ReadMatrix<Dimension>(m_ArgsInfo.matrix_arg);
 213             if (m_Verbose) std::cout << "Reading the matrix..." << std::endl;
 214           }
 215         else {
 216           if (m_ArgsInfo.elastix_given) {
 217             std::string filename(m_ArgsInfo.elastix_arg);
 218             matrix = createMatrixFromElastixFile<Dimension>(filename, m_Verbose);
 219           }
 220           else
 221             matrix.SetIdentity();
 222         }
 223       }
 224     if (m_Verbose)
 225       std::cout << "Using the following matrix:" << std::endl
 226                 << matrix << std::endl;
 227     typename itk::Matrix<double, Dimension, Dimension> rotationMatrix = clitk::GetRotationalPartMatrix(matrix);
 228     typename itk::Vector<double, Dimension> translationPart = clitk::GetTranslationPartMatrix(matrix);
 229
 230     // Transform
 231     typedef itk::AffineTransform<double, Dimension> AffineTransformType;
 232     typename AffineTransformType::Pointer affineTransform=AffineTransformType::New();
 233     affineTransform->SetMatrix(rotationMatrix);
 234     affineTransform->SetTranslation(translationPart);
 235
 236     // Interp
 237     typedef clitk::GenericInterpolator<args_info_type, InputImageType, double> GenericInterpolatorType;
 238     typename GenericInterpolatorType::Pointer genericInterpolator=GenericInterpolatorType::New();
 239     genericInterpolator->SetArgsInfo(m_ArgsInfo);
 240
 241     // Properties
 242     if (m_ArgsInfo.like_given) {
 243       typename InputReaderType::Pointer likeReader=InputReaderType::New();
 244       likeReader->SetFileName(m_ArgsInfo.like_arg);
 245       likeReader->Update();
 246       resampler->SetOutputParametersFromImage(likeReader->GetOutput());
 247       resampler->SetOutputDirection(likeReader->GetOutput()->GetDirection());
 248       if (autoGaussEnabled) { // Automated sigma when downsample
 249         for(unsigned int i=0; i<Dimension; i++) {
 250           if (likeReader->GetOutput()->GetSpacing()[i] > input->GetSpacing()[i]) { // downsample
 251             gaussianSigma[i] = 0.5*likeReader->GetOutput()->GetSpacing()[i];// / inputSpacing[i]);
 252           }
 253           else gaussianSigma[i] = 0; // will be ignore after
 254         }
 255       }
 256     } else if(m_ArgsInfo.transform_grid_flag) {
 257       typename itk::Matrix<double, Dimension+1, Dimension+1> invMatrix( matrix.GetInverse() );
 258       typename itk::Matrix<double, Dimension, Dimension> invRotMatrix( clitk::GetRotationalPartMatrix(invMatrix) );
 259       typename itk::Vector<double,Dimension> invTrans =  clitk::GetTranslationPartMatrix(invMatrix);
 260
 261       // Display warning
 262       if (m_ArgsInfo.spacing_given)
 263         std::cout << "Warning --spacing ignored (because --transform_grid_flag)" << std::endl;
 264       if (m_ArgsInfo.origin_given)
 265         std::cout << "Warning --origin ignored (because --transform_grid_flag)" << std::endl;
 266
 267       // Spacing is influenced by affine transform matrix and input direction
 268       typename InputImageType::SpacingType outputSpacing;
 269       outputSpacing = invRotMatrix *
 270         input->GetDirection() *
 271         input->GetSpacing();
 272       if (autoGaussEnabled) { // Automated sigma when downsample
 273         for(unsigned int i=0; i<Dimension; i++) {
 274           if (outputSpacing[i] > input->GetSpacing()[i]) { // downsample
 275             gaussianSigma[i] = 0.5*outputSpacing[i];// / inputSpacing[i]);
 276           }
 277           else gaussianSigma[i] = 0; // will be ignore after
 278         }
 279       }
 280
 281       // Origin is influenced by translation but not by input direction
 282       typename InputImageType::PointType outputOrigin;
 283       outputOrigin = invRotMatrix *
 284         input->GetOrigin() +
 285         invTrans;
 286
 287       // Size is influenced by affine transform matrix and input direction
 288       // Size is converted to double, transformed and converted back to size type.
 289       vnl_vector<double> vnlOutputSize(Dimension);
 290       for(unsigned int i=0; i< Dimension; i++) {
 291         vnlOutputSize[i] = input->GetLargestPossibleRegion().GetSize()[i];
 292       }
 293       vnlOutputSize = invRotMatrix *
 294         input->GetDirection().GetVnlMatrix() *
 295         vnlOutputSize;
 296       typename OutputImageType::SizeType outputSize;
 297       for(unsigned int i=0; i< Dimension; i++) {
 298         // If the size is negative, we have a flip and we must modify
 299         // the origin and the spacing accordingly.
 300         if(vnlOutputSize[i]<0.) {
 301           vnlOutputSize[i] *= -1.;
 302           outputOrigin[i] = outputOrigin[i] + outputSpacing[i] * (vnlOutputSize[i]-1);
 303           outputSpacing[i] *= -1.;
 304         }
 305         outputSize[i] = lrint(vnlOutputSize[i]);
 306       }
 307       resampler->SetSize( outputSize );
 308       resampler->SetOutputSpacing( outputSpacing );
 309       resampler->SetOutputOrigin( outputOrigin );
 310     } else {
 311       //Size
 312       typename OutputImageType::SizeType outputSize;
 313       if (m_ArgsInfo.size_given) {
 314         for(unsigned int i=0; i< Dimension; i++)
 315           outputSize[i]=m_ArgsInfo.size_arg[i];
 316       } else outputSize=input->GetLargestPossibleRegion().GetSize();
 317
 318       //Spacing
 319       typename OutputImageType::SpacingType outputSpacing;
 320       if (m_ArgsInfo.spacing_given) {
 321         for(unsigned int i=0; i< Dimension; i++)
 322           outputSpacing[i]=m_ArgsInfo.spacing_arg[i];
 323       } else outputSpacing=input->GetSpacing();
 324       if (autoGaussEnabled) { // Automated sigma when downsample
 325         for(unsigned int i=0; i<Dimension; i++) {
 326           if (outputSpacing[i] > input->GetSpacing()[i]) { // downsample
 327             gaussianSigma[i] = 0.5*outputSpacing[i];// / inputSpacing[i]);
 328           }
 329           else gaussianSigma[i] = 0; // will be ignore after
 330         }
 331       }
 332
 333       //Origin
 334       typename OutputImageType::PointType outputOrigin;
 335       if (m_ArgsInfo.origin_given) {
 336         for(unsigned int i=0; i< Dimension; i++)
 337           outputOrigin[i]=m_ArgsInfo.origin_arg[i];
 338       } else outputOrigin=input->GetOrigin();
 339
 340       //Direction
 341       typename OutputImageType::DirectionType outputDirection;
 342       if (m_ArgsInfo.direction_given) {
 343         for(unsigned int j=0; j< Dimension; j++)
 344             for(unsigned int i=0; i< Dimension; i++)
 345                 outputDirection[j][i]=m_ArgsInfo.direction_arg[i+Dimension*j];
 346       } else outputDirection=input->GetDirection();
 347
 348       // Set
 349       resampler->SetSize( outputSize );
 350       resampler->SetOutputSpacing( outputSpacing );
 351       resampler->SetOutputOrigin(  outputOrigin );
 352       resampler->SetOutputDirection( outputDirection );
 353
 354     }
 355
 356     if (m_ArgsInfo.spacinglike_given) {
 357       typename InputReaderType::Pointer likeReader=InputReaderType::New();
 358       likeReader->SetFileName(m_ArgsInfo.spacinglike_arg);
 359       likeReader->Update();
 360
 361       // set the support like the image
 362       if (m_ArgsInfo.like_given) {
 363         typename OutputImageType::SizeType outputSize;
 364         outputSize[0] = ceil(resampler->GetSize()[0]*resampler->GetOutputSpacing()[0]
 365                              /likeReader->GetOutput()->GetSpacing()[0]);
 366         outputSize[1] = ceil(resampler->GetSize()[1]*resampler->GetOutputSpacing()[1]
 367                              /likeReader->GetOutput()->GetSpacing()[1]);
 368         outputSize[2] = ceil(resampler->GetSize()[2]*resampler->GetOutputSpacing()[2]
 369                              /likeReader->GetOutput()->GetSpacing()[2]);
 370         if (m_ArgsInfo.verbose_flag) {
 371           std::cout << "Compute the number of pixels such as the support is like " << m_ArgsInfo.like_arg << std::endl;
 372         }
 373         resampler->SetSize( outputSize );
 374       }
 375
 376       resampler->SetOutputSpacing( likeReader->GetOutput()->GetSpacing() );
 377     }
 378
 379     if (m_ArgsInfo.verbose_flag) {
 380       std::cout << "Setting the output size to " << resampler->GetSize() << "..." << std::endl;
 381       std::cout << "Setting the output spacing to " << resampler->GetOutputSpacing() << "..." << std::endl;
 382       std::cout << "Setting the output origin to " << resampler->GetOutputOrigin() << "..." << std::endl;
 383       std::cout << "Setting the output direction to " << resampler->GetOutputDirection() << "..." << std::endl;
 384     }
 385
 386     typedef itk::RecursiveGaussianImageFilter<InputImageType, InputImageType> GaussianFilterType;
 387     std::vector<typename GaussianFilterType::Pointer> gaussianFilters;
 388     if (gaussianFilteringEnabled || autoGaussEnabled) {
 389       for(unsigned int i=0; i<Dimension; i++) {
 390         if (gaussianSigma[i] != 0) {
 391           gaussianFilters.push_back(GaussianFilterType::New());
 392           gaussianFilters[i]->SetDirection(i);
 393           gaussianFilters[i]->SetOrder(GaussianFilterType::ZeroOrder);
 394           gaussianFilters[i]->SetNormalizeAcrossScale(false);
 395           gaussianFilters[i]->SetSigma(gaussianSigma[i]); // in millimeter !
 396           if (gaussianFilters.size() == 1) { // first
 397             gaussianFilters[0]->SetInput(input);
 398           } else {
 399             gaussianFilters[i]->SetInput(gaussianFilters[i-1]->GetOutput());
 400           }
 401         }
 402       }
 403       if (gaussianFilters.size() > 0) {
 404         resampler->SetInput(gaussianFilters[gaussianFilters.size()-1]->GetOutput());
 405       } else resampler->SetInput(input);
 406     } else resampler->SetInput(input);
 407
 408     resampler->SetTransform( affineTransform );
 409     resampler->SetInterpolator( genericInterpolator->GetInterpolatorPointer());
 410     resampler->SetDefaultPixelValue( static_cast<PixelType>(m_ArgsInfo.pad_arg) );
 411
 412     try {
 413       resampler->Update();
 414     } catch(itk::ExceptionObject) {
 415       std::cerr<<"Error resampling the image"<<std::endl;
 416     }
 417
 418     typename OutputImageType::Pointer output = resampler->GetOutput();
 419
 420     // Output
 421     typedef itk::ImageFileWriter<OutputImageType> WriterType;
 422     typename WriterType::Pointer writer = WriterType::New();
 423     writer->SetFileName(m_ArgsInfo.output_arg);
 424     writer->SetInput(output);
 425     writer->Update();
 426
 427   }
 428   //-------------------------------------------------------------------
 429
 430
 431   //-------------------------------------------------------------------
 432   // Update with the number of dimensions and the pixeltype (components)
 433   //-------------------------------------------------------------------
 434   template<class args_info_type>
 435   template<unsigned int Dimension, class PixelType>
 436   void AffineTransformGenericFilter<args_info_type>::UpdateWithDimAndVectorType()
 437   {
 438     // ImageTypes
 439     typedef itk::Image<PixelType, Dimension> InputImageType;
 440     typedef itk::Image<PixelType, Dimension> OutputImageType;
 441
 442     // Read the input
 443     typedef itk::ImageFileReader<InputImageType> InputReaderType;
 444     typename InputReaderType::Pointer reader = InputReaderType::New();
 445     reader->SetFileName( m_InputFileName);
 446     reader->Update();
 447     typename InputImageType::Pointer input= reader->GetOutput();
 448
 449     //Gaussian pre-filtering
 450     typename itk::Vector<double, Dimension> gaussianSigma;
 451     gaussianSigma.Fill(0);
 452     bool gaussianFilteringEnabled(false);
 453     bool autoGaussEnabled(false);
 454     if (m_ArgsInfo.autogauss_given) { // Gaussian filter auto
 455       autoGaussEnabled = m_ArgsInfo.autogauss_flag;
 456     }
 457     if (m_ArgsInfo.gauss_given) { // Gaussian filter set by user
 458       gaussianFilteringEnabled = true;
 459       if (m_ArgsInfo.gauss_given == 1)
 460       {
 461         for (unsigned int i=0; i<Dimension; i++)
 462         {
 463           gaussianSigma[i] = m_ArgsInfo.gauss_arg[0];
 464         }
 465       }
 466       else if (m_ArgsInfo.gauss_given == Dimension)
 467       {
 468         for (unsigned int i=0; i<Dimension; i++)
 469         {
 470           gaussianSigma[i] = m_ArgsInfo.gauss_arg[i];
 471         }
 472       }
 473       else
 474       {
 475         std::cerr << "Gaussian sigma dimension is incorrect" << std::endl;
 476         return;
 477       }
 478     }
 479
 480     //Filter
 481     typedef  itk::VectorResampleImageFilter< InputImageType,OutputImageType, double >  ResampleFilterType;
 482     typename ResampleFilterType::Pointer resampler = ResampleFilterType::New();
 483
 484     // Matrix
 485     typename itk::Matrix<double, Dimension+1, Dimension+1> matrix;
 486     if (m_ArgsInfo.rotate_given || m_ArgsInfo.translate_given)
 487       {
 488         if (m_ArgsInfo.matrix_given)
 489           {
 490             std::cerr << "You must use either rotate/translate or matrix options" << std::endl;
 491             return;
 492           }
 493         itk::Array<double> transformParameters(2 * Dimension);
 494         transformParameters.Fill(0.0);
 495         if (m_ArgsInfo.rotate_given)
 496           {
 497             if (Dimension == 2)
 498               transformParameters[0] = m_ArgsInfo.rotate_arg[0];
 499             else
 500               for (unsigned int i = 0; i < 3; i++)
 501                 transformParameters[i] = m_ArgsInfo.rotate_arg[i];
 502           }
 503         if (m_ArgsInfo.translate_given)
 504           {
 505             int pos = 3;
 506             if (Dimension == 2)
 507               pos = 1;
 508             for (unsigned int i = 0; i < Dimension && i < 3; i++)
 509               transformParameters[pos++] = m_ArgsInfo.translate_arg[i];
 510           }
 511         if (Dimension == 4)
 512           {
 513             matrix.SetIdentity();
 514             itk::Matrix<double, 4, 4> tmp = GetForwardAffineMatrix3D(transformParameters);
 515             for (unsigned int i = 0; i < 3; ++i)
 516               for (unsigned int j = 0; j < 3; ++j)
 517                 matrix[i][j] = tmp[i][j];
 518             for (unsigned int i = 0; i < 3; ++i)
 519               matrix[i][4] = tmp[i][3];
 520           }
 521         else
 522           matrix = GetForwardAffineMatrix<Dimension>(transformParameters);
 523       }
 524     else
 525       {
 526         if (m_ArgsInfo.matrix_given)
 527           {
 528             matrix= clitk::ReadMatrix<Dimension>(m_ArgsInfo.matrix_arg);
 529             if (m_Verbose) std::cout << "Reading the matrix..." << std::endl;
 530           }
 531         else
 532           matrix.SetIdentity();
 533       }
 534     if (m_Verbose)
 535       std::cout << "Using the following matrix:" << std::endl
 536                 << matrix << std::endl;
 537     typename itk::Matrix<double, Dimension, Dimension> rotationMatrix = clitk::GetRotationalPartMatrix(matrix);
 538     typename itk::Vector<double, Dimension> translationPart = clitk::GetTranslationPartMatrix(matrix);
 539
 540     // Transform
 541     typedef itk::AffineTransform<double, Dimension> AffineTransformType;
 542     typename AffineTransformType::Pointer affineTransform=AffineTransformType::New();
 543     affineTransform->SetMatrix(rotationMatrix);
 544     affineTransform->SetTranslation(translationPart);
 545
 546     // Interp
 547     typedef clitk::GenericVectorInterpolator<args_info_type, InputImageType, double> GenericInterpolatorType;
 548     typename GenericInterpolatorType::Pointer genericInterpolator=GenericInterpolatorType::New();
 549     genericInterpolator->SetArgsInfo(m_ArgsInfo);
 550
 551     // Properties
 552     if (m_ArgsInfo.like_given) {
 553       typename InputReaderType::Pointer likeReader=InputReaderType::New();
 554       likeReader->SetFileName(m_ArgsInfo.like_arg);
 555       likeReader->Update();
 556       resampler->SetSize( likeReader->GetOutput()->GetLargestPossibleRegion().GetSize() );
 557       resampler->SetOutputSpacing( likeReader->GetOutput()->GetSpacing() );
 558       resampler->SetOutputOrigin(  likeReader->GetOutput()->GetOrigin() );
 559       resampler->SetOutputDirection( likeReader->GetOutput()->GetDirection() );
 560       if (autoGaussEnabled) { // Automated sigma when downsample
 561         for(unsigned int i=0; i<Dimension; i++) {
 562           if (likeReader->GetOutput()->GetSpacing()[i] > input->GetSpacing()[i]) { // downsample
 563             gaussianSigma[i] = 0.5*likeReader->GetOutput()->GetSpacing()[i];// / inputSpacing[i]);
 564           }
 565           else gaussianSigma[i] = 0; // will be ignore after
 566         }
 567       }
 568     } else {
 569       //Size
 570       typename OutputImageType::SizeType outputSize;
 571       if (m_ArgsInfo.size_given) {
 572         for(unsigned int i=0; i< Dimension; i++)
 573           outputSize[i]=m_ArgsInfo.size_arg[i];
 574       } else outputSize=input->GetLargestPossibleRegion().GetSize();
 575       std::cout<<"Setting the size to "<<outputSize<<"..."<<std::endl;
 576
 577       //Spacing
 578       typename OutputImageType::SpacingType outputSpacing;
 579       if (m_ArgsInfo.spacing_given) {
 580         for(unsigned int i=0; i< Dimension; i++)
 581           outputSpacing[i]=m_ArgsInfo.spacing_arg[i];
 582       } else outputSpacing=input->GetSpacing();
 583       if (autoGaussEnabled) { // Automated sigma when downsample
 584         for(unsigned int i=0; i<Dimension; i++) {
 585           if (outputSpacing[i] > input->GetSpacing()[i]) { // downsample
 586             gaussianSigma[i] = 0.5*outputSpacing[i];// / inputSpacing[i]);
 587           }
 588           else gaussianSigma[i] = 0; // will be ignore after
 589         }
 590       }
 591       std::cout<<"Setting the spacing to "<<outputSpacing<<"..."<<std::endl;
 592
 593       //Origin
 594       typename OutputImageType::PointType outputOrigin;
 595       if (m_ArgsInfo.origin_given) {
 596         for(unsigned int i=0; i< Dimension; i++)
 597           outputOrigin[i]=m_ArgsInfo.origin_arg[i];
 598       } else outputOrigin=input->GetOrigin();
 599       std::cout<<"Setting the origin to "<<outputOrigin<<"..."<<std::endl;
 600
 601       //Direction
 602       typename OutputImageType::DirectionType outputDirection;
 603       if (m_ArgsInfo.direction_given) {
 604         for(unsigned int j=0; j< Dimension; j++)
 605             for(unsigned int i=0; i< Dimension; i++)
 606                 outputDirection[j][i]=m_ArgsInfo.direction_arg[i+Dimension*j];
 607       } else outputDirection=input->GetDirection();
 608       std::cout<<"Setting the direction to "<<outputDirection<<"..."<<std::endl;
 609
 610       // Set
 611       resampler->SetSize( outputSize );
 612       resampler->SetOutputSpacing( outputSpacing );
 613       resampler->SetOutputOrigin(  outputOrigin );
 614       resampler->SetOutputDirection( outputDirection );
 615
 616     }
 617
 618     typedef itk::RecursiveGaussianImageFilter<InputImageType, InputImageType> GaussianFilterType;
 619     std::vector<typename GaussianFilterType::Pointer> gaussianFilters;
 620     if (gaussianFilteringEnabled || autoGaussEnabled) {
 621       for(unsigned int i=0; i<Dimension; i++) {
 622         if (gaussianSigma[i] != 0) {
 623           gaussianFilters.push_back(GaussianFilterType::New());
 624           gaussianFilters[i]->SetDirection(i);
 625           gaussianFilters[i]->SetOrder(GaussianFilterType::ZeroOrder);
 626           gaussianFilters[i]->SetNormalizeAcrossScale(false);
 627           gaussianFilters[i]->SetSigma(gaussianSigma[i]); // in millimeter !
 628           if (gaussianFilters.size() == 1) { // first
 629             gaussianFilters[0]->SetInput(input);
 630           } else {
 631             gaussianFilters[i]->SetInput(gaussianFilters[i-1]->GetOutput());
 632           }
 633         }
 634       }
 635       if (gaussianFilters.size() > 0) {
 636         resampler->SetInput(gaussianFilters[gaussianFilters.size()-1]->GetOutput());
 637       } else resampler->SetInput(input);
 638     } else resampler->SetInput(input);
 639
 640     resampler->SetInput( input );
 641     resampler->SetTransform( affineTransform );
 642     resampler->SetInterpolator( genericInterpolator->GetInterpolatorPointer());
 643     resampler->SetDefaultPixelValue( static_cast<PixelType>(m_ArgsInfo.pad_arg) );
 644
 645     try {
 646       resampler->Update();
 647     } catch(itk::ExceptionObject) {
 648       std::cerr<<"Error resampling the image"<<std::endl;
 649     }
 650
 651     typename OutputImageType::Pointer output = resampler->GetOutput();
 652
 653     // Output
 654     typedef itk::ImageFileWriter<OutputImageType> WriterType;
 655     typename WriterType::Pointer writer = WriterType::New();
 656     writer->SetFileName(m_ArgsInfo.output_arg);
 657     writer->SetInput(output);
 658     writer->Update();
 659
 660   }
 661   //-------------------------------------------------------------------
 662
 663 } //end clitk
 664
 665 #endif //#define clitkAffineTransformGenericFilter_txx