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