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 itk::Matrix<double, 4, 4> ReadMatrix3D(std::string fileName);
260 inline vtkMatrix4x4* ReadVTKMatrix3D(std::string fileName) {
263 openFileForReading(is, fileName);
264 std::vector<double> nb;
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++]);
283 inline itk::Matrix<double, 3, 3> ReadMatrix2D(std::string fileName)
287 openFileForReading(is, fileName);
288 std::vector<double> nb;
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++];
307 template <unsigned int Dimension > inline itk::Matrix<double, Dimension+1 , Dimension+1> ReadMatrix(std::string fileName)
312 openFileForReading(is, fileName);
313 std::vector<double> nb;
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++];
333 // template<> inline itk::Matrix<double, 3, 3> ReadMatrix<2> (std::string fileName)
335 // return ReadMatrix2D(fileName);
337 // template<> inline itk::Matrix<double, 4, 4> ReadMatrix<3> (std::string fileName)
339 // return ReadMatrix3D(fileName);
341 // template<> inline itk::Matrix<double, 5, 5> ReadMatrix<4> (std::string fileName)
343 // return ReadMatrix4D(fileName);
347 //========================================================================================
348 inline itk::Matrix<double, 4, 4> GetRotationalPartMatrix4D(itk::Matrix<double, 5, 5> input)
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];
374 inline itk::Matrix<double, 3, 3> GetRotationalPartMatrix3D(itk::Matrix<double, 4, 4> input)
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];
389 inline itk::Matrix<double, 2, 2> GetRotationalPartMatrix2D(itk::Matrix<double, 3, 3> input)
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];
400 inline itk::Matrix<double, 4, 4> GetRotationalPartMatrix(itk::Matrix<double, 5, 5> input)
402 return GetRotationalPartMatrix4D(input);
405 inline itk::Matrix<double, 3, 3> GetRotationalPartMatrix(itk::Matrix<double, 4, 4> input)
407 return GetRotationalPartMatrix3D(input);
410 inline itk::Matrix<double, 2, 2> GetRotationalPartMatrix(itk::Matrix<double, 3, 3> input)
412 return GetRotationalPartMatrix2D(input);
416 //========================================================================================
417 inline itk::Vector<double,4> GetTranslationPartMatrix4D(itk::Matrix<double, 5, 5> input)
419 itk::Vector<double,4> vec;
428 inline itk::Vector<double,3> GetTranslationPartMatrix3D(itk::Matrix<double, 4, 4> input)
430 itk::Vector<double,3> vec;
437 inline itk::Vector<double,2> GetTranslationPartMatrix2D(itk::Matrix<double, 3, 3> input)
439 itk::Vector<double,2> vec;
446 inline itk::Vector<double,4> GetTranslationPartMatrix(itk::Matrix<double, 5, 5> input)
449 return GetTranslationPartMatrix4D(input);
452 inline itk::Vector<double,3> GetTranslationPartMatrix(itk::Matrix<double, 4, 4> input)
455 return GetTranslationPartMatrix3D(input);
458 inline itk::Vector<double,2> GetTranslationPartMatrix(itk::Matrix<double, 3, 3> input)
461 return GetTranslationPartMatrix2D(input);
465 #endif //#define CLITKTRANSFORMUTILITIES_H