tools/clitkAffineTransformGenericFilter.txx

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