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