common/clitkTransformUtilities.h

   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
  19 #ifndef CLITKTRANSFORMUTILITIES_H
  20 #define CLITKTRANSFORMUTILITIES_H
  21 #include "itkMatrix.h"
  22 #include "itkArray.h"
  23 #include "itkPoint.h"
  24 #include "clitkImageCommon.h"
  25 #include "clitkCommon.h"
  26 #include <vtkMatrix4x4.h>
  27 #include <vtkSmartPointer.h>
  28
  29 namespace clitk
  30 {
  31   //============================================================================
  32   //Declarations
  33   //============================================================================
  34   itk::Matrix<double, 3, 3> GetForwardAffineMatrix2D(itk::Array<double> transformParameters);
  35   itk::Matrix<double, 4, 4> GetForwardAffineMatrix3D(itk::Array<double> transformParameters);
  36   template <unsigned int Dimension > itk::Matrix<double, Dimension+1, Dimension+1> GetForwardAffineMatrix(itk::Array<double> transformParameters);
  37
  38   itk::Matrix<double, 3, 3> GetBackwardAffineMatrix2D(itk::Array<double> transformParameters);
  39   itk::Matrix<double, 4, 4> GetBackwardAffineMatrix3D(itk::Array<double> transformParameters);
  40   template <unsigned int Dimension > itk::Matrix<double, Dimension+1, Dimension+1> GetBackwardAffineMatrix(itk::Array<double> transformParameters);
  41
  42   itk::Matrix<double, 3, 3> GetRotationMatrix3D(itk::Array<double> rotationParameters);
  43   itk::Matrix<double, 2, 2> GetRotationMatrix2D(itk::Array<double> rotationParameters);
  44   template <unsigned int Dimension> itk::Matrix<double, Dimension, Dimension> GetRotationMatrix(itk::Array<double> rotationParameters);
  45   itk::Point<double, 3> GetRotatedPoint3D(itk::Array<double> rotationParameters, itk::Point<double, 3> input);
  46   itk::Matrix<double, 4, 4> GetCenteredRotationMatrix3D(itk::Array<double> rotationParameters,itk::Point<double,3> centerOfRotation);
  47   //   itk::Matrix<double, 4, 4> GetComposedMatrix3D(itk::Matrix<double, 4, 4> firstTransform, itk::Matrix<double, 4, 4> secondTransform);
  48
  49   itk::Matrix<double, 5, 5> ReadMatrix4D(std::string fileName);
  50   itk::Matrix<double, 4, 4> ReadMatrix3D(std::string fileName);
  51   itk::Matrix<double, 3, 3> ReadMatrix2D(std::string fileName);
  52   template <unsigned int Dimension > itk::Matrix<double, Dimension+1 , Dimension+1> ReadMatrix(std::string fileName);
  53
  54   itk::Matrix<double, 3, 3> GetRotationalPartMatrix3D(itk::Matrix<double, 4, 4> input);
  55   itk::Matrix<double, 3, 3> GetRotationalPartMatrix(itk::Matrix<double, 4, 4> input);
  56   itk::Matrix<double, 2, 2> GetRotationalPartMatrix2D(itk::Matrix<double, 3, 3> input);
  57   itk::Matrix<double, 2, 2> GetRotationalPartMatrix(itk::Matrix<double, 3, 3> input);
  58
  59   itk::Vector<double,3> GetTranslationPartMatrix3D(itk::Matrix<double, 4, 4> input);
  60   itk::Vector<double,3> GetTranslationPartMatrix(itk::Matrix<double, 4, 4> input);
  61   itk::Vector<double,2> GetTranslationPartMatrix2D(itk::Matrix<double, 3, 3> input);
  62   itk::Vector<double,2> GetTranslationPartMatrix(itk::Matrix<double, 3, 3> input);
  63
  64
  65   //============================================================================
  66   //Inline functions definition in header file, otherwise linker errors
  67   //============================================================================
  68
  69   //========================================================================================
  70   inline  itk::Matrix<double, 3, 3> GetForwardAffineMatrix2D(itk::Array<double> transformParameters)
  71   {
  72     itk::Matrix<double, 3, 3> matrix;
  73     //rotation part
  74     matrix[0][0]=cos(transformParameters[0]);
  75     matrix[0][1]=-sin(transformParameters[0]);
  76     matrix[1][0]=sin(transformParameters[0]);
  77     matrix[1][1]=cos(transformParameters[0]);
  78     //translation part
  79     matrix[0][2]=transformParameters[1];
  80     matrix[1][2]=transformParameters[2];
  81     //homogenize
  82     matrix[2][0]=0.;
  83     matrix[2][1]=0.;
  84     matrix[2][2]=1.;
  85     return matrix;
  86   }
  87
  88   inline  itk::Matrix<double, 4, 4> GetForwardAffineMatrix3D(itk::Array<double> transformParameters)
  89   {
  90     itk::Matrix<double, 4, 4> matrix;
  91     //rotational part
  92     matrix[0][0]= cos(transformParameters[1])*cos(transformParameters[2]);
  93     matrix[0][1]= sin(transformParameters[0])*sin(transformParameters[1])*cos(transformParameters[2])+ sin(transformParameters[2])*cos(transformParameters[0]);
  94     matrix[0][2]= -cos(transformParameters[0])*sin(transformParameters[1])*cos(transformParameters[2])+sin(transformParameters[0])*sin(transformParameters[2]);
  95     matrix[1][0]= -cos(transformParameters[1])*sin(transformParameters[2]);
  96     matrix[1][1]= -sin(transformParameters[0])*sin(transformParameters[1])*sin(transformParameters[2])+cos(transformParameters[0])*cos(transformParameters[2]);
  97     matrix[1][2]= cos(transformParameters[0])*sin(transformParameters[1])*sin(transformParameters[2])+sin(transformParameters[0])*cos(transformParameters[2]);
  98     matrix[2][0]= sin(transformParameters[1]);
  99     matrix[2][1]= -sin(transformParameters[0])*cos(transformParameters[1]);
 100     matrix[2][2]= cos(transformParameters[0])*cos(transformParameters[1]);
 101     //translational part
 102     matrix[0][3]=transformParameters[3];
 103     matrix[1][3]=transformParameters[4];
 104     matrix[2][3]=transformParameters[5];
 105     //homogenize
 106     matrix[3][0]=0.;
 107     matrix[3][1]=0.;
 108     matrix[3][2]=0.;
 109     matrix[3][3]=1.;
 110     return matrix;
 111   }
 112
 113
 114   inline  itk::Matrix<double, 3, 3> GetBackwardAffineMatrix2D(itk::Array<double> transformParameters)
 115   {
 116     itk::Matrix<double, 3, 3> matrix;
 117     //rotation part
 118     matrix[0][0]=cos(transformParameters[0]);
 119     matrix[0][1]=sin(transformParameters[0]);
 120     matrix[1][0]=-sin(transformParameters[0]);
 121     matrix[1][1]=cos(transformParameters[0]);
 122     //translation part
 123     matrix[0][2]=transformParameters[1];
 124     matrix[1][2]=transformParameters[2];
 125     //homogenize
 126     matrix[2][0]=0.;
 127     matrix[2][1]=0.;
 128     matrix[2][2]=1.;
 129     return matrix;
 130   }
 131
 132
 133   inline  itk::Matrix<double, 4, 4> GetBackwardAffineMatrix3D(itk::Array<double> transformParameters)
 134   {
 135     itk::Matrix<double, 4, 4> matrix;
 136     //rotational part
 137     matrix[0][0]= cos(transformParameters[1])*cos(transformParameters[2]);
 138     matrix[0][1]= sin(transformParameters[0])*sin(transformParameters[1])*cos(transformParameters[2])- sin(transformParameters[2])*cos(transformParameters[0]);
 139     matrix[0][2]= cos(transformParameters[0])*sin(transformParameters[1])*cos(transformParameters[2])+sin(transformParameters[0])*sin(transformParameters[2]);
 140     matrix[1][0]= cos(transformParameters[1])*sin(transformParameters[2]);
 141     matrix[1][1]= sin(transformParameters[0])*sin(transformParameters[1])*sin(transformParameters[2])+cos(transformParameters[0])*cos(transformParameters[2]);
 142     matrix[1][2]= cos(transformParameters[0])*sin(transformParameters[1])*sin(transformParameters[2])-sin(transformParameters[0])*cos(transformParameters[2]);
 143     matrix[2][0]= -sin(transformParameters[1]);
 144     matrix[2][1]= sin(transformParameters[0])*cos(transformParameters[1]);
 145     matrix[2][2]= cos(transformParameters[0])*cos(transformParameters[1]);
 146     //translational part
 147     matrix[0][3]=transformParameters[3];
 148     matrix[1][3]=transformParameters[4];
 149     matrix[2][3]=transformParameters[5];
 150     //homogenize
 151     matrix[3][0]=0.;
 152     matrix[3][1]=0.;
 153     matrix[3][2]=0.;
 154     matrix[3][3]=1.;
 155     return matrix;
 156   }
 157
 158   inline itk::Matrix<double, 3, 3> GetRotationMatrix3D(itk::Array<double> rotationParameters)
 159   {
 160     itk::Matrix<double, 3, 3> matrix;
 161     //rotational part
 162     matrix[0][0]= cos(rotationParameters[1])*cos(rotationParameters[2]);
 163     matrix[0][1]= sin(rotationParameters[0])*sin(rotationParameters[1])*cos(rotationParameters[2])+ sin(rotationParameters[2])*cos(rotationParameters[0]);
 164     matrix[0][2]= -cos(rotationParameters[0])*sin(rotationParameters[1])*cos(rotationParameters[2])+sin(rotationParameters[0])*sin(rotationParameters[2]);
 165     matrix[1][0]= -cos(rotationParameters[1])*sin(rotationParameters[2]);
 166     matrix[1][1]= -sin(rotationParameters[0])*sin(rotationParameters[1])*sin(rotationParameters[2])+cos(rotationParameters[0])*cos(rotationParameters[2]);
 167     matrix[1][2]= cos(rotationParameters[0])*sin(rotationParameters[1])*sin(rotationParameters[2])+sin(rotationParameters[0])*cos(rotationParameters[2]);
 168     matrix[2][0]= sin(rotationParameters[1]);
 169     matrix[2][1]= -sin(rotationParameters[0])*cos(rotationParameters[1]);
 170     matrix[2][2]= cos(rotationParameters[0])*cos(rotationParameters[1]);
 171     return matrix;
 172   }
 173
 174   inline itk::Matrix<double, 2, 2> GetRotationMatrix2D(itk::Array<double> rotationParameters)
 175   {
 176     itk::Matrix<double, 2, 2> matrix;
 177     matrix[0][0] = cos(rotationParameters[0]);
 178     matrix[1][0] = sin(rotationParameters[0]);
 179     matrix[0][1] = -matrix[1][0];
 180     matrix[1][1] = matrix[0][0];
 181     return matrix;
 182   }
 183
 184   //========================================================================================
 185   inline  itk::Point<double, 3> GetRotatedPoint3D(itk::Array<double> rotationParameters, itk::Point<double, 3> input)
 186   {
 187     itk::Matrix<double, 3, 3> matrix = GetRotationMatrix3D(rotationParameters);
 188     itk::Point<double, 3> output;
 189     for (unsigned int i=0;i<3;i++)
 190       {
 191         output[i]=0.0;
 192         for (unsigned int j=0;j<3;j++)
 193           output[i]+=matrix(i,j)*input[j];
 194       }
 195     return output;
 196   }
 197
 198
 199   inline itk::Matrix<double, 4, 4> GetCenteredRotationMatrix3D(itk::Array<double> rotationParameters,itk::Point<double,3> centerOfRotation )
 200   {
 201     //rotational part is identical as affine matrix, translations change
 202     itk::Array<double> parameters(6);
 203     for(unsigned int i=0; i<3;i++) parameters[i]=rotationParameters[i];
 204     for(unsigned int i=3; i<6;i++) parameters[i]=centerOfRotation[i-3];
 205     itk::Matrix<double, 4, 4> matrix=GetForwardAffineMatrix3D(parameters);
 206
 207     //Get the rotation of the centerOfRotation
 208     itk::Matrix<double,3,3> rotation = GetRotationalPartMatrix3D(matrix);
 209     itk::Point<double,3> rotatedCenter=rotation*centerOfRotation; //GetRotatedPoint3D(rotationParameters, centerOfRotation);
 210
 211     //Substract this point to the translational part
 212     matrix(0,3)-=rotatedCenter[0];
 213     matrix(1,3)-=rotatedCenter[1];
 214     matrix(2,3)-=rotatedCenter[2];
 215     return matrix;
 216   }
 217
 218
 219   //   inline  itk::Matrix<double, 4, 4> GetComposedMatrix3D(itk::Matrix<double, 4, 4> firstAppliedTransform, itk::Matrix<double, 4, 4> secondAppliedTransform)
 220   //   {
 221   //     itk::Matrix<double, 4, 4> matrix;
 222   //     for (unsigned int i=0;i<4;i++)
 223   //       for (unsigned int j=0;j<4;j++)
 224   //    {
 225   //      matrix[i][j]=0.0;
 226   //      for (unsigned int k=0;k<4;k++)
 227   //        matrix[i][j]+=firstAppliedTransform[i][k]*secondAppliedTransform[k][j];
 228   //    }
 229   //     return matrix;
 230   //   }
 231
 232
 233   //========================================================================================
 234   inline itk::Matrix<double, 5, 5> ReadMatrix4D(std::string fileName)
 235   {
 236     // read input matrix
 237     std::ifstream is;
 238     openFileForReading(is, fileName);
 239     std::vector<double> nb;
 240     double x;
 241     skipComment(is);
 242     is >> x;
 243     while (!is.eof()) {
 244       nb.push_back(x);
 245       skipComment(is);
 246       is >> x;
 247     }
 248
 249     //copy it to the matrix
 250     itk::Matrix<double, 5, 5> matrix;
 251     unsigned int index=0;
 252     for (unsigned int i=0;i<5;i++)
 253       for (unsigned int j=0;j<5;j++)
 254         matrix[i][j]=nb[index++];
 255     return matrix;
 256   }
 257
 258   inline itk::Matrix<double, 4, 4> ReadMatrix3D(std::string fileName)
 259   {
 260     // read input matrix
 261     std::ifstream is;
 262     openFileForReading(is, fileName);
 263     std::vector<double> nb;
 264     double x;
 265     skipComment(is);
 266     is >> x;
 267     while (!is.eof()) {
 268       nb.push_back(x);
 269       skipComment(is);
 270       is >> x;
 271     }
 272
 273     //copy it to the matrix
 274     itk::Matrix<double, 4, 4> matrix;
 275     unsigned int index=0;
 276     for (unsigned int i=0;i<4;i++)
 277       for (unsigned int j=0;j<4;j++)
 278         matrix[i][j]=nb[index++];
 279     return matrix;
 280   }
 281
 282   inline vtkMatrix4x4* ReadVTKMatrix3D(std::string fileName) {
 283     // read input matrix
 284     std::ifstream is;
 285     openFileForReading(is, fileName);
 286     std::vector<double> nb;
 287     double x;
 288     skipComment(is);
 289     is >> x;
 290     while (!is.eof()) {
 291       nb.push_back(x);
 292       skipComment(is);
 293       is >> x;
 294     }
 295
 296     vtkSmartPointer<vtkMatrix4x4> matrix = vtkSmartPointer<vtkMatrix4x4>::New();
 297     unsigned int index=0;
 298     for (unsigned int i=0;i<4;i++)
 299       for (unsigned int j=0;j<4;j++)
 300          matrix->SetElement(i,j, nb[index++]);
 301
 302     return matrix;
 303   }
 304
 305   inline itk::Matrix<double, 3, 3> ReadMatrix2D(std::string fileName)
 306   {
 307     // read input matrix
 308     std::ifstream is;
 309     openFileForReading(is, fileName);
 310     std::vector<double> nb;
 311     double x;
 312     skipComment(is);
 313     is >> x;
 314     while (!is.eof()) {
 315       nb.push_back(x);
 316       skipComment(is);
 317       is >> x;
 318     }
 319
 320     //copy it to the matrix
 321     itk::Matrix<double, 3, 3> matrix;
 322     unsigned int index=0;
 323     for (unsigned int i=0;i<3;i++)
 324       for (unsigned int j=0;j<3;j++)
 325         matrix[i][j]=nb[index++];
 326     return matrix;
 327   }
 328
 329   template <unsigned int Dimension > inline itk::Matrix<double, Dimension+1 , Dimension+1> ReadMatrix(std::string fileName)
 330   {
 331
 332     // read input matrix
 333     std::ifstream is;
 334     openFileForReading(is, fileName);
 335     std::vector<double> nb;
 336     double x;
 337     skipComment(is);
 338     is >> x;
 339     while (!is.eof()) {
 340       nb.push_back(x);
 341       skipComment(is);
 342       is >> x;
 343     }
 344
 345     //copy it to the matrix
 346     itk::Matrix<double, Dimension+1, Dimension+1> matrix;
 347     unsigned int index=0;
 348     for (unsigned int i=0;i<Dimension+1;i++)
 349       for (unsigned int j=0;j<Dimension+1;j++)
 350         matrix[i][j]=nb[index++];
 351     return matrix;
 352   }
 353
 354
 355   //    template<> inline itk::Matrix<double, 3, 3>  ReadMatrix<2> (std::string fileName)
 356   //    {
 357   //      return  ReadMatrix2D(fileName);
 358   //    }
 359   //    template<> inline itk::Matrix<double, 4, 4>  ReadMatrix<3> (std::string fileName)
 360   //    {
 361   //      return  ReadMatrix3D(fileName);
 362   //    }
 363   //  template<> inline itk::Matrix<double, 5, 5>  ReadMatrix<4> (std::string fileName)
 364   //    {
 365   //      return  ReadMatrix4D(fileName);
 366   //    }
 367
 368
 369   //========================================================================================
 370   inline  itk::Matrix<double, 4, 4> GetRotationalPartMatrix4D(itk::Matrix<double, 5, 5> input)
 371   {
 372     itk::Matrix<double,4,4> matrix;
 373     matrix[0][0]= input[0][0];
 374     matrix[0][1]= input[0][1];
 375     matrix[0][2]= input[0][2];
 376     matrix[0][3]= input[0][3];
 377     matrix[1][0]= input[1][0];
 378     matrix[1][1]= input[1][1];
 379     matrix[1][2]= input[1][2];
 380     matrix[1][3]= input[1][3];
 381     matrix[2][0]= input[2][0];
 382     matrix[2][1]= input[2][1];
 383     matrix[2][2]= input[2][2];
 384     matrix[2][2]= input[2][2];
 385     matrix[2][3]= input[2][3];
 386     matrix[3][0]= input[3][0];
 387     matrix[3][1]= input[3][1];
 388     matrix[3][2]= input[3][2];
 389     matrix[3][2]= input[3][2];
 390     matrix[3][3]= input[3][3];
 391
 392     return matrix;
 393   }
 394
 395
 396   inline  itk::Matrix<double, 3, 3> GetRotationalPartMatrix3D(itk::Matrix<double, 4, 4> input)
 397   {
 398     itk::Matrix<double,3,3> matrix;
 399     matrix[0][0]= input[0][0];
 400     matrix[0][1]= input[0][1];
 401     matrix[0][2]= input[0][2];
 402     matrix[1][0]= input[1][0];
 403     matrix[1][1]= input[1][1];
 404     matrix[1][2]= input[1][2];
 405     matrix[2][0]= input[2][0];
 406     matrix[2][1]= input[2][1];
 407     matrix[2][2]= input[2][2];
 408     return matrix;
 409   }
 410
 411   inline  itk::Matrix<double, 2, 2> GetRotationalPartMatrix2D(itk::Matrix<double, 3, 3> input)
 412   {
 413     itk::Matrix<double,2,2> matrix;
 414     matrix[0][0]= input[0][0];
 415     matrix[0][1]= input[0][1];
 416     matrix[0][2]= input[0][2];
 417     matrix[1][0]= input[1][0];
 418     matrix[1][1]= input[1][1];
 419     return matrix;
 420   }
 421
 422   inline itk::Matrix<double, 4, 4> GetRotationalPartMatrix(itk::Matrix<double, 5, 5> input)
 423   {
 424     return GetRotationalPartMatrix4D(input);
 425   }
 426
 427   inline itk::Matrix<double, 3, 3> GetRotationalPartMatrix(itk::Matrix<double, 4, 4> input)
 428   {
 429     return GetRotationalPartMatrix3D(input);
 430   }
 431
 432   inline itk::Matrix<double, 2, 2> GetRotationalPartMatrix(itk::Matrix<double, 3, 3> input)
 433   {
 434     return GetRotationalPartMatrix2D(input);
 435   }
 436
 437
 438   //========================================================================================
 439   inline  itk::Vector<double,4> GetTranslationPartMatrix4D(itk::Matrix<double, 5, 5> input)
 440   {
 441     itk::Vector<double,4> vec;
 442     vec[0]= input[0][4];
 443     vec[1]= input[1][4];
 444     vec[2]= input[2][4];
 445     vec[3]= input[3][4];
 446     return vec;
 447
 448   }
 449
 450   inline  itk::Vector<double,3> GetTranslationPartMatrix3D(itk::Matrix<double, 4, 4> input)
 451   {
 452     itk::Vector<double,3> vec;
 453     vec[0]= input[0][3];
 454     vec[1]= input[1][3];
 455     vec[2]= input[2][3];
 456     return vec;
 457
 458   }
 459   inline  itk::Vector<double,2> GetTranslationPartMatrix2D(itk::Matrix<double, 3, 3> input)
 460   {
 461     itk::Vector<double,2> vec;
 462     vec[0]= input[0][2];
 463     vec[1]= input[1][2];
 464     return vec;
 465
 466   }
 467
 468   inline itk::Vector<double,4> GetTranslationPartMatrix(itk::Matrix<double, 5, 5> input)
 469   {
 470
 471     return GetTranslationPartMatrix4D(input);
 472   }
 473
 474   inline itk::Vector<double,3> GetTranslationPartMatrix(itk::Matrix<double, 4, 4> input)
 475   {
 476
 477     return GetTranslationPartMatrix3D(input);
 478   }
 479
 480   inline itk::Vector<double,2> GetTranslationPartMatrix(itk::Matrix<double, 3, 3> input)
 481   {
 482
 483     return GetTranslationPartMatrix2D(input);
 484   }
 485 }
 486
 487 #endif //#define CLITKTRANSFORMUTILITIES_H