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"
30 //============================================================================
32 //============================================================================
33 itk::Matrix<double, 3, 3> GetForwardAffineMatrix2D(itk::Array<double> transformParameters);
34 itk::Matrix<double, 4, 4> GetForwardAffineMatrix3D(itk::Array<double> transformParameters);
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);
40 itk::Matrix<double, 3, 3> GetRotationMatrix3D(itk::Array<double> rotationParameters);
41 itk::Matrix<double, 2, 2> GetRotationMatrix2D(itk::Array<double> rotationParameters);
42 template <unsigned int Dimension> itk::Matrix<double, Dimension, Dimension> GetRotationMatrix(itk::Array<double> rotationParameters);
43 itk::Point<double, 3> GetRotatedPoint3D(itk::Array<double> rotationParameters, itk::Point<double, 3> input);
44 itk::Matrix<double, 4, 4> GetCenteredRotationMatrix3D(itk::Array<double> rotationParameters,itk::Point<double,3> centerOfRotation);
45 // itk::Matrix<double, 4, 4> GetComposedMatrix3D(itk::Matrix<double, 4, 4> firstTransform, itk::Matrix<double, 4, 4> secondTransform);
47 itk::Matrix<double, 5, 5> ReadMatrix4D(std::string fileName);
48 itk::Matrix<double, 4, 4> ReadMatrix3D(std::string fileName);
49 itk::Matrix<double, 3, 3> ReadMatrix2D(std::string fileName);
50 template <unsigned int Dimension > itk::Matrix<double, Dimension+1 , Dimension+1> ReadMatrix(std::string fileName);
52 itk::Matrix<double, 3, 3> GetRotationalPartMatrix3D(itk::Matrix<double, 4, 4> input);
53 itk::Matrix<double, 3, 3> GetRotationalPartMatrix(itk::Matrix<double, 4, 4> input);
54 itk::Matrix<double, 2, 2> GetRotationalPartMatrix2D(itk::Matrix<double, 3, 3> input);
55 itk::Matrix<double, 2, 2> GetRotationalPartMatrix(itk::Matrix<double, 3, 3> input);
57 itk::Vector<double,3> GetTranslationPartMatrix3D(itk::Matrix<double, 4, 4> input);
58 itk::Vector<double,3> GetTranslationPartMatrix(itk::Matrix<double, 4, 4> input);
59 itk::Vector<double,2> GetTranslationPartMatrix2D(itk::Matrix<double, 3, 3> input);
60 itk::Vector<double,2> GetTranslationPartMatrix(itk::Matrix<double, 3, 3> input);
63 //============================================================================
64 //Inline functions definition in header file, otherwise linker errors
65 //============================================================================
67 //========================================================================================
68 inline itk::Matrix<double, 3, 3> GetForwardAffineMatrix2D(itk::Array<double> transformParameters)
70 itk::Matrix<double, 3, 3> matrix;
72 matrix[0][0]=cos(transformParameters[0]);
73 matrix[0][1]=-sin(transformParameters[0]);
74 matrix[1][0]=sin(transformParameters[0]);
75 matrix[1][1]=cos(transformParameters[0]);
77 matrix[0][2]=transformParameters[1];
78 matrix[1][2]=transformParameters[2];
86 inline itk::Matrix<double, 4, 4> GetForwardAffineMatrix3D(itk::Array<double> transformParameters)
88 itk::Matrix<double, 4, 4> matrix;
90 matrix[0][0]= cos(transformParameters[1])*cos(transformParameters[2]);
91 matrix[0][1]= sin(transformParameters[0])*sin(transformParameters[1])*cos(transformParameters[2])+ sin(transformParameters[2])*cos(transformParameters[0]);
92 matrix[0][2]= -cos(transformParameters[0])*sin(transformParameters[1])*cos(transformParameters[2])+sin(transformParameters[0])*sin(transformParameters[2]);
93 matrix[1][0]= -cos(transformParameters[1])*sin(transformParameters[2]);
94 matrix[1][1]= -sin(transformParameters[0])*sin(transformParameters[1])*sin(transformParameters[2])+cos(transformParameters[0])*cos(transformParameters[2]);
95 matrix[1][2]= cos(transformParameters[0])*sin(transformParameters[1])*sin(transformParameters[2])+sin(transformParameters[0])*cos(transformParameters[2]);
96 matrix[2][0]= sin(transformParameters[1]);
97 matrix[2][1]= -sin(transformParameters[0])*cos(transformParameters[1]);
98 matrix[2][2]= cos(transformParameters[0])*cos(transformParameters[1]);
100 matrix[0][3]=transformParameters[3];
101 matrix[1][3]=transformParameters[4];
102 matrix[2][3]=transformParameters[5];
112 inline itk::Matrix<double, 3, 3> GetBackwardAffineMatrix2D(itk::Array<double> transformParameters)
114 itk::Matrix<double, 3, 3> matrix;
116 matrix[0][0]=cos(transformParameters[0]);
117 matrix[0][1]=sin(transformParameters[0]);
118 matrix[1][0]=-sin(transformParameters[0]);
119 matrix[1][1]=cos(transformParameters[0]);
121 matrix[0][2]=transformParameters[1];
122 matrix[1][2]=transformParameters[2];
131 inline itk::Matrix<double, 4, 4> GetBackwardAffineMatrix3D(itk::Array<double> transformParameters)
133 itk::Matrix<double, 4, 4> matrix;
135 matrix[0][0]= cos(transformParameters[1])*cos(transformParameters[2]);
136 matrix[0][1]= sin(transformParameters[0])*sin(transformParameters[1])*cos(transformParameters[2])- sin(transformParameters[2])*cos(transformParameters[0]);
137 matrix[0][2]= cos(transformParameters[0])*sin(transformParameters[1])*cos(transformParameters[2])+sin(transformParameters[0])*sin(transformParameters[2]);
138 matrix[1][0]= cos(transformParameters[1])*sin(transformParameters[2]);
139 matrix[1][1]= sin(transformParameters[0])*sin(transformParameters[1])*sin(transformParameters[2])+cos(transformParameters[0])*cos(transformParameters[2]);
140 matrix[1][2]= cos(transformParameters[0])*sin(transformParameters[1])*sin(transformParameters[2])-sin(transformParameters[0])*cos(transformParameters[2]);
141 matrix[2][0]= -sin(transformParameters[1]);
142 matrix[2][1]= sin(transformParameters[0])*cos(transformParameters[1]);
143 matrix[2][2]= cos(transformParameters[0])*cos(transformParameters[1]);
145 matrix[0][3]=transformParameters[3];
146 matrix[1][3]=transformParameters[4];
147 matrix[2][3]=transformParameters[5];
156 inline itk::Matrix<double, 3, 3> GetRotationMatrix3D(itk::Array<double> rotationParameters)
158 itk::Matrix<double, 3, 3> matrix;
160 matrix[0][0]= cos(rotationParameters[1])*cos(rotationParameters[2]);
161 matrix[0][1]= sin(rotationParameters[0])*sin(rotationParameters[1])*cos(rotationParameters[2])+ sin(rotationParameters[2])*cos(rotationParameters[0]);
162 matrix[0][2]= -cos(rotationParameters[0])*sin(rotationParameters[1])*cos(rotationParameters[2])+sin(rotationParameters[0])*sin(rotationParameters[2]);
163 matrix[1][0]= -cos(rotationParameters[1])*sin(rotationParameters[2]);
164 matrix[1][1]= -sin(rotationParameters[0])*sin(rotationParameters[1])*sin(rotationParameters[2])+cos(rotationParameters[0])*cos(rotationParameters[2]);
165 matrix[1][2]= cos(rotationParameters[0])*sin(rotationParameters[1])*sin(rotationParameters[2])+sin(rotationParameters[0])*cos(rotationParameters[2]);
166 matrix[2][0]= sin(rotationParameters[1]);
167 matrix[2][1]= -sin(rotationParameters[0])*cos(rotationParameters[1]);
168 matrix[2][2]= cos(rotationParameters[0])*cos(rotationParameters[1]);
172 inline itk::Matrix<double, 2, 2> GetRotationMatrix2D(itk::Array<double> rotationParameters)
174 itk::Matrix<double, 2, 2> matrix;
175 matrix[0][0] = cos(rotationParameters[0]);
176 matrix[1][0] = sin(rotationParameters[0]);
177 matrix[0][1] = -matrix[1][0];
178 matrix[1][1] = matrix[0][0];
182 //========================================================================================
183 inline itk::Point<double, 3> GetRotatedPoint3D(itk::Array<double> rotationParameters, itk::Point<double, 3> input)
185 itk::Matrix<double, 3, 3> matrix = GetRotationMatrix3D(rotationParameters);
186 itk::Point<double, 3> output;
187 for (unsigned int i=0;i<3;i++)
190 for (unsigned int j=0;j<3;j++)
191 output[i]+=matrix(i,j)*input[j];
197 inline itk::Matrix<double, 4, 4> GetCenteredRotationMatrix3D(itk::Array<double> rotationParameters,itk::Point<double,3> centerOfRotation )
199 //rotational part is identical as affine matrix, translations change
200 itk::Array<double> parameters(6);
201 for(unsigned int i=0; i<3;i++) parameters[i]=rotationParameters[i];
202 for(unsigned int i=3; i<6;i++) parameters[i]=centerOfRotation[i-3];
203 itk::Matrix<double, 4, 4> matrix=GetForwardAffineMatrix3D(parameters);
205 //Get the rotation of the centerOfRotation
206 itk::Matrix<double,3,3> rotation = GetRotationalPartMatrix3D(matrix);
207 itk::Point<double,3> rotatedCenter=rotation*centerOfRotation; //GetRotatedPoint3D(rotationParameters, centerOfRotation);
209 //Substract this point to the translational part
210 matrix(0,3)-=rotatedCenter[0];
211 matrix(1,3)-=rotatedCenter[1];
212 matrix(2,3)-=rotatedCenter[2];
217 // inline itk::Matrix<double, 4, 4> GetComposedMatrix3D(itk::Matrix<double, 4, 4> firstAppliedTransform, itk::Matrix<double, 4, 4> secondAppliedTransform)
219 // itk::Matrix<double, 4, 4> matrix;
220 // for (unsigned int i=0;i<4;i++)
221 // for (unsigned int j=0;j<4;j++)
224 // for (unsigned int k=0;k<4;k++)
225 // matrix[i][j]+=firstAppliedTransform[i][k]*secondAppliedTransform[k][j];
231 //========================================================================================
232 inline itk::Matrix<double, 5, 5> ReadMatrix4D(std::string fileName)
236 openFileForReading(is, fileName);
237 std::vector<double> nb;
247 //copy it to the matrix
248 itk::Matrix<double, 5, 5> matrix;
249 unsigned int index=0;
250 for (unsigned int i=0;i<5;i++)
251 for (unsigned int j=0;j<5;j++)
252 matrix[i][j]=nb[index++];
256 inline itk::Matrix<double, 4, 4> ReadMatrix3D(std::string fileName)
260 openFileForReading(is, fileName);
261 std::vector<double> nb;
271 //copy it to the matrix
272 itk::Matrix<double, 4, 4> matrix;
273 unsigned int index=0;
274 for (unsigned int i=0;i<4;i++)
275 for (unsigned int j=0;j<4;j++)
276 matrix[i][j]=nb[index++];
280 inline itk::Matrix<double, 3, 3> ReadMatrix2D(std::string fileName)
284 openFileForReading(is, fileName);
285 std::vector<double> nb;
295 //copy it to the matrix
296 itk::Matrix<double, 3, 3> matrix;
297 unsigned int index=0;
298 for (unsigned int i=0;i<3;i++)
299 for (unsigned int j=0;j<3;j++)
300 matrix[i][j]=nb[index++];
304 template <unsigned int Dimension > inline itk::Matrix<double, Dimension+1 , Dimension+1> ReadMatrix(std::string fileName)
309 openFileForReading(is, fileName);
310 std::vector<double> nb;
320 //copy it to the matrix
321 itk::Matrix<double, Dimension+1, Dimension+1> matrix;
322 unsigned int index=0;
323 for (unsigned int i=0;i<Dimension+1;i++)
324 for (unsigned int j=0;j<Dimension+1;j++)
325 matrix[i][j]=nb[index++];
330 // template<> inline itk::Matrix<double, 3, 3> ReadMatrix<2> (std::string fileName)
332 // return ReadMatrix2D(fileName);
334 // template<> inline itk::Matrix<double, 4, 4> ReadMatrix<3> (std::string fileName)
336 // return ReadMatrix3D(fileName);
338 // template<> inline itk::Matrix<double, 5, 5> ReadMatrix<4> (std::string fileName)
340 // return ReadMatrix4D(fileName);
344 //========================================================================================
345 inline itk::Matrix<double, 4, 4> GetRotationalPartMatrix4D(itk::Matrix<double, 5, 5> input)
347 itk::Matrix<double,4,4> matrix;
348 matrix[0][0]= input[0][0];
349 matrix[0][1]= input[0][1];
350 matrix[0][2]= input[0][2];
351 matrix[0][3]= input[0][3];
352 matrix[1][0]= input[1][0];
353 matrix[1][1]= input[1][1];
354 matrix[1][2]= input[1][2];
355 matrix[1][3]= input[1][3];
356 matrix[2][0]= input[2][0];
357 matrix[2][1]= input[2][1];
358 matrix[2][2]= input[2][2];
359 matrix[2][2]= input[2][2];
360 matrix[2][3]= input[2][3];
361 matrix[3][0]= input[3][0];
362 matrix[3][1]= input[3][1];
363 matrix[3][2]= input[3][2];
364 matrix[3][2]= input[3][2];
365 matrix[3][3]= input[3][3];
371 inline itk::Matrix<double, 3, 3> GetRotationalPartMatrix3D(itk::Matrix<double, 4, 4> input)
373 itk::Matrix<double,3,3> matrix;
374 matrix[0][0]= input[0][0];
375 matrix[0][1]= input[0][1];
376 matrix[0][2]= input[0][2];
377 matrix[1][0]= input[1][0];
378 matrix[1][1]= input[1][1];
379 matrix[1][2]= input[1][2];
380 matrix[2][0]= input[2][0];
381 matrix[2][1]= input[2][1];
382 matrix[2][2]= input[2][2];
386 inline itk::Matrix<double, 2, 2> GetRotationalPartMatrix2D(itk::Matrix<double, 3, 3> input)
388 itk::Matrix<double,2,2> matrix;
389 matrix[0][0]= input[0][0];
390 matrix[0][1]= input[0][1];
391 matrix[0][2]= input[0][2];
392 matrix[1][0]= input[1][0];
393 matrix[1][1]= input[1][1];
397 inline itk::Matrix<double, 4, 4> GetRotationalPartMatrix(itk::Matrix<double, 5, 5> input)
399 return GetRotationalPartMatrix4D(input);
402 inline itk::Matrix<double, 3, 3> GetRotationalPartMatrix(itk::Matrix<double, 4, 4> input)
404 return GetRotationalPartMatrix3D(input);
407 inline itk::Matrix<double, 2, 2> GetRotationalPartMatrix(itk::Matrix<double, 3, 3> input)
409 return GetRotationalPartMatrix2D(input);
413 //========================================================================================
414 inline itk::Vector<double,4> GetTranslationPartMatrix4D(itk::Matrix<double, 5, 5> input)
416 itk::Vector<double,4> vec;
425 inline itk::Vector<double,3> GetTranslationPartMatrix3D(itk::Matrix<double, 4, 4> input)
427 itk::Vector<double,3> vec;
434 inline itk::Vector<double,2> GetTranslationPartMatrix2D(itk::Matrix<double, 3, 3> input)
436 itk::Vector<double,2> vec;
443 inline itk::Vector<double,4> GetTranslationPartMatrix(itk::Matrix<double, 5, 5> input)
446 return GetTranslationPartMatrix4D(input);
449 inline itk::Vector<double,3> GetTranslationPartMatrix(itk::Matrix<double, 4, 4> input)
452 return GetTranslationPartMatrix3D(input);
455 inline itk::Vector<double,2> GetTranslationPartMatrix(itk::Matrix<double, 3, 3> input)
458 return GetTranslationPartMatrix2D(input);
462 #endif //#define CLITKTRANSFORMUTILITIES_H