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