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