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