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