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