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