1 #ifndef CLITKTRANSFORMUTILITIES_H
2 #define CLITKTRANSFORMUTILITIES_H
7 #include "clitkIOCommon.h"
8 #include "clitkCommon.h"
13 //============================================================================
15 //============================================================================
16 itk::Matrix<double, 3, 3> GetForwardAffineMatrix2D(itk::Array<double> transformParameters);
17 itk::Matrix<double, 3, 3> GetBackwardAffineMatrix2D(itk::Array<double> transformParameters);
18 itk::Matrix<double, 4, 4> GetForwardAffineMatrix3D(itk::Array<double> transformParameters);
19 itk::Matrix<double, 4, 4> GetBackwardAffineMatrix3D(itk::Array<double> transformParameters);
20 itk::Matrix<double, 3, 3> GetRotationMatrix3D(itk::Array<double> rotationParameters);
21 itk::Point<double, 3> GetRotatedPoint3D(itk::Array<double> rotationParameters, itk::Point<double, 3> input);
22 itk::Matrix<double, 4, 4> GetCenteredRotationMatrix3D(itk::Array<double> rotationParameters,itk::Point<double,3> centerOfRotation);
23 // itk::Matrix<double, 4, 4> GetComposedMatrix3D(itk::Matrix<double, 4, 4> firstTransform, itk::Matrix<double, 4, 4> secondTransform);
25 itk::Matrix<double, 5, 5> ReadMatrix4D(std::string fileName);
26 itk::Matrix<double, 4, 4> ReadMatrix3D(std::string fileName);
27 itk::Matrix<double, 3, 3> ReadMatrix2D(std::string fileName);
28 template <unsigned int Dimension > itk::Matrix<double, Dimension+1 , Dimension+1> ReadMatrix(std::string fileName);
30 itk::Matrix<double, 3, 3> GetRotationalPartMatrix3D(itk::Matrix<double, 4, 4> input);
31 itk::Matrix<double, 3, 3> GetRotationalPartMatrix(itk::Matrix<double, 4, 4> input);
32 itk::Matrix<double, 2, 2> GetRotationalPartMatrix2D(itk::Matrix<double, 3, 3> input);
33 itk::Matrix<double, 2, 2> GetRotationalPartMatrix(itk::Matrix<double, 3, 3> input);
35 itk::Vector<double,3> GetTranslationPartMatrix3D(itk::Matrix<double, 4, 4> input);
36 itk::Vector<double,3> GetTranslationPartMatrix(itk::Matrix<double, 4, 4> input);
37 itk::Vector<double,2> GetTranslationPartMatrix2D(itk::Matrix<double, 3, 3> input);
38 itk::Vector<double,2> GetTranslationPartMatrix(itk::Matrix<double, 3, 3> input);
41 //============================================================================
42 //Inline functions definition in header file, otherwise linker errors
43 //============================================================================
45 //========================================================================================
46 inline itk::Matrix<double, 3, 3> GetForwardAffineMatrix2D(itk::Array<double> transformParameters)
48 itk::Matrix<double, 3, 3> matrix;
50 matrix[0][0]=cos(transformParameters[0]);
51 matrix[0][1]=-sin(transformParameters[0]);
52 matrix[1][0]=sin(transformParameters[0]);
53 matrix[1][1]=cos(transformParameters[0]);
55 matrix[0][2]=transformParameters[1];
56 matrix[1][2]=transformParameters[2];
64 inline itk::Matrix<double, 3, 3> GetBackwardAffineMatrix2D(itk::Array<double> transformParameters)
66 itk::Matrix<double, 3, 3> matrix;
68 matrix[0][0]=cos(transformParameters[0]);
69 matrix[0][1]=sin(transformParameters[0]);
70 matrix[1][0]=-sin(transformParameters[0]);
71 matrix[1][1]=cos(transformParameters[0]);
73 matrix[0][2]=transformParameters[1];
74 matrix[1][2]=transformParameters[2];
83 inline itk::Matrix<double, 4, 4> GetForwardAffineMatrix3D(itk::Array<double> transformParameters)
85 itk::Matrix<double, 4, 4> matrix;
87 matrix[0][0]= cos(transformParameters[1])*cos(transformParameters[2]);
88 matrix[0][1]= sin(transformParameters[0])*sin(transformParameters[1])*cos(transformParameters[2])+ sin(transformParameters[2])*cos(transformParameters[0]);
89 matrix[0][2]= -cos(transformParameters[0])*sin(transformParameters[1])*cos(transformParameters[2])+sin(transformParameters[0])*sin(transformParameters[2]);
90 matrix[1][0]= -cos(transformParameters[1])*sin(transformParameters[2]);
91 matrix[1][1]= -sin(transformParameters[0])*sin(transformParameters[1])*sin(transformParameters[2])+cos(transformParameters[0])*cos(transformParameters[2]);
92 matrix[1][2]= cos(transformParameters[0])*sin(transformParameters[1])*sin(transformParameters[2])+sin(transformParameters[0])*cos(transformParameters[2]);
93 matrix[2][0]= sin(transformParameters[1]);
94 matrix[2][1]= -sin(transformParameters[0])*cos(transformParameters[1]);
95 matrix[2][2]= cos(transformParameters[0])*cos(transformParameters[1]);
97 matrix[0][3]=transformParameters[3];
98 matrix[1][3]=transformParameters[4];
99 matrix[2][3]=transformParameters[5];
109 inline itk::Matrix<double, 4, 4> GetBackwardAffineMatrix3D(itk::Array<double> transformParameters)
111 itk::Matrix<double, 4, 4> matrix;
113 matrix[0][0]= cos(transformParameters[1])*cos(transformParameters[2]);
114 matrix[0][1]= sin(transformParameters[0])*sin(transformParameters[1])*cos(transformParameters[2])- sin(transformParameters[2])*cos(transformParameters[0]);
115 matrix[0][2]= cos(transformParameters[0])*sin(transformParameters[1])*cos(transformParameters[2])+sin(transformParameters[0])*sin(transformParameters[2]);
116 matrix[1][0]= cos(transformParameters[1])*sin(transformParameters[2]);
117 matrix[1][1]= sin(transformParameters[0])*sin(transformParameters[1])*sin(transformParameters[2])+cos(transformParameters[0])*cos(transformParameters[2]);
118 matrix[1][2]= cos(transformParameters[0])*sin(transformParameters[1])*sin(transformParameters[2])-sin(transformParameters[0])*cos(transformParameters[2]);
119 matrix[2][0]= -sin(transformParameters[1]);
120 matrix[2][1]= sin(transformParameters[0])*cos(transformParameters[1]);
121 matrix[2][2]= cos(transformParameters[0])*cos(transformParameters[1]);
123 matrix[0][3]=transformParameters[3];
124 matrix[1][3]=transformParameters[4];
125 matrix[2][3]=transformParameters[5];
134 inline itk::Matrix<double, 3, 3> GetRotationMatrix3D(itk::Array<double> rotationParameters)
136 itk::Matrix<double, 3, 3> matrix;
138 matrix[0][0]= cos(rotationParameters[1])*cos(rotationParameters[2]);
139 matrix[0][1]= sin(rotationParameters[0])*sin(rotationParameters[1])*cos(rotationParameters[2])+ sin(rotationParameters[2])*cos(rotationParameters[0]);
140 matrix[0][2]= -cos(rotationParameters[0])*sin(rotationParameters[1])*cos(rotationParameters[2])+sin(rotationParameters[0])*sin(rotationParameters[2]);
141 matrix[1][0]= -cos(rotationParameters[1])*sin(rotationParameters[2]);
142 matrix[1][1]= -sin(rotationParameters[0])*sin(rotationParameters[1])*sin(rotationParameters[2])+cos(rotationParameters[0])*cos(rotationParameters[2]);
143 matrix[1][2]= cos(rotationParameters[0])*sin(rotationParameters[1])*sin(rotationParameters[2])+sin(rotationParameters[0])*cos(rotationParameters[2]);
144 matrix[2][0]= sin(rotationParameters[1]);
145 matrix[2][1]= -sin(rotationParameters[0])*cos(rotationParameters[1]);
146 matrix[2][2]= cos(rotationParameters[0])*cos(rotationParameters[1]);
152 //========================================================================================
153 inline itk::Point<double, 3> GetRotatedPoint3D(itk::Array<double> rotationParameters, itk::Point<double, 3> input)
155 itk::Matrix<double, 3, 3> matrix = GetRotationMatrix3D(rotationParameters);
156 itk::Point<double, 3> output;
157 for (unsigned int i=0;i<3;i++)
160 for (unsigned int j=0;j<3;j++)
161 output[i]+=matrix(i,j)*input[j];
167 inline itk::Matrix<double, 4, 4> GetCenteredRotationMatrix3D(itk::Array<double> rotationParameters,itk::Point<double,3> centerOfRotation )
169 //rotational part is identical as affine matrix, translations change
170 itk::Array<double> parameters(6);
171 for(unsigned int i=0; i<3;i++) parameters[i]=rotationParameters[i];
172 for(unsigned int i=3; i<6;i++) parameters[i]=centerOfRotation[i-3];
173 itk::Matrix<double, 4, 4> matrix=GetForwardAffineMatrix3D(parameters);
175 //Get the rotation of the centerOfRotation
176 itk::Matrix<double,3,3> rotation = GetRotationalPartMatrix3D(matrix);
177 itk::Point<double,3> rotatedCenter=rotation*centerOfRotation; //GetRotatedPoint3D(rotationParameters, centerOfRotation);
179 //Substract this point to the translational part
180 matrix(0,3)-=rotatedCenter[0];
181 matrix(1,3)-=rotatedCenter[1];
182 matrix(2,3)-=rotatedCenter[2];
187 // inline itk::Matrix<double, 4, 4> GetComposedMatrix3D(itk::Matrix<double, 4, 4> firstAppliedTransform, itk::Matrix<double, 4, 4> secondAppliedTransform)
189 // itk::Matrix<double, 4, 4> matrix;
190 // for (unsigned int i=0;i<4;i++)
191 // for (unsigned int j=0;j<4;j++)
194 // for (unsigned int k=0;k<4;k++)
195 // matrix[i][j]+=firstAppliedTransform[i][k]*secondAppliedTransform[k][j];
201 //========================================================================================
202 inline itk::Matrix<double, 5, 5> ReadMatrix4D(std::string fileName)
206 openFileForReading(is, fileName);
207 std::vector<double> nb;
217 //copy it to the matrix
218 itk::Matrix<double, 5, 5> matrix;
219 unsigned int index=0;
220 for (unsigned int i=0;i<5;i++)
221 for (unsigned int j=0;j<5;j++)
222 matrix[i][j]=nb[index++];
226 inline itk::Matrix<double, 4, 4> ReadMatrix3D(std::string fileName)
230 openFileForReading(is, fileName);
231 std::vector<double> nb;
241 //copy it to the matrix
242 itk::Matrix<double, 4, 4> matrix;
243 unsigned int index=0;
244 for (unsigned int i=0;i<4;i++)
245 for (unsigned int j=0;j<4;j++)
246 matrix[i][j]=nb[index++];
250 inline itk::Matrix<double, 3, 3> ReadMatrix2D(std::string fileName)
254 openFileForReading(is, fileName);
255 std::vector<double> nb;
265 //copy it to the matrix
266 itk::Matrix<double, 3, 3> matrix;
267 unsigned int index=0;
268 for (unsigned int i=0;i<3;i++)
269 for (unsigned int j=0;j<3;j++)
270 matrix[i][j]=nb[index++];
274 template <unsigned int Dimension > inline itk::Matrix<double, Dimension+1 , Dimension+1> ReadMatrix(std::string fileName)
279 openFileForReading(is, fileName);
280 std::vector<double> nb;
290 //copy it to the matrix
291 itk::Matrix<double, Dimension+1, Dimension+1> matrix;
292 unsigned int index=0;
293 for (unsigned int i=0;i<Dimension+1;i++)
294 for (unsigned int j=0;j<Dimension+1;j++)
295 matrix[i][j]=nb[index++];
300 // template<> inline itk::Matrix<double, 3, 3> ReadMatrix<2> (std::string fileName)
302 // return ReadMatrix2D(fileName);
304 // template<> inline itk::Matrix<double, 4, 4> ReadMatrix<3> (std::string fileName)
306 // return ReadMatrix3D(fileName);
308 // template<> inline itk::Matrix<double, 5, 5> ReadMatrix<4> (std::string fileName)
310 // return ReadMatrix4D(fileName);
314 //========================================================================================
315 inline itk::Matrix<double, 4, 4> GetRotationalPartMatrix4D(itk::Matrix<double, 5, 5> input)
317 itk::Matrix<double,4,4> matrix;
318 matrix[0][0]= input[0][0];
319 matrix[0][1]= input[0][1];
320 matrix[0][2]= input[0][2];
321 matrix[0][3]= input[0][3];
322 matrix[1][0]= input[1][0];
323 matrix[1][1]= input[1][1];
324 matrix[1][2]= input[1][2];
325 matrix[1][3]= input[1][3];
326 matrix[2][0]= input[2][0];
327 matrix[2][1]= input[2][1];
328 matrix[2][2]= input[2][2];
329 matrix[2][2]= input[2][2];
330 matrix[2][3]= input[2][3];
331 matrix[3][0]= input[3][0];
332 matrix[3][1]= input[3][1];
333 matrix[3][2]= input[3][2];
334 matrix[3][2]= input[3][2];
335 matrix[3][3]= input[3][3];
341 inline itk::Matrix<double, 3, 3> GetRotationalPartMatrix3D(itk::Matrix<double, 4, 4> input)
343 itk::Matrix<double,3,3> matrix;
344 matrix[0][0]= input[0][0];
345 matrix[0][1]= input[0][1];
346 matrix[0][2]= input[0][2];
347 matrix[1][0]= input[1][0];
348 matrix[1][1]= input[1][1];
349 matrix[1][2]= input[1][2];
350 matrix[2][0]= input[2][0];
351 matrix[2][1]= input[2][1];
352 matrix[2][2]= input[2][2];
356 inline itk::Matrix<double, 2, 2> GetRotationalPartMatrix2D(itk::Matrix<double, 3, 3> input)
358 itk::Matrix<double,2,2> matrix;
359 matrix[0][0]= input[0][0];
360 matrix[0][1]= input[0][1];
361 matrix[0][2]= input[0][2];
362 matrix[1][0]= input[1][0];
363 matrix[1][1]= input[1][1];
367 inline itk::Matrix<double, 4, 4> GetRotationalPartMatrix(itk::Matrix<double, 5, 5> input)
369 return GetRotationalPartMatrix4D(input);
372 inline itk::Matrix<double, 3, 3> GetRotationalPartMatrix(itk::Matrix<double, 4, 4> input)
374 return GetRotationalPartMatrix3D(input);
377 inline itk::Matrix<double, 2, 2> GetRotationalPartMatrix(itk::Matrix<double, 3, 3> input)
379 return GetRotationalPartMatrix2D(input);
383 //========================================================================================
384 inline itk::Vector<double,4> GetTranslationPartMatrix4D(itk::Matrix<double, 5, 5> input)
386 itk::Vector<double,4> vec;
395 inline itk::Vector<double,3> GetTranslationPartMatrix3D(itk::Matrix<double, 4, 4> input)
397 itk::Vector<double,3> vec;
404 inline itk::Vector<double,2> GetTranslationPartMatrix2D(itk::Matrix<double, 3, 3> input)
406 itk::Vector<double,2> vec;
413 inline itk::Vector<double,4> GetTranslationPartMatrix(itk::Matrix<double, 5, 5> input)
416 return GetTranslationPartMatrix4D(input);
419 inline itk::Vector<double,3> GetTranslationPartMatrix(itk::Matrix<double, 4, 4> input)
422 return GetTranslationPartMatrix3D(input);
425 inline itk::Vector<double,2> GetTranslationPartMatrix(itk::Matrix<double, 3, 3> input)
428 return GetTranslationPartMatrix2D(input);
432 #endif //#define CLITKTRANSFORMUTILITIES_H