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