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   itk::Matrix<double, 4, 4> ReadMatrix3D(std::string fileName);
 259
 260   inline vtkMatrix4x4* ReadVTKMatrix3D(std::string fileName) {
 261     // read input matrix
 262     std::ifstream is;
 263     openFileForReading(is, fileName);
 264     std::vector<double> nb;
 265     double x;
 266     skipComment(is);
 267     is >> x;
 268     while (!is.eof()) {
 269       nb.push_back(x);
 270       skipComment(is);
 271       is >> x;
 272     }
 273
 274     vtkSmartPointer<vtkMatrix4x4> matrix = vtkSmartPointer<vtkMatrix4x4>::New();
 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->SetElement(i,j, nb[index++]);
 279
 280     return matrix;
 281   }
 282
 283   inline itk::Matrix<double, 3, 3> ReadMatrix2D(std::string fileName)
 284   {
 285     // read input matrix
 286     std::ifstream is;
 287     openFileForReading(is, fileName);
 288     std::vector<double> nb;
 289     double x;
 290     skipComment(is);
 291     is >> x;
 292     while (!is.eof()) {
 293       nb.push_back(x);
 294       skipComment(is);
 295       is >> x;
 296     }
 297
 298     //copy it to the matrix
 299     itk::Matrix<double, 3, 3> matrix;
 300     unsigned int index=0;
 301     for (unsigned int i=0;i<3;i++)
 302       for (unsigned int j=0;j<3;j++)
 303         matrix[i][j]=nb[index++];
 304     return matrix;
 305   }
 306
 307   template <unsigned int Dimension > inline itk::Matrix<double, Dimension+1 , Dimension+1> ReadMatrix(std::string fileName)
 308   {
 309
 310     // read input matrix
 311     std::ifstream is;
 312     openFileForReading(is, fileName);
 313     std::vector<double> nb;
 314     double x;
 315     skipComment(is);
 316     is >> x;
 317     while (!is.eof()) {
 318       nb.push_back(x);
 319       skipComment(is);
 320       is >> x;
 321     }
 322
 323     //copy it to the matrix
 324     itk::Matrix<double, Dimension+1, Dimension+1> matrix;
 325     unsigned int index=0;
 326     for (unsigned int i=0;i<Dimension+1;i++)
 327       for (unsigned int j=0;j<Dimension+1;j++)
 328         matrix[i][j]=nb[index++];
 329     return matrix;
 330   }
 331
 332
 333   //    template<> inline itk::Matrix<double, 3, 3>  ReadMatrix<2> (std::string fileName)
 334   //    {
 335   //      return  ReadMatrix2D(fileName);
 336   //    }
 337   //    template<> inline itk::Matrix<double, 4, 4>  ReadMatrix<3> (std::string fileName)
 338   //    {
 339   //      return  ReadMatrix3D(fileName);
 340   //    }
 341   //  template<> inline itk::Matrix<double, 5, 5>  ReadMatrix<4> (std::string fileName)
 342   //    {
 343   //      return  ReadMatrix4D(fileName);
 344   //    }
 345
 346
 347   //========================================================================================
 348   inline  itk::Matrix<double, 4, 4> GetRotationalPartMatrix4D(itk::Matrix<double, 5, 5> input)
 349   {
 350     itk::Matrix<double,4,4> matrix;
 351     matrix[0][0]= input[0][0];
 352     matrix[0][1]= input[0][1];
 353     matrix[0][2]= input[0][2];
 354     matrix[0][3]= input[0][3];
 355     matrix[1][0]= input[1][0];
 356     matrix[1][1]= input[1][1];
 357     matrix[1][2]= input[1][2];
 358     matrix[1][3]= input[1][3];
 359     matrix[2][0]= input[2][0];
 360     matrix[2][1]= input[2][1];
 361     matrix[2][2]= input[2][2];
 362     matrix[2][2]= input[2][2];
 363     matrix[2][3]= input[2][3];
 364     matrix[3][0]= input[3][0];
 365     matrix[3][1]= input[3][1];
 366     matrix[3][2]= input[3][2];
 367     matrix[3][2]= input[3][2];
 368     matrix[3][3]= input[3][3];
 369
 370     return matrix;
 371   }
 372
 373
 374   inline  itk::Matrix<double, 3, 3> GetRotationalPartMatrix3D(itk::Matrix<double, 4, 4> input)
 375   {
 376     itk::Matrix<double,3,3> matrix;
 377     matrix[0][0]= input[0][0];
 378     matrix[0][1]= input[0][1];
 379     matrix[0][2]= input[0][2];
 380     matrix[1][0]= input[1][0];
 381     matrix[1][1]= input[1][1];
 382     matrix[1][2]= input[1][2];
 383     matrix[2][0]= input[2][0];
 384     matrix[2][1]= input[2][1];
 385     matrix[2][2]= input[2][2];
 386     return matrix;
 387   }
 388
 389   inline  itk::Matrix<double, 2, 2> GetRotationalPartMatrix2D(itk::Matrix<double, 3, 3> input)
 390   {
 391     itk::Matrix<double,2,2> matrix;
 392     matrix[0][0]= input[0][0];
 393     matrix[0][1]= input[0][1];
 394     matrix[0][2]= input[0][2];
 395     matrix[1][0]= input[1][0];
 396     matrix[1][1]= input[1][1];
 397     return matrix;
 398   }
 399
 400   inline itk::Matrix<double, 4, 4> GetRotationalPartMatrix(itk::Matrix<double, 5, 5> input)
 401   {
 402     return GetRotationalPartMatrix4D(input);
 403   }
 404
 405   inline itk::Matrix<double, 3, 3> GetRotationalPartMatrix(itk::Matrix<double, 4, 4> input)
 406   {
 407     return GetRotationalPartMatrix3D(input);
 408   }
 409
 410   inline itk::Matrix<double, 2, 2> GetRotationalPartMatrix(itk::Matrix<double, 3, 3> input)
 411   {
 412     return GetRotationalPartMatrix2D(input);
 413   }
 414
 415
 416   //========================================================================================
 417   inline  itk::Vector<double,4> GetTranslationPartMatrix4D(itk::Matrix<double, 5, 5> input)
 418   {
 419     itk::Vector<double,4> vec;
 420     vec[0]= input[0][4];
 421     vec[1]= input[1][4];
 422     vec[2]= input[2][4];
 423     vec[3]= input[3][4];
 424     return vec;
 425
 426   }
 427
 428   inline  itk::Vector<double,3> GetTranslationPartMatrix3D(itk::Matrix<double, 4, 4> input)
 429   {
 430     itk::Vector<double,3> vec;
 431     vec[0]= input[0][3];
 432     vec[1]= input[1][3];
 433     vec[2]= input[2][3];
 434     return vec;
 435
 436   }
 437   inline  itk::Vector<double,2> GetTranslationPartMatrix2D(itk::Matrix<double, 3, 3> input)
 438   {
 439     itk::Vector<double,2> vec;
 440     vec[0]= input[0][2];
 441     vec[1]= input[1][2];
 442     return vec;
 443
 444   }
 445
 446   inline itk::Vector<double,4> GetTranslationPartMatrix(itk::Matrix<double, 5, 5> input)
 447   {
 448
 449     return GetTranslationPartMatrix4D(input);
 450   }
 451
 452   inline itk::Vector<double,3> GetTranslationPartMatrix(itk::Matrix<double, 4, 4> input)
 453   {
 454
 455     return GetTranslationPartMatrix3D(input);
 456   }
 457
 458   inline itk::Vector<double,2> GetTranslationPartMatrix(itk::Matrix<double, 3, 3> input)
 459   {
 460
 461     return GetTranslationPartMatrix2D(input);
 462   }
 463 }
 464
 465 #endif //#define CLITKTRANSFORMUTILITIES_H