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             std::vector<std::string> s;
 185             for(uint i=0; i<m_ArgsInfo.elastix_given; i++) s.push_back(m_ArgsInfo.elastix_arg[i]);
 186             matrix = createMatrixFromElastixFile<Dimension,PixelType>(s);
 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
 496   //-------------------------------------------------------------------
 497   template<class args_info_type>
 498   template<unsigned int Dimension, class PixelType>
 499   typename itk::Matrix<double, Dimension+1, Dimension+1>
 500                                                            AffineTransformGenericFilter<args_info_type>::createMatrixFromElastixFile(std::vector<std::string> & filename)
 501   {
 502     if (Dimension != 3) {
 503       FATAL("Only 3D yet" << std::endl);
 504     }
 505     typename itk::Matrix<double, Dimension+1, Dimension+1> matrix;
 506
 507     itk::CenteredEuler3DTransform<double>::Pointer mat = itk::CenteredEuler3DTransform<double>::New();
 508     itk::CenteredEuler3DTransform<double>::Pointer previous;
 509     for(uint j=0; j<filename.size(); j++) {
 510
 511       // Open file
 512       if (m_Verbose) std::cout << "Read elastix parameters in " << filename[j] << std::endl;
 513       std::ifstream is;
 514       clitk::openFileForReading(is, filename[j]);
 515
 516       // Check Transform
 517       std::string s;
 518       bool b = GetElastixValueFromTag(is, "Transform ", s);
 519       if (!b) {
 520         FATAL("Error must read 'Transform' in " << filename[j] << std::endl);
 521       }
 522       if (s != "EulerTransform") {
 523         FATAL("Sorry only 'EulerTransform'" << std::endl);
 524       }
 525
 526       // FIXME check
 527       //    (InitialTransformParametersFilename[j] "NoInitialTransform")
 528
 529       // Get CenterOfRotationPoint
 530       GetElastixValueFromTag(is, "CenterOfRotationPoint ", s); // space is needed
 531       if (!b) {
 532         FATAL("Error must read 'CenterOfRotationPoint' in " << filename[j] << std::endl);
 533       }
 534       std::vector<std::string> cor;
 535       GetValuesFromValue(s, cor);
 536
 537       // Get Transformparameters
 538       GetElastixValueFromTag(is, "TransformParameters ", s); // space is needed
 539       if (!b) {
 540         FATAL("Error must read 'TransformParameters' in " << filename[j] << std::endl);
 541       }
 542       std::vector<std::string> results;
 543       GetValuesFromValue(s, results);
 544
 545       // construct a stream from the string
 546       itk::CenteredEuler3DTransform<double>::ParametersType p;
 547       p.SetSize(9);
 548       for(uint i=0; i<3; i++)
 549         p[i] = atof(results[i].c_str()); // Rotation
 550       for(uint i=0; i<3; i++)
 551         p[i+3] = atof(cor[i].c_str()); // Centre of rotation
 552       for(uint i=0; i<3; i++)
 553         p[i+6] = atof(results[i+3].c_str()); // Translation
 554       mat->SetParameters(p);
 555
 556       if (m_Verbose) {
 557         std::cout << "Rotation      (deg) : " << rad2deg(p[0]) << " " << rad2deg(p[1]) << " " << rad2deg(p[2]) << std::endl;
 558         std::cout << "Center of rot (phy) : " << p[3] << " " << p[4] << " " << p[5] << std::endl;
 559         std::cout << "Translation   (phy) : " << p[6] << " " << p[7] << " " << p[8] << std::endl;
 560       }
 561
 562       // Compose with previous if needed
 563       if (j!=0) {
 564         mat->Compose(previous);
 565         if (m_Verbose) {
 566           std::cout << "Composed rotation      (deg) : " << rad2deg(mat->GetAngleX()) << " " << rad2deg(mat->GetAngleY()) << " " << rad2deg(mat->GetAngleZ()) << std::endl;
 567           std::cout << "Composed center of rot (phy) : " << mat->GetCenter() << std::endl;
 568           std::cout << "Compsoed translation   (phy) : " << mat->GetTranslation() << std::endl;
 569         }
 570       }
 571       // previous = mat->Clone(); // ITK4
 572       previous = itk::CenteredEuler3DTransform<double>::New();
 573       previous->SetParameters(mat->GetParameters());
 574     }
 575
 576     mat = previous;
 577     for(uint i=0; i<3; i++)
 578       for(uint j=0; j<3; j++)
 579         matrix[i][j] = mat->GetMatrix()[i][j];
 580     // Offset is -Rc + t + c
 581     matrix[0][3] = mat->GetOffset()[0];
 582     matrix[1][3] = mat->GetOffset()[1];
 583     matrix[2][3] = mat->GetOffset()[2];
 584     matrix[3][3] = 1;
 585
 586     return matrix;
 587   }
 588
 589   //-------------------------------------------------------------------
 590   template<class args_info_type>
 591   bool
 592   AffineTransformGenericFilter<args_info_type>::GetElastixValueFromTag(std::ifstream & is,
 593                                                                        std::string tag,
 594                                                                        std::string & value)
 595   {
 596     std::string line;
 597     is.seekg (0, is.beg);
 598     while(std::getline(is, line))   {
 599       unsigned pos = line.find(tag);
 600       if (pos<line.size()) {
 601         value=line.substr(pos+tag.size(),line.size()-2);// remove last ')'
 602         value.erase (std::remove (value.begin(), value.end(), '"'), value.end());
 603         value.erase (std::remove (value.begin(), value.end(), ')'), value.end());
 604         return true;
 605       }
 606     }
 607     return false;
 608   }
 609   //-------------------------------------------------------------------
 610
 611
 612   //-------------------------------------------------------------------
 613   template<class args_info_type>
 614   void
 615   AffineTransformGenericFilter<args_info_type>::GetValuesFromValue(const std::string & s,
 616                                                                    std::vector<std::string> & values)
 617   {
 618     std::stringstream strstr(s);
 619     std::istream_iterator<std::string> it(strstr);
 620     std::istream_iterator<std::string> end;
 621     std::vector<std::string> results(it, end);
 622     values.clear();
 623     values.resize(results.size());
 624     for(uint i=0; i<results.size(); i++) values[i] = results[i];
 625   }
 626   //-------------------------------------------------------------------
 627
 628
 629 } //end clitk
 630
 631 #endif //#define clitkAffineTransformGenericFilter_txx