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://oncora1.lyon.fnclcc.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 ======================================================================-====*/
18 #ifndef CLITKTRANSFORMUTILITIES_H
19 #define CLITKTRANSFORMUTILITIES_H
20 #include "itkMatrix.h"
23 #include "clitkImageCommon.h"
24 #include "clitkCommon.h"
29 //============================================================================
31 //============================================================================
32 itk::Matrix<double, 3, 3> GetForwardAffineMatrix2D(itk::Array<double> transformParameters);
33 itk::Matrix<double, 3, 3> GetBackwardAffineMatrix2D(itk::Array<double> transformParameters);
34 itk::Matrix<double, 4, 4> GetForwardAffineMatrix3D(itk::Array<double> transformParameters);
35 itk::Matrix<double, 4, 4> GetBackwardAffineMatrix3D(itk::Array<double> transformParameters);
36 itk::Matrix<double, 3, 3> GetRotationMatrix3D(itk::Array<double> rotationParameters);
37 itk::Point<double, 3> GetRotatedPoint3D(itk::Array<double> rotationParameters, itk::Point<double, 3> input);
38 itk::Matrix<double, 4, 4> GetCenteredRotationMatrix3D(itk::Array<double> rotationParameters,itk::Point<double,3> centerOfRotation);
39 // itk::Matrix<double, 4, 4> GetComposedMatrix3D(itk::Matrix<double, 4, 4> firstTransform, itk::Matrix<double, 4, 4> secondTransform);
41 itk::Matrix<double, 5, 5> ReadMatrix4D(std::string fileName);
42 itk::Matrix<double, 4, 4> ReadMatrix3D(std::string fileName);
43 itk::Matrix<double, 3, 3> ReadMatrix2D(std::string fileName);
44 template <unsigned int Dimension > itk::Matrix<double, Dimension+1 , Dimension+1> ReadMatrix(std::string fileName);
46 itk::Matrix<double, 3, 3> GetRotationalPartMatrix3D(itk::Matrix<double, 4, 4> input);
47 itk::Matrix<double, 3, 3> GetRotationalPartMatrix(itk::Matrix<double, 4, 4> input);
48 itk::Matrix<double, 2, 2> GetRotationalPartMatrix2D(itk::Matrix<double, 3, 3> input);
49 itk::Matrix<double, 2, 2> GetRotationalPartMatrix(itk::Matrix<double, 3, 3> input);
51 itk::Vector<double,3> GetTranslationPartMatrix3D(itk::Matrix<double, 4, 4> input);
52 itk::Vector<double,3> GetTranslationPartMatrix(itk::Matrix<double, 4, 4> input);
53 itk::Vector<double,2> GetTranslationPartMatrix2D(itk::Matrix<double, 3, 3> input);
54 itk::Vector<double,2> GetTranslationPartMatrix(itk::Matrix<double, 3, 3> input);
57 //============================================================================
58 //Inline functions definition in header file, otherwise linker errors
59 //============================================================================
61 //========================================================================================
62 inline itk::Matrix<double, 3, 3> GetForwardAffineMatrix2D(itk::Array<double> transformParameters)
64 itk::Matrix<double, 3, 3> matrix;
66 matrix[0][0]=cos(transformParameters[0]);
67 matrix[0][1]=-sin(transformParameters[0]);
68 matrix[1][0]=sin(transformParameters[0]);
69 matrix[1][1]=cos(transformParameters[0]);
71 matrix[0][2]=transformParameters[1];
72 matrix[1][2]=transformParameters[2];
80 inline itk::Matrix<double, 3, 3> GetBackwardAffineMatrix2D(itk::Array<double> transformParameters)
82 itk::Matrix<double, 3, 3> matrix;
84 matrix[0][0]=cos(transformParameters[0]);
85 matrix[0][1]=sin(transformParameters[0]);
86 matrix[1][0]=-sin(transformParameters[0]);
87 matrix[1][1]=cos(transformParameters[0]);
89 matrix[0][2]=transformParameters[1];
90 matrix[1][2]=transformParameters[2];
99 inline itk::Matrix<double, 4, 4> GetForwardAffineMatrix3D(itk::Array<double> transformParameters)
101 itk::Matrix<double, 4, 4> matrix;
103 matrix[0][0]= cos(transformParameters[1])*cos(transformParameters[2]);
104 matrix[0][1]= sin(transformParameters[0])*sin(transformParameters[1])*cos(transformParameters[2])+ sin(transformParameters[2])*cos(transformParameters[0]);
105 matrix[0][2]= -cos(transformParameters[0])*sin(transformParameters[1])*cos(transformParameters[2])+sin(transformParameters[0])*sin(transformParameters[2]);
106 matrix[1][0]= -cos(transformParameters[1])*sin(transformParameters[2]);
107 matrix[1][1]= -sin(transformParameters[0])*sin(transformParameters[1])*sin(transformParameters[2])+cos(transformParameters[0])*cos(transformParameters[2]);
108 matrix[1][2]= cos(transformParameters[0])*sin(transformParameters[1])*sin(transformParameters[2])+sin(transformParameters[0])*cos(transformParameters[2]);
109 matrix[2][0]= sin(transformParameters[1]);
110 matrix[2][1]= -sin(transformParameters[0])*cos(transformParameters[1]);
111 matrix[2][2]= cos(transformParameters[0])*cos(transformParameters[1]);
113 matrix[0][3]=transformParameters[3];
114 matrix[1][3]=transformParameters[4];
115 matrix[2][3]=transformParameters[5];
125 inline itk::Matrix<double, 4, 4> GetBackwardAffineMatrix3D(itk::Array<double> transformParameters)
127 itk::Matrix<double, 4, 4> matrix;
129 matrix[0][0]= cos(transformParameters[1])*cos(transformParameters[2]);
130 matrix[0][1]= sin(transformParameters[0])*sin(transformParameters[1])*cos(transformParameters[2])- sin(transformParameters[2])*cos(transformParameters[0]);
131 matrix[0][2]= cos(transformParameters[0])*sin(transformParameters[1])*cos(transformParameters[2])+sin(transformParameters[0])*sin(transformParameters[2]);
132 matrix[1][0]= cos(transformParameters[1])*sin(transformParameters[2]);
133 matrix[1][1]= sin(transformParameters[0])*sin(transformParameters[1])*sin(transformParameters[2])+cos(transformParameters[0])*cos(transformParameters[2]);
134 matrix[1][2]= cos(transformParameters[0])*sin(transformParameters[1])*sin(transformParameters[2])-sin(transformParameters[0])*cos(transformParameters[2]);
135 matrix[2][0]= -sin(transformParameters[1]);
136 matrix[2][1]= sin(transformParameters[0])*cos(transformParameters[1]);
137 matrix[2][2]= cos(transformParameters[0])*cos(transformParameters[1]);
139 matrix[0][3]=transformParameters[3];
140 matrix[1][3]=transformParameters[4];
141 matrix[2][3]=transformParameters[5];
150 inline itk::Matrix<double, 3, 3> GetRotationMatrix3D(itk::Array<double> rotationParameters)
152 itk::Matrix<double, 3, 3> matrix;
154 matrix[0][0]= cos(rotationParameters[1])*cos(rotationParameters[2]);
155 matrix[0][1]= sin(rotationParameters[0])*sin(rotationParameters[1])*cos(rotationParameters[2])+ sin(rotationParameters[2])*cos(rotationParameters[0]);
156 matrix[0][2]= -cos(rotationParameters[0])*sin(rotationParameters[1])*cos(rotationParameters[2])+sin(rotationParameters[0])*sin(rotationParameters[2]);
157 matrix[1][0]= -cos(rotationParameters[1])*sin(rotationParameters[2]);
158 matrix[1][1]= -sin(rotationParameters[0])*sin(rotationParameters[1])*sin(rotationParameters[2])+cos(rotationParameters[0])*cos(rotationParameters[2]);
159 matrix[1][2]= cos(rotationParameters[0])*sin(rotationParameters[1])*sin(rotationParameters[2])+sin(rotationParameters[0])*cos(rotationParameters[2]);
160 matrix[2][0]= sin(rotationParameters[1]);
161 matrix[2][1]= -sin(rotationParameters[0])*cos(rotationParameters[1]);
162 matrix[2][2]= cos(rotationParameters[0])*cos(rotationParameters[1]);
168 //========================================================================================
169 inline itk::Point<double, 3> GetRotatedPoint3D(itk::Array<double> rotationParameters, itk::Point<double, 3> input)
171 itk::Matrix<double, 3, 3> matrix = GetRotationMatrix3D(rotationParameters);
172 itk::Point<double, 3> output;
173 for (unsigned int i=0;i<3;i++)
176 for (unsigned int j=0;j<3;j++)
177 output[i]+=matrix(i,j)*input[j];
183 inline itk::Matrix<double, 4, 4> GetCenteredRotationMatrix3D(itk::Array<double> rotationParameters,itk::Point<double,3> centerOfRotation )
185 //rotational part is identical as affine matrix, translations change
186 itk::Array<double> parameters(6);
187 for(unsigned int i=0; i<3;i++) parameters[i]=rotationParameters[i];
188 for(unsigned int i=3; i<6;i++) parameters[i]=centerOfRotation[i-3];
189 itk::Matrix<double, 4, 4> matrix=GetForwardAffineMatrix3D(parameters);
191 //Get the rotation of the centerOfRotation
192 itk::Matrix<double,3,3> rotation = GetRotationalPartMatrix3D(matrix);
193 itk::Point<double,3> rotatedCenter=rotation*centerOfRotation; //GetRotatedPoint3D(rotationParameters, centerOfRotation);
195 //Substract this point to the translational part
196 matrix(0,3)-=rotatedCenter[0];
197 matrix(1,3)-=rotatedCenter[1];
198 matrix(2,3)-=rotatedCenter[2];
203 // inline itk::Matrix<double, 4, 4> GetComposedMatrix3D(itk::Matrix<double, 4, 4> firstAppliedTransform, itk::Matrix<double, 4, 4> secondAppliedTransform)
205 // itk::Matrix<double, 4, 4> matrix;
206 // for (unsigned int i=0;i<4;i++)
207 // for (unsigned int j=0;j<4;j++)
210 // for (unsigned int k=0;k<4;k++)
211 // matrix[i][j]+=firstAppliedTransform[i][k]*secondAppliedTransform[k][j];
217 //========================================================================================
218 inline itk::Matrix<double, 5, 5> ReadMatrix4D(std::string fileName)
222 openFileForReading(is, fileName);
223 std::vector<double> nb;
233 //copy it to the matrix
234 itk::Matrix<double, 5, 5> matrix;
235 unsigned int index=0;
236 for (unsigned int i=0;i<5;i++)
237 for (unsigned int j=0;j<5;j++)
238 matrix[i][j]=nb[index++];
242 inline itk::Matrix<double, 4, 4> ReadMatrix3D(std::string fileName)
246 openFileForReading(is, fileName);
247 std::vector<double> nb;
257 //copy it to the matrix
258 itk::Matrix<double, 4, 4> matrix;
259 unsigned int index=0;
260 for (unsigned int i=0;i<4;i++)
261 for (unsigned int j=0;j<4;j++)
262 matrix[i][j]=nb[index++];
266 inline itk::Matrix<double, 3, 3> ReadMatrix2D(std::string fileName)
270 openFileForReading(is, fileName);
271 std::vector<double> nb;
281 //copy it to the matrix
282 itk::Matrix<double, 3, 3> matrix;
283 unsigned int index=0;
284 for (unsigned int i=0;i<3;i++)
285 for (unsigned int j=0;j<3;j++)
286 matrix[i][j]=nb[index++];
290 template <unsigned int Dimension > inline itk::Matrix<double, Dimension+1 , Dimension+1> ReadMatrix(std::string fileName)
295 openFileForReading(is, fileName);
296 std::vector<double> nb;
306 //copy it to the matrix
307 itk::Matrix<double, Dimension+1, Dimension+1> matrix;
308 unsigned int index=0;
309 for (unsigned int i=0;i<Dimension+1;i++)
310 for (unsigned int j=0;j<Dimension+1;j++)
311 matrix[i][j]=nb[index++];
316 // template<> inline itk::Matrix<double, 3, 3> ReadMatrix<2> (std::string fileName)
318 // return ReadMatrix2D(fileName);
320 // template<> inline itk::Matrix<double, 4, 4> ReadMatrix<3> (std::string fileName)
322 // return ReadMatrix3D(fileName);
324 // template<> inline itk::Matrix<double, 5, 5> ReadMatrix<4> (std::string fileName)
326 // return ReadMatrix4D(fileName);
330 //========================================================================================
331 inline itk::Matrix<double, 4, 4> GetRotationalPartMatrix4D(itk::Matrix<double, 5, 5> input)
333 itk::Matrix<double,4,4> matrix;
334 matrix[0][0]= input[0][0];
335 matrix[0][1]= input[0][1];
336 matrix[0][2]= input[0][2];
337 matrix[0][3]= input[0][3];
338 matrix[1][0]= input[1][0];
339 matrix[1][1]= input[1][1];
340 matrix[1][2]= input[1][2];
341 matrix[1][3]= input[1][3];
342 matrix[2][0]= input[2][0];
343 matrix[2][1]= input[2][1];
344 matrix[2][2]= input[2][2];
345 matrix[2][2]= input[2][2];
346 matrix[2][3]= input[2][3];
347 matrix[3][0]= input[3][0];
348 matrix[3][1]= input[3][1];
349 matrix[3][2]= input[3][2];
350 matrix[3][2]= input[3][2];
351 matrix[3][3]= input[3][3];
357 inline itk::Matrix<double, 3, 3> GetRotationalPartMatrix3D(itk::Matrix<double, 4, 4> input)
359 itk::Matrix<double,3,3> matrix;
360 matrix[0][0]= input[0][0];
361 matrix[0][1]= input[0][1];
362 matrix[0][2]= input[0][2];
363 matrix[1][0]= input[1][0];
364 matrix[1][1]= input[1][1];
365 matrix[1][2]= input[1][2];
366 matrix[2][0]= input[2][0];
367 matrix[2][1]= input[2][1];
368 matrix[2][2]= input[2][2];
372 inline itk::Matrix<double, 2, 2> GetRotationalPartMatrix2D(itk::Matrix<double, 3, 3> input)
374 itk::Matrix<double,2,2> 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];
383 inline itk::Matrix<double, 4, 4> GetRotationalPartMatrix(itk::Matrix<double, 5, 5> input)
385 return GetRotationalPartMatrix4D(input);
388 inline itk::Matrix<double, 3, 3> GetRotationalPartMatrix(itk::Matrix<double, 4, 4> input)
390 return GetRotationalPartMatrix3D(input);
393 inline itk::Matrix<double, 2, 2> GetRotationalPartMatrix(itk::Matrix<double, 3, 3> input)
395 return GetRotationalPartMatrix2D(input);
399 //========================================================================================
400 inline itk::Vector<double,4> GetTranslationPartMatrix4D(itk::Matrix<double, 5, 5> input)
402 itk::Vector<double,4> vec;
411 inline itk::Vector<double,3> GetTranslationPartMatrix3D(itk::Matrix<double, 4, 4> input)
413 itk::Vector<double,3> vec;
420 inline itk::Vector<double,2> GetTranslationPartMatrix2D(itk::Matrix<double, 3, 3> input)
422 itk::Vector<double,2> vec;
429 inline itk::Vector<double,4> GetTranslationPartMatrix(itk::Matrix<double, 5, 5> input)
432 return GetTranslationPartMatrix4D(input);
435 inline itk::Vector<double,3> GetTranslationPartMatrix(itk::Matrix<double, 4, 4> input)
438 return GetTranslationPartMatrix3D(input);
441 inline itk::Vector<double,2> GetTranslationPartMatrix(itk::Matrix<double, 3, 3> input)
444 return GetTranslationPartMatrix2D(input);
448 #endif //#define CLITKTRANSFORMUTILITIES_H