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
  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   template <unsigned int Dimension >
 155   inline itk::Matrix<double, Dimension+1, Dimension+1>
 156   GetBackwardAffineMatrix(itk::Array<double> transformParameters)
 157   {
 158   }
 159
 160   inline itk::Matrix<double, 3, 3> GetRotationMatrix3D(itk::Array<double> rotationParameters)
 161   {
 162     itk::Matrix<double, 3, 3> matrix;
 163     //rotational part
 164     matrix[0][0]= cos(rotationParameters[1])*cos(rotationParameters[2]);
 165     matrix[0][1]= sin(rotationParameters[0])*sin(rotationParameters[1])*cos(rotationParameters[2])+ sin(rotationParameters[2])*cos(rotationParameters[0]);
 166     matrix[0][2]= -cos(rotationParameters[0])*sin(rotationParameters[1])*cos(rotationParameters[2])+sin(rotationParameters[0])*sin(rotationParameters[2]);
 167     matrix[1][0]= -cos(rotationParameters[1])*sin(rotationParameters[2]);
 168     matrix[1][1]= -sin(rotationParameters[0])*sin(rotationParameters[1])*sin(rotationParameters[2])+cos(rotationParameters[0])*cos(rotationParameters[2]);
 169     matrix[1][2]= cos(rotationParameters[0])*sin(rotationParameters[1])*sin(rotationParameters[2])+sin(rotationParameters[0])*cos(rotationParameters[2]);
 170     matrix[2][0]= sin(rotationParameters[1]);
 171     matrix[2][1]= -sin(rotationParameters[0])*cos(rotationParameters[1]);
 172     matrix[2][2]= cos(rotationParameters[0])*cos(rotationParameters[1]);
 173     return matrix;
 174   }
 175
 176
 177
 178   //========================================================================================
 179   inline  itk::Point<double, 3> GetRotatedPoint3D(itk::Array<double> rotationParameters, itk::Point<double, 3> input)
 180   {
 181     itk::Matrix<double, 3, 3> matrix = GetRotationMatrix3D(rotationParameters);
 182     itk::Point<double, 3> output;
 183     for (unsigned int i=0;i<3;i++)
 184       {
 185         output[i]=0.0;
 186         for (unsigned int j=0;j<3;j++)
 187           output[i]+=matrix(i,j)*input[j];
 188       }
 189     return output;
 190   }
 191
 192
 193   inline itk::Matrix<double, 4, 4> GetCenteredRotationMatrix3D(itk::Array<double> rotationParameters,itk::Point<double,3> centerOfRotation )
 194   {
 195     //rotational part is identical as affine matrix, translations change
 196     itk::Array<double> parameters(6);
 197     for(unsigned int i=0; i<3;i++) parameters[i]=rotationParameters[i];
 198     for(unsigned int i=3; i<6;i++) parameters[i]=centerOfRotation[i-3];
 199     itk::Matrix<double, 4, 4> matrix=GetForwardAffineMatrix3D(parameters);
 200
 201     //Get the rotation of the centerOfRotation
 202     itk::Matrix<double,3,3> rotation = GetRotationalPartMatrix3D(matrix);
 203     itk::Point<double,3> rotatedCenter=rotation*centerOfRotation; //GetRotatedPoint3D(rotationParameters, centerOfRotation);
 204
 205     //Substract this point to the translational part
 206     matrix(0,3)-=rotatedCenter[0];
 207     matrix(1,3)-=rotatedCenter[1];
 208     matrix(2,3)-=rotatedCenter[2];
 209     return matrix;
 210   }
 211
 212
 213   //   inline  itk::Matrix<double, 4, 4> GetComposedMatrix3D(itk::Matrix<double, 4, 4> firstAppliedTransform, itk::Matrix<double, 4, 4> secondAppliedTransform)
 214   //   {
 215   //     itk::Matrix<double, 4, 4> matrix;
 216   //     for (unsigned int i=0;i<4;i++)
 217   //       for (unsigned int j=0;j<4;j++)
 218   //    {
 219   //      matrix[i][j]=0.0;
 220   //      for (unsigned int k=0;k<4;k++)
 221   //        matrix[i][j]+=firstAppliedTransform[i][k]*secondAppliedTransform[k][j];
 222   //    }
 223   //     return matrix;
 224   //   }
 225
 226
 227   //========================================================================================
 228   inline itk::Matrix<double, 5, 5> ReadMatrix4D(std::string fileName)
 229   {
 230     // read input matrix
 231     std::ifstream is;
 232     openFileForReading(is, fileName);
 233     std::vector<double> nb;
 234     double x;
 235     skipComment(is);
 236     is >> x;
 237     while (!is.eof()) {
 238       nb.push_back(x);
 239       skipComment(is);
 240       is >> x;
 241     }
 242
 243     //copy it to the matrix
 244     itk::Matrix<double, 5, 5> matrix;
 245     unsigned int index=0;
 246     for (unsigned int i=0;i<5;i++)
 247       for (unsigned int j=0;j<5;j++)
 248         matrix[i][j]=nb[index++];
 249     return matrix;
 250   }
 251
 252   inline itk::Matrix<double, 4, 4> ReadMatrix3D(std::string fileName)
 253   {
 254     // read input matrix
 255     std::ifstream is;
 256     openFileForReading(is, fileName);
 257     std::vector<double> nb;
 258     double x;
 259     skipComment(is);
 260     is >> x;
 261     while (!is.eof()) {
 262       nb.push_back(x);
 263       skipComment(is);
 264       is >> x;
 265     }
 266
 267     //copy it to the matrix
 268     itk::Matrix<double, 4, 4> matrix;
 269     unsigned int index=0;
 270     for (unsigned int i=0;i<4;i++)
 271       for (unsigned int j=0;j<4;j++)
 272         matrix[i][j]=nb[index++];
 273     return matrix;
 274   }
 275
 276   inline itk::Matrix<double, 3, 3> ReadMatrix2D(std::string fileName)
 277   {
 278     // read input matrix
 279     std::ifstream is;
 280     openFileForReading(is, fileName);
 281     std::vector<double> nb;
 282     double x;
 283     skipComment(is);
 284     is >> x;
 285     while (!is.eof()) {
 286       nb.push_back(x);
 287       skipComment(is);
 288       is >> x;
 289     }
 290
 291     //copy it to the matrix
 292     itk::Matrix<double, 3, 3> matrix;
 293     unsigned int index=0;
 294     for (unsigned int i=0;i<3;i++)
 295       for (unsigned int j=0;j<3;j++)
 296         matrix[i][j]=nb[index++];
 297     return matrix;
 298   }
 299
 300   template <unsigned int Dimension > inline itk::Matrix<double, Dimension+1 , Dimension+1> ReadMatrix(std::string fileName)
 301   {
 302
 303     // read input matrix
 304     std::ifstream is;
 305     openFileForReading(is, fileName);
 306     std::vector<double> nb;
 307     double x;
 308     skipComment(is);
 309     is >> x;
 310     while (!is.eof()) {
 311       nb.push_back(x);
 312       skipComment(is);
 313       is >> x;
 314     }
 315
 316     //copy it to the matrix
 317     itk::Matrix<double, Dimension+1, Dimension+1> matrix;
 318     unsigned int index=0;
 319     for (unsigned int i=0;i<Dimension+1;i++)
 320       for (unsigned int j=0;j<Dimension+1;j++)
 321         matrix[i][j]=nb[index++];
 322     return matrix;
 323   }
 324
 325
 326   //    template<> inline itk::Matrix<double, 3, 3>  ReadMatrix<2> (std::string fileName)
 327   //    {
 328   //      return  ReadMatrix2D(fileName);
 329   //    }
 330   //    template<> inline itk::Matrix<double, 4, 4>  ReadMatrix<3> (std::string fileName)
 331   //    {
 332   //      return  ReadMatrix3D(fileName);
 333   //    }
 334   //  template<> inline itk::Matrix<double, 5, 5>  ReadMatrix<4> (std::string fileName)
 335   //    {
 336   //      return  ReadMatrix4D(fileName);
 337   //    }
 338
 339
 340   //========================================================================================
 341   inline  itk::Matrix<double, 4, 4> GetRotationalPartMatrix4D(itk::Matrix<double, 5, 5> input)
 342   {
 343     itk::Matrix<double,4,4> matrix;
 344     matrix[0][0]= input[0][0];
 345     matrix[0][1]= input[0][1];
 346     matrix[0][2]= input[0][2];
 347     matrix[0][3]= input[0][3];
 348     matrix[1][0]= input[1][0];
 349     matrix[1][1]= input[1][1];
 350     matrix[1][2]= input[1][2];
 351     matrix[1][3]= input[1][3];
 352     matrix[2][0]= input[2][0];
 353     matrix[2][1]= input[2][1];
 354     matrix[2][2]= input[2][2];
 355     matrix[2][2]= input[2][2];
 356     matrix[2][3]= input[2][3];
 357     matrix[3][0]= input[3][0];
 358     matrix[3][1]= input[3][1];
 359     matrix[3][2]= input[3][2];
 360     matrix[3][2]= input[3][2];
 361     matrix[3][3]= input[3][3];
 362
 363     return matrix;
 364   }
 365
 366
 367   inline  itk::Matrix<double, 3, 3> GetRotationalPartMatrix3D(itk::Matrix<double, 4, 4> input)
 368   {
 369     itk::Matrix<double,3,3> matrix;
 370     matrix[0][0]= input[0][0];
 371     matrix[0][1]= input[0][1];
 372     matrix[0][2]= input[0][2];
 373     matrix[1][0]= input[1][0];
 374     matrix[1][1]= input[1][1];
 375     matrix[1][2]= input[1][2];
 376     matrix[2][0]= input[2][0];
 377     matrix[2][1]= input[2][1];
 378     matrix[2][2]= input[2][2];
 379     return matrix;
 380   }
 381
 382   inline  itk::Matrix<double, 2, 2> GetRotationalPartMatrix2D(itk::Matrix<double, 3, 3> input)
 383   {
 384     itk::Matrix<double,2,2> matrix;
 385     matrix[0][0]= input[0][0];
 386     matrix[0][1]= input[0][1];
 387     matrix[0][2]= input[0][2];
 388     matrix[1][0]= input[1][0];
 389     matrix[1][1]= input[1][1];
 390     return matrix;
 391   }
 392
 393   inline itk::Matrix<double, 4, 4> GetRotationalPartMatrix(itk::Matrix<double, 5, 5> input)
 394   {
 395     return GetRotationalPartMatrix4D(input);
 396   }
 397
 398   inline itk::Matrix<double, 3, 3> GetRotationalPartMatrix(itk::Matrix<double, 4, 4> input)
 399   {
 400     return GetRotationalPartMatrix3D(input);
 401   }
 402
 403   inline itk::Matrix<double, 2, 2> GetRotationalPartMatrix(itk::Matrix<double, 3, 3> input)
 404   {
 405     return GetRotationalPartMatrix2D(input);
 406   }
 407
 408
 409   //========================================================================================
 410   inline  itk::Vector<double,4> GetTranslationPartMatrix4D(itk::Matrix<double, 5, 5> input)
 411   {
 412     itk::Vector<double,4> vec;
 413     vec[0]= input[0][4];
 414     vec[1]= input[1][4];
 415     vec[2]= input[2][4];
 416     vec[3]= input[3][4];
 417     return vec;
 418
 419   }
 420
 421   inline  itk::Vector<double,3> GetTranslationPartMatrix3D(itk::Matrix<double, 4, 4> input)
 422   {
 423     itk::Vector<double,3> vec;
 424     vec[0]= input[0][3];
 425     vec[1]= input[1][3];
 426     vec[2]= input[2][3];
 427     return vec;
 428
 429   }
 430   inline  itk::Vector<double,2> GetTranslationPartMatrix2D(itk::Matrix<double, 3, 3> input)
 431   {
 432     itk::Vector<double,2> vec;
 433     vec[0]= input[0][2];
 434     vec[1]= input[1][2];
 435     return vec;
 436
 437   }
 438
 439   inline itk::Vector<double,4> GetTranslationPartMatrix(itk::Matrix<double, 5, 5> input)
 440   {
 441
 442     return GetTranslationPartMatrix4D(input);
 443   }
 444
 445   inline itk::Vector<double,3> GetTranslationPartMatrix(itk::Matrix<double, 4, 4> input)
 446   {
 447
 448     return GetTranslationPartMatrix3D(input);
 449   }
 450
 451   inline itk::Vector<double,2> GetTranslationPartMatrix(itk::Matrix<double, 3, 3> input)
 452   {
 453
 454     return GetTranslationPartMatrix2D(input);
 455   }
 456 }
 457
 458 #endif //#define CLITKTRANSFORMUTILITIES_H