]> Creatis software - clitk.git/blob - common/clitkTransformUtilities.h
itk4.3 compatibility
[clitk.git] / common / clitkTransformUtilities.h
1 /*=========================================================================
2   Program:   vv                     http://www.creatis.insa-lyon.fr/rio/vv
3
4   Authors belong to: 
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
8
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.
12
13   It is distributed under dual licence
14
15   - BSD        See included LICENSE.txt file
16   - CeCILL-B   http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
17 ===========================================================================**/
18
19 #ifndef CLITKTRANSFORMUTILITIES_H
20 #define CLITKTRANSFORMUTILITIES_H
21 #include "itkMatrix.h"
22 #include "itkArray.h"
23 #include "itkPoint.h"
24 #include "clitkImageCommon.h"
25 #include "clitkCommon.h"
26  
27  
28 namespace clitk
29 {
30   //============================================================================
31   //Declarations
32   //============================================================================
33   itk::Matrix<double, 3, 3> GetForwardAffineMatrix2D(itk::Array<double> transformParameters);
34   itk::Matrix<double, 4, 4> GetForwardAffineMatrix3D(itk::Array<double> transformParameters);
35   template <unsigned int Dimension > itk::Matrix<double, Dimension+1, Dimension+1> GetForwardAffineMatrix(itk::Array<double> transformParameters);
36
37   itk::Matrix<double, 3, 3> GetBackwardAffineMatrix2D(itk::Array<double> transformParameters);
38   itk::Matrix<double, 4, 4> GetBackwardAffineMatrix3D(itk::Array<double> transformParameters);
39   template <unsigned int Dimension > itk::Matrix<double, Dimension+1, Dimension+1> GetBackwardAffineMatrix(itk::Array<double> transformParameters);
40
41   itk::Matrix<double, 3, 3> GetRotationMatrix3D(itk::Array<double> rotationParameters);
42   itk::Matrix<double, 2, 2> GetRotationMatrix2D(itk::Array<double> rotationParameters);
43   template <unsigned int Dimension> itk::Matrix<double, Dimension, Dimension> GetRotationMatrix(itk::Array<double> rotationParameters);
44   itk::Point<double, 3> GetRotatedPoint3D(itk::Array<double> rotationParameters, itk::Point<double, 3> input);
45   itk::Matrix<double, 4, 4> GetCenteredRotationMatrix3D(itk::Array<double> rotationParameters,itk::Point<double,3> centerOfRotation);
46   //   itk::Matrix<double, 4, 4> GetComposedMatrix3D(itk::Matrix<double, 4, 4> firstTransform, itk::Matrix<double, 4, 4> secondTransform);
47   
48   itk::Matrix<double, 5, 5> ReadMatrix4D(std::string fileName);
49   itk::Matrix<double, 4, 4> ReadMatrix3D(std::string fileName);
50   itk::Matrix<double, 3, 3> ReadMatrix2D(std::string fileName);
51   template <unsigned int Dimension > itk::Matrix<double, Dimension+1 , Dimension+1> ReadMatrix(std::string fileName);
52    
53   itk::Matrix<double, 3, 3> GetRotationalPartMatrix3D(itk::Matrix<double, 4, 4> input);
54   itk::Matrix<double, 3, 3> GetRotationalPartMatrix(itk::Matrix<double, 4, 4> input);
55   itk::Matrix<double, 2, 2> GetRotationalPartMatrix2D(itk::Matrix<double, 3, 3> input);
56   itk::Matrix<double, 2, 2> GetRotationalPartMatrix(itk::Matrix<double, 3, 3> input);
57      
58   itk::Vector<double,3> GetTranslationPartMatrix3D(itk::Matrix<double, 4, 4> input);
59   itk::Vector<double,3> GetTranslationPartMatrix(itk::Matrix<double, 4, 4> input);
60   itk::Vector<double,2> GetTranslationPartMatrix2D(itk::Matrix<double, 3, 3> input);
61   itk::Vector<double,2> GetTranslationPartMatrix(itk::Matrix<double, 3, 3> input);
62     
63  
64   //============================================================================
65   //Inline functions definition in header file, otherwise linker errors
66   //============================================================================
67  
68   //========================================================================================
69   inline  itk::Matrix<double, 3, 3> GetForwardAffineMatrix2D(itk::Array<double> transformParameters)
70   {
71     itk::Matrix<double, 3, 3> matrix;
72     //rotation part
73     matrix[0][0]=cos(transformParameters[0]);
74     matrix[0][1]=-sin(transformParameters[0]);
75     matrix[1][0]=sin(transformParameters[0]);
76     matrix[1][1]=cos(transformParameters[0]);
77     //translation part
78     matrix[0][2]=transformParameters[1];
79     matrix[1][2]=transformParameters[2];
80     //homogenize
81     matrix[2][0]=0.;
82     matrix[2][1]=0.;
83     matrix[2][2]=1.;
84     return matrix;
85   }
86   
87   inline  itk::Matrix<double, 4, 4> GetForwardAffineMatrix3D(itk::Array<double> transformParameters)
88   {
89     itk::Matrix<double, 4, 4> matrix;
90     //rotational part
91     matrix[0][0]= cos(transformParameters[1])*cos(transformParameters[2]);
92     matrix[0][1]= sin(transformParameters[0])*sin(transformParameters[1])*cos(transformParameters[2])+ sin(transformParameters[2])*cos(transformParameters[0]);
93     matrix[0][2]= -cos(transformParameters[0])*sin(transformParameters[1])*cos(transformParameters[2])+sin(transformParameters[0])*sin(transformParameters[2]);
94     matrix[1][0]= -cos(transformParameters[1])*sin(transformParameters[2]);
95     matrix[1][1]= -sin(transformParameters[0])*sin(transformParameters[1])*sin(transformParameters[2])+cos(transformParameters[0])*cos(transformParameters[2]);
96     matrix[1][2]= cos(transformParameters[0])*sin(transformParameters[1])*sin(transformParameters[2])+sin(transformParameters[0])*cos(transformParameters[2]);
97     matrix[2][0]= sin(transformParameters[1]);
98     matrix[2][1]= -sin(transformParameters[0])*cos(transformParameters[1]);
99     matrix[2][2]= cos(transformParameters[0])*cos(transformParameters[1]);
100     //translational part
101     matrix[0][3]=transformParameters[3];
102     matrix[1][3]=transformParameters[4];
103     matrix[2][3]=transformParameters[5];
104     //homogenize
105     matrix[3][0]=0.;
106     matrix[3][1]=0.;
107     matrix[3][2]=0.;  
108     matrix[3][3]=1.;
109     return matrix;
110   }
111  
112  
113   inline  itk::Matrix<double, 3, 3> GetBackwardAffineMatrix2D(itk::Array<double> transformParameters)
114   {
115     itk::Matrix<double, 3, 3> matrix;
116     //rotation part
117     matrix[0][0]=cos(transformParameters[0]);
118     matrix[0][1]=sin(transformParameters[0]);
119     matrix[1][0]=-sin(transformParameters[0]);
120     matrix[1][1]=cos(transformParameters[0]);
121     //translation part
122     matrix[0][2]=transformParameters[1];
123     matrix[1][2]=transformParameters[2];
124     //homogenize
125     matrix[2][0]=0.;
126     matrix[2][1]=0.;
127     matrix[2][2]=1.;
128     return matrix;
129   }
130
131
132   inline  itk::Matrix<double, 4, 4> GetBackwardAffineMatrix3D(itk::Array<double> transformParameters)
133   {
134     itk::Matrix<double, 4, 4> matrix;
135     //rotational part
136     matrix[0][0]= cos(transformParameters[1])*cos(transformParameters[2]);
137     matrix[0][1]= sin(transformParameters[0])*sin(transformParameters[1])*cos(transformParameters[2])- sin(transformParameters[2])*cos(transformParameters[0]);
138     matrix[0][2]= cos(transformParameters[0])*sin(transformParameters[1])*cos(transformParameters[2])+sin(transformParameters[0])*sin(transformParameters[2]);
139     matrix[1][0]= cos(transformParameters[1])*sin(transformParameters[2]);
140     matrix[1][1]= sin(transformParameters[0])*sin(transformParameters[1])*sin(transformParameters[2])+cos(transformParameters[0])*cos(transformParameters[2]);
141     matrix[1][2]= cos(transformParameters[0])*sin(transformParameters[1])*sin(transformParameters[2])-sin(transformParameters[0])*cos(transformParameters[2]);
142     matrix[2][0]= -sin(transformParameters[1]);
143     matrix[2][1]= sin(transformParameters[0])*cos(transformParameters[1]);
144     matrix[2][2]= cos(transformParameters[0])*cos(transformParameters[1]);
145     //translational part
146     matrix[0][3]=transformParameters[3];
147     matrix[1][3]=transformParameters[4];
148     matrix[2][3]=transformParameters[5];
149     //homogenize
150     matrix[3][0]=0.;
151     matrix[3][1]=0.;
152     matrix[3][2]=0.;  
153     matrix[3][3]=1.;
154     return matrix;
155   }
156   
157   inline itk::Matrix<double, 3, 3> GetRotationMatrix3D(itk::Array<double> rotationParameters)
158   {
159     itk::Matrix<double, 3, 3> matrix;
160     //rotational part
161     matrix[0][0]= cos(rotationParameters[1])*cos(rotationParameters[2]);
162     matrix[0][1]= sin(rotationParameters[0])*sin(rotationParameters[1])*cos(rotationParameters[2])+ sin(rotationParameters[2])*cos(rotationParameters[0]);
163     matrix[0][2]= -cos(rotationParameters[0])*sin(rotationParameters[1])*cos(rotationParameters[2])+sin(rotationParameters[0])*sin(rotationParameters[2]);
164     matrix[1][0]= -cos(rotationParameters[1])*sin(rotationParameters[2]);
165     matrix[1][1]= -sin(rotationParameters[0])*sin(rotationParameters[1])*sin(rotationParameters[2])+cos(rotationParameters[0])*cos(rotationParameters[2]);
166     matrix[1][2]= cos(rotationParameters[0])*sin(rotationParameters[1])*sin(rotationParameters[2])+sin(rotationParameters[0])*cos(rotationParameters[2]);
167     matrix[2][0]= sin(rotationParameters[1]);
168     matrix[2][1]= -sin(rotationParameters[0])*cos(rotationParameters[1]);
169     matrix[2][2]= cos(rotationParameters[0])*cos(rotationParameters[1]);
170     return matrix;
171   }
172
173   inline itk::Matrix<double, 2, 2> GetRotationMatrix2D(itk::Array<double> rotationParameters)
174   {
175     itk::Matrix<double, 2, 2> matrix;
176     matrix[0][0] = cos(rotationParameters[0]);
177     matrix[1][0] = sin(rotationParameters[0]);
178     matrix[0][1] = -matrix[1][0];
179     matrix[1][1] = matrix[0][0];
180     return matrix;
181   }
182
183   //========================================================================================
184   inline  itk::Point<double, 3> GetRotatedPoint3D(itk::Array<double> rotationParameters, itk::Point<double, 3> input)
185   {
186     itk::Matrix<double, 3, 3> matrix = GetRotationMatrix3D(rotationParameters);
187     itk::Point<double, 3> output;
188     for (unsigned int i=0;i<3;i++)
189       {
190         output[i]=0.0;
191         for (unsigned int j=0;j<3;j++)
192           output[i]+=matrix(i,j)*input[j];
193       }
194     return output;
195   }
196    
197    
198   inline itk::Matrix<double, 4, 4> GetCenteredRotationMatrix3D(itk::Array<double> rotationParameters,itk::Point<double,3> centerOfRotation )
199   {
200     //rotational part is identical as affine matrix, translations change
201     itk::Array<double> parameters(6);
202     for(unsigned int i=0; i<3;i++) parameters[i]=rotationParameters[i];
203     for(unsigned int i=3; i<6;i++) parameters[i]=centerOfRotation[i-3];
204     itk::Matrix<double, 4, 4> matrix=GetForwardAffineMatrix3D(parameters);
205      
206     //Get the rotation of the centerOfRotation
207     itk::Matrix<double,3,3> rotation = GetRotationalPartMatrix3D(matrix);
208     itk::Point<double,3> rotatedCenter=rotation*centerOfRotation; //GetRotatedPoint3D(rotationParameters, centerOfRotation);
209      
210     //Substract this point to the translational part 
211     matrix(0,3)-=rotatedCenter[0];
212     matrix(1,3)-=rotatedCenter[1];
213     matrix(2,3)-=rotatedCenter[2];
214     return matrix;
215   }
216  
217     
218   //   inline  itk::Matrix<double, 4, 4> GetComposedMatrix3D(itk::Matrix<double, 4, 4> firstAppliedTransform, itk::Matrix<double, 4, 4> secondAppliedTransform)
219   //   {
220   //     itk::Matrix<double, 4, 4> matrix;
221   //     for (unsigned int i=0;i<4;i++)
222   //       for (unsigned int j=0;j<4;j++)
223   //    {
224   //      matrix[i][j]=0.0;
225   //      for (unsigned int k=0;k<4;k++)
226   //        matrix[i][j]+=firstAppliedTransform[i][k]*secondAppliedTransform[k][j];
227   //    }
228   //     return matrix;
229   //   }
230  
231  
232   //========================================================================================
233   inline itk::Matrix<double, 5, 5> ReadMatrix4D(std::string fileName)
234   {
235     // read input matrix
236     std::ifstream is;
237     openFileForReading(is, fileName);
238     std::vector<double> nb;
239     double x;
240     skipComment(is);
241     is >> x;
242     while (!is.eof()) {
243       nb.push_back(x);
244       skipComment(is);
245       is >> x;  
246     }
247     
248     //copy it to the matrix
249     itk::Matrix<double, 5, 5> matrix;
250     unsigned int index=0;
251     for (unsigned int i=0;i<5;i++)
252       for (unsigned int j=0;j<5;j++)
253         matrix[i][j]=nb[index++];
254     return matrix; 
255   }
256    
257   inline itk::Matrix<double, 4, 4> ReadMatrix3D(std::string fileName)
258   {
259     // read input matrix
260     std::ifstream is;
261     openFileForReading(is, fileName);
262     std::vector<double> nb;
263     double x;
264     skipComment(is);
265     is >> x;
266     while (!is.eof()) {
267       nb.push_back(x);
268       skipComment(is);
269       is >> x;  
270     }
271      
272     //copy it to the matrix
273     itk::Matrix<double, 4, 4> matrix;
274     unsigned int index=0;
275     for (unsigned int i=0;i<4;i++)
276       for (unsigned int j=0;j<4;j++)
277         matrix[i][j]=nb[index++];
278     return matrix; 
279   }
280  
281   inline itk::Matrix<double, 3, 3> ReadMatrix2D(std::string fileName)
282   {
283     // read input matrix
284     std::ifstream is;
285     openFileForReading(is, fileName);
286     std::vector<double> nb;
287     double x;
288     skipComment(is);
289     is >> x;
290     while (!is.eof()) {
291       nb.push_back(x);
292       skipComment(is);
293       is >> x;  
294     }
295      
296     //copy it to the matrix
297     itk::Matrix<double, 3, 3> matrix;
298     unsigned int index=0;
299     for (unsigned int i=0;i<3;i++)
300       for (unsigned int j=0;j<3;j++)
301         matrix[i][j]=nb[index++];
302     return matrix; 
303   }
304    
305   template <unsigned int Dimension > inline itk::Matrix<double, Dimension+1 , Dimension+1> ReadMatrix(std::string fileName)
306   {
307
308     // read input matrix
309     std::ifstream is;
310     openFileForReading(is, fileName);
311     std::vector<double> nb;
312     double x;
313     skipComment(is);
314     is >> x;
315     while (!is.eof()) {
316       nb.push_back(x);
317       skipComment(is);
318       is >> x;  
319     }
320      
321     //copy it to the matrix
322     itk::Matrix<double, Dimension+1, Dimension+1> matrix;
323     unsigned int index=0;
324     for (unsigned int i=0;i<Dimension+1;i++)
325       for (unsigned int j=0;j<Dimension+1;j++)
326         matrix[i][j]=nb[index++];
327     return matrix; 
328   }
329    
330    
331   //    template<> inline itk::Matrix<double, 3, 3>  ReadMatrix<2> (std::string fileName)
332   //    {
333   //      return  ReadMatrix2D(fileName);
334   //    }
335   //    template<> inline itk::Matrix<double, 4, 4>  ReadMatrix<3> (std::string fileName)
336   //    {
337   //      return  ReadMatrix3D(fileName);
338   //    }
339   //  template<> inline itk::Matrix<double, 5, 5>  ReadMatrix<4> (std::string fileName)
340   //    {
341   //      return  ReadMatrix4D(fileName);
342   //    }
343
344
345   //========================================================================================
346   inline  itk::Matrix<double, 4, 4> GetRotationalPartMatrix4D(itk::Matrix<double, 5, 5> input)
347   {
348     itk::Matrix<double,4,4> matrix;
349     matrix[0][0]= input[0][0];
350     matrix[0][1]= input[0][1];
351     matrix[0][2]= input[0][2];
352     matrix[0][3]= input[0][3];
353     matrix[1][0]= input[1][0];
354     matrix[1][1]= input[1][1];
355     matrix[1][2]= input[1][2];
356     matrix[1][3]= input[1][3];
357     matrix[2][0]= input[2][0];
358     matrix[2][1]= input[2][1];
359     matrix[2][2]= input[2][2];
360     matrix[2][2]= input[2][2];
361     matrix[2][3]= input[2][3];
362     matrix[3][0]= input[3][0];
363     matrix[3][1]= input[3][1];
364     matrix[3][2]= input[3][2];
365     matrix[3][2]= input[3][2];
366     matrix[3][3]= input[3][3];
367  
368     return matrix;
369   }
370  
371  
372   inline  itk::Matrix<double, 3, 3> GetRotationalPartMatrix3D(itk::Matrix<double, 4, 4> input)
373   {
374     itk::Matrix<double,3,3> 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];
380     matrix[1][2]= input[1][2];
381     matrix[2][0]= input[2][0];
382     matrix[2][1]= input[2][1];
383     matrix[2][2]= input[2][2];
384     return matrix;
385   }
386    
387   inline  itk::Matrix<double, 2, 2> GetRotationalPartMatrix2D(itk::Matrix<double, 3, 3> input)
388   {
389     itk::Matrix<double,2,2> matrix;
390     matrix[0][0]= input[0][0];
391     matrix[0][1]= input[0][1];
392     matrix[0][2]= input[0][2];
393     matrix[1][0]= input[1][0];
394     matrix[1][1]= input[1][1];
395     return matrix;
396   }
397    
398   inline itk::Matrix<double, 4, 4> GetRotationalPartMatrix(itk::Matrix<double, 5, 5> input)
399   {
400     return GetRotationalPartMatrix4D(input);
401   }
402  
403   inline itk::Matrix<double, 3, 3> GetRotationalPartMatrix(itk::Matrix<double, 4, 4> input)
404   {
405     return GetRotationalPartMatrix3D(input);
406   }
407  
408   inline itk::Matrix<double, 2, 2> GetRotationalPartMatrix(itk::Matrix<double, 3, 3> input)
409   {
410     return GetRotationalPartMatrix2D(input);
411   }
412  
413  
414   //========================================================================================
415   inline  itk::Vector<double,4> GetTranslationPartMatrix4D(itk::Matrix<double, 5, 5> input)
416   {
417     itk::Vector<double,4> vec;
418     vec[0]= input[0][4];
419     vec[1]= input[1][4];
420     vec[2]= input[2][4];
421     vec[3]= input[3][4];
422     return vec;
423  
424   }
425  
426   inline  itk::Vector<double,3> GetTranslationPartMatrix3D(itk::Matrix<double, 4, 4> input)
427   {
428     itk::Vector<double,3> vec;
429     vec[0]= input[0][3];
430     vec[1]= input[1][3];
431     vec[2]= input[2][3];
432     return vec;
433  
434   }
435   inline  itk::Vector<double,2> GetTranslationPartMatrix2D(itk::Matrix<double, 3, 3> input)
436   {
437     itk::Vector<double,2> vec;
438     vec[0]= input[0][2];
439     vec[1]= input[1][2];
440     return vec;
441  
442   }
443  
444   inline itk::Vector<double,4> GetTranslationPartMatrix(itk::Matrix<double, 5, 5> input)
445   {
446  
447     return GetTranslationPartMatrix4D(input);
448   }
449    
450   inline itk::Vector<double,3> GetTranslationPartMatrix(itk::Matrix<double, 4, 4> input)
451   {
452  
453     return GetTranslationPartMatrix3D(input);
454   }
455    
456   inline itk::Vector<double,2> GetTranslationPartMatrix(itk::Matrix<double, 3, 3> input)
457   {
458  
459     return GetTranslationPartMatrix2D(input);
460   }
461 }
462  
463 #endif //#define CLITKTRANSFORMUTILITIES_H