1 /*=========================================================================
2 Program: vv http://www.creatis.insa-lyon.fr/rio/vv
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
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.
13 It is distributed under dual licence
15 - BSD See included LICENSE.txt file
16 - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
17 ===========================================================================**/
19 #ifndef CLITKTRANSFORMUTILITIES_H
20 #define CLITKTRANSFORMUTILITIES_H
21 #include "itkMatrix.h"
24 #include "clitkImageCommon.h"
25 #include "clitkCommon.h"
26 #include <vtkMatrix4x4.h>
27 #include <vtkSmartPointer.h>
31 //============================================================================
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);
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);
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);
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);
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);
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);
65 //============================================================================
66 //Inline functions definition in header file, otherwise linker errors
67 //============================================================================
69 //========================================================================================
70 inline itk::Matrix<double, 3, 3> GetForwardAffineMatrix2D(itk::Array<double> transformParameters)
72 itk::Matrix<double, 3, 3> matrix;
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]);
79 matrix[0][2]=transformParameters[1];
80 matrix[1][2]=transformParameters[2];
88 inline itk::Matrix<double, 4, 4> GetForwardAffineMatrix3D(itk::Array<double> transformParameters)
90 itk::Matrix<double, 4, 4> matrix;
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]);
102 matrix[0][3]=transformParameters[3];
103 matrix[1][3]=transformParameters[4];
104 matrix[2][3]=transformParameters[5];
114 inline itk::Matrix<double, 3, 3> GetBackwardAffineMatrix2D(itk::Array<double> transformParameters)
116 itk::Matrix<double, 3, 3> matrix;
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]);
123 matrix[0][2]=transformParameters[1];
124 matrix[1][2]=transformParameters[2];
133 inline itk::Matrix<double, 4, 4> GetBackwardAffineMatrix3D(itk::Array<double> transformParameters)
135 itk::Matrix<double, 4, 4> matrix;
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]);
147 matrix[0][3]=transformParameters[3];
148 matrix[1][3]=transformParameters[4];
149 matrix[2][3]=transformParameters[5];
158 inline itk::Matrix<double, 3, 3> GetRotationMatrix3D(itk::Array<double> rotationParameters)
160 itk::Matrix<double, 3, 3> matrix;
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]);
174 inline itk::Matrix<double, 2, 2> GetRotationMatrix2D(itk::Array<double> rotationParameters)
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];
184 //========================================================================================
185 inline itk::Point<double, 3> GetRotatedPoint3D(itk::Array<double> rotationParameters, itk::Point<double, 3> input)
187 itk::Matrix<double, 3, 3> matrix = GetRotationMatrix3D(rotationParameters);
188 itk::Point<double, 3> output;
189 for (unsigned int i=0;i<3;i++)
192 for (unsigned int j=0;j<3;j++)
193 output[i]+=matrix(i,j)*input[j];
199 inline itk::Matrix<double, 4, 4> GetCenteredRotationMatrix3D(itk::Array<double> rotationParameters,itk::Point<double,3> centerOfRotation )
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);
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);
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];
219 // inline itk::Matrix<double, 4, 4> GetComposedMatrix3D(itk::Matrix<double, 4, 4> firstAppliedTransform, itk::Matrix<double, 4, 4> secondAppliedTransform)
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++)
226 // for (unsigned int k=0;k<4;k++)
227 // matrix[i][j]+=firstAppliedTransform[i][k]*secondAppliedTransform[k][j];
233 //========================================================================================
234 inline itk::Matrix<double, 5, 5> ReadMatrix4D(std::string fileName)
238 openFileForReading(is, fileName);
239 std::vector<double> nb;
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++];
258 inline itk::Matrix<double, 4, 4> ReadMatrix3D(std::string fileName)
262 openFileForReading(is, fileName);
263 std::vector<double> nb;
273 //copy it to the matrix
274 itk::Matrix<double, 4, 4> matrix;
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[i][j]=nb[index++];
282 inline vtkMatrix4x4* ReadVTKMatrix3D(std::string fileName) {
285 openFileForReading(is, fileName);
286 std::vector<double> nb;
296 vtkSmartPointer<vtkMatrix4x4> matrix = vtkSmartPointer<vtkMatrix4x4>::New();
297 unsigned int index=0;
298 for (unsigned int i=0;i<4;i++)
299 for (unsigned int j=0;j<4;j++)
300 matrix->SetElement(i,j, nb[index++]);
305 inline itk::Matrix<double, 3, 3> ReadMatrix2D(std::string fileName)
309 openFileForReading(is, fileName);
310 std::vector<double> nb;
320 //copy it to the matrix
321 itk::Matrix<double, 3, 3> matrix;
322 unsigned int index=0;
323 for (unsigned int i=0;i<3;i++)
324 for (unsigned int j=0;j<3;j++)
325 matrix[i][j]=nb[index++];
329 template <unsigned int Dimension > inline itk::Matrix<double, Dimension+1 , Dimension+1> ReadMatrix(std::string fileName)
334 openFileForReading(is, fileName);
335 std::vector<double> nb;
345 //copy it to the matrix
346 itk::Matrix<double, Dimension+1, Dimension+1> matrix;
347 unsigned int index=0;
348 for (unsigned int i=0;i<Dimension+1;i++)
349 for (unsigned int j=0;j<Dimension+1;j++)
350 matrix[i][j]=nb[index++];
355 // template<> inline itk::Matrix<double, 3, 3> ReadMatrix<2> (std::string fileName)
357 // return ReadMatrix2D(fileName);
359 // template<> inline itk::Matrix<double, 4, 4> ReadMatrix<3> (std::string fileName)
361 // return ReadMatrix3D(fileName);
363 // template<> inline itk::Matrix<double, 5, 5> ReadMatrix<4> (std::string fileName)
365 // return ReadMatrix4D(fileName);
369 //========================================================================================
370 inline itk::Matrix<double, 4, 4> GetRotationalPartMatrix4D(itk::Matrix<double, 5, 5> input)
372 itk::Matrix<double,4,4> matrix;
373 matrix[0][0]= input[0][0];
374 matrix[0][1]= input[0][1];
375 matrix[0][2]= input[0][2];
376 matrix[0][3]= input[0][3];
377 matrix[1][0]= input[1][0];
378 matrix[1][1]= input[1][1];
379 matrix[1][2]= input[1][2];
380 matrix[1][3]= input[1][3];
381 matrix[2][0]= input[2][0];
382 matrix[2][1]= input[2][1];
383 matrix[2][2]= input[2][2];
384 matrix[2][2]= input[2][2];
385 matrix[2][3]= input[2][3];
386 matrix[3][0]= input[3][0];
387 matrix[3][1]= input[3][1];
388 matrix[3][2]= input[3][2];
389 matrix[3][2]= input[3][2];
390 matrix[3][3]= input[3][3];
396 inline itk::Matrix<double, 3, 3> GetRotationalPartMatrix3D(itk::Matrix<double, 4, 4> input)
398 itk::Matrix<double,3,3> matrix;
399 matrix[0][0]= input[0][0];
400 matrix[0][1]= input[0][1];
401 matrix[0][2]= input[0][2];
402 matrix[1][0]= input[1][0];
403 matrix[1][1]= input[1][1];
404 matrix[1][2]= input[1][2];
405 matrix[2][0]= input[2][0];
406 matrix[2][1]= input[2][1];
407 matrix[2][2]= input[2][2];
411 inline itk::Matrix<double, 2, 2> GetRotationalPartMatrix2D(itk::Matrix<double, 3, 3> input)
413 itk::Matrix<double,2,2> matrix;
414 matrix[0][0]= input[0][0];
415 matrix[0][1]= input[0][1];
416 matrix[0][2]= input[0][2];
417 matrix[1][0]= input[1][0];
418 matrix[1][1]= input[1][1];
422 inline itk::Matrix<double, 4, 4> GetRotationalPartMatrix(itk::Matrix<double, 5, 5> input)
424 return GetRotationalPartMatrix4D(input);
427 inline itk::Matrix<double, 3, 3> GetRotationalPartMatrix(itk::Matrix<double, 4, 4> input)
429 return GetRotationalPartMatrix3D(input);
432 inline itk::Matrix<double, 2, 2> GetRotationalPartMatrix(itk::Matrix<double, 3, 3> input)
434 return GetRotationalPartMatrix2D(input);
438 //========================================================================================
439 inline itk::Vector<double,4> GetTranslationPartMatrix4D(itk::Matrix<double, 5, 5> input)
441 itk::Vector<double,4> vec;
450 inline itk::Vector<double,3> GetTranslationPartMatrix3D(itk::Matrix<double, 4, 4> input)
452 itk::Vector<double,3> vec;
459 inline itk::Vector<double,2> GetTranslationPartMatrix2D(itk::Matrix<double, 3, 3> input)
461 itk::Vector<double,2> vec;
468 inline itk::Vector<double,4> GetTranslationPartMatrix(itk::Matrix<double, 5, 5> input)
471 return GetTranslationPartMatrix4D(input);
474 inline itk::Vector<double,3> GetTranslationPartMatrix(itk::Matrix<double, 4, 4> input)
477 return GetTranslationPartMatrix3D(input);
480 inline itk::Vector<double,2> GetTranslationPartMatrix(itk::Matrix<double, 3, 3> input)
483 return GetTranslationPartMatrix2D(input);
487 #endif //#define CLITKTRANSFORMUTILITIES_H