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