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