Cleanup: remove redundant doxygen \file argument
[blender.git] / intern / elbeem / intern / ntl_matrices.h
1 /** \file \ingroup elbeem
2  */
3
4 /******************************************************************************
5  *
6  * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method
7  * Copyright 2003-2006 Nils Thuerey
8  *
9  * Basic matrix utility include file
10  *
11  *****************************************************************************/
12 #ifndef NTL_MATRICES_H
13
14 #include "ntl_vector3dim.h"
15
16 #ifdef WITH_CXX_GUARDEDALLOC
17 #  include "MEM_guardedalloc.h"
18 #endif
19
20 // The basic vector class
21 template<class Scalar>
22 class ntlMatrix4x4
23 {
24 public:
25   // Constructor
26   inline ntlMatrix4x4(void );
27   // Copy-Constructor
28   inline ntlMatrix4x4(const ntlMatrix4x4<Scalar> &v );
29   // construct a vector from one Scalar
30   inline ntlMatrix4x4(Scalar);
31   // construct a vector from three Scalars
32   inline ntlMatrix4x4(Scalar, Scalar, Scalar);
33
34   // Assignment operator
35   inline const ntlMatrix4x4<Scalar>& operator=  (const ntlMatrix4x4<Scalar>& v);
36   // Assignment operator
37   inline const ntlMatrix4x4<Scalar>& operator=  (Scalar s);
38   // Assign and add operator
39   inline const ntlMatrix4x4<Scalar>& operator+= (const ntlMatrix4x4<Scalar>& v);
40   // Assign and add operator
41   inline const ntlMatrix4x4<Scalar>& operator+= (Scalar s);
42   // Assign and sub operator
43   inline const ntlMatrix4x4<Scalar>& operator-= (const ntlMatrix4x4<Scalar>& v);
44   // Assign and sub operator
45   inline const ntlMatrix4x4<Scalar>& operator-= (Scalar s);
46   // Assign and mult operator
47   inline const ntlMatrix4x4<Scalar>& operator*= (const ntlMatrix4x4<Scalar>& v);
48   // Assign and mult operator
49   inline const ntlMatrix4x4<Scalar>& operator*= (Scalar s);
50   // Assign and div operator
51   inline const ntlMatrix4x4<Scalar>& operator/= (const ntlMatrix4x4<Scalar>& v);
52   // Assign and div operator
53   inline const ntlMatrix4x4<Scalar>& operator/= (Scalar s);
54
55   
56   // unary operator
57   inline ntlMatrix4x4<Scalar> operator- () const;
58
59   // binary operator add
60   inline ntlMatrix4x4<Scalar> operator+ (const ntlMatrix4x4<Scalar>&) const;
61   // binary operator add
62   inline ntlMatrix4x4<Scalar> operator+ (Scalar) const;
63   // binary operator sub
64   inline ntlMatrix4x4<Scalar> operator- (const ntlMatrix4x4<Scalar>&) const;
65   // binary operator sub
66   inline ntlMatrix4x4<Scalar> operator- (Scalar) const;
67   // binary operator mult
68   inline ntlMatrix4x4<Scalar> operator* (const ntlMatrix4x4<Scalar>&) const;
69   // binary operator mult
70   inline ntlVector3Dim<Scalar> operator* (const ntlVector3Dim<Scalar>&) const;
71   // binary operator mult
72   inline ntlMatrix4x4<Scalar> operator* (Scalar) const;
73   // binary operator div
74   inline ntlMatrix4x4<Scalar> operator/ (Scalar) const;
75
76         // init function
77         //! init identity matrix
78         inline void initId();
79         //! init rotation matrix
80         inline void initTranslation(Scalar x, Scalar y, Scalar z);
81         //! init rotation matrix
82         inline void initRotationX(Scalar rot);
83         inline void initRotationY(Scalar rot);
84         inline void initRotationZ(Scalar rot);
85         inline void initRotationXYZ(Scalar rotx,Scalar roty, Scalar rotz);
86         //! init scaling matrix
87         inline void initScaling(Scalar scale);
88         inline void initScaling(Scalar x, Scalar y, Scalar z);
89         //! from 16 value array (init id if all 0)
90         inline void initArrayCheck(Scalar *array);
91
92         //! decompose matrix
93         void decompose(ntlVector3Dim<Scalar> &trans,ntlVector3Dim<Scalar> &scale,ntlVector3Dim<Scalar> &rot,ntlVector3Dim<Scalar> &shear);
94
95         //! public to avoid [][] operators
96   Scalar value[4][4];  //< Storage of vector values
97
98
99 protected:
100
101 private:
102 #ifdef WITH_CXX_GUARDEDALLOC
103         MEM_CXX_CLASS_ALLOC_FUNCS("ELBEEM:ntlMatrix4x4")
104 #endif
105 };
106
107
108
109 //------------------------------------------------------------------------------
110 // TYPEDEFS
111 //------------------------------------------------------------------------------
112
113 // a 3D vector for graphics output, typically float?
114 //typedef ntlMatrix4x4<float>  ntlVec3Gfx; 
115
116 //typedef ntlMatrix4x4<double>  ntlMat4d; 
117 typedef ntlMatrix4x4<double>  ntlMat4d; 
118
119 // a 3D vector with single precision
120 typedef ntlMatrix4x4<float>   ntlMat4f; 
121
122 // a 3D vector with grafix precision
123 typedef ntlMatrix4x4<gfxReal>   ntlMat4Gfx;
124
125 // a 3D integer vector
126 typedef ntlMatrix4x4<int>     ntlMat4i; 
127
128
129
130
131
132 //------------------------------------------------------------------------------
133 // STREAM FUNCTIONS
134 //------------------------------------------------------------------------------
135
136
137
138 /*************************************************************************
139   Outputs the object in human readable form using the format
140   [x,y,z]
141   */
142 template<class Scalar>
143 std::ostream&
144 operator<<( std::ostream& os, const ntlMatrix4x4<Scalar>& m )
145 {
146         for(int i=0; i<4; i++) {
147         os << '|' << m.value[i][0] << ", " << m.value[i][1] << ", " << m.value[i][2] << ", " << m.value[i][3] << '|';
148         }
149   return os;
150 }
151
152
153
154 /*************************************************************************
155   Reads the contents of the object from a stream using the same format
156   as the output operator.
157   */
158 template<class Scalar>
159 std::istream&
160 operator>>( std::istream& is, ntlMatrix4x4<Scalar>& m )
161 {
162   char c;
163   char dummy[3];
164   
165         for(int i=0; i<4; i++) {
166         is >> c >> m.value[i][0] >> dummy >> m.value[i][1] >> dummy >> m.value[i][2] >> dummy >> m.value[i][3] >> c;
167         }
168   return is;
169 }
170
171
172 //------------------------------------------------------------------------------
173 // VECTOR inline FUNCTIONS
174 //------------------------------------------------------------------------------
175
176
177
178 /*************************************************************************
179   Constructor.
180   */
181 template<class Scalar>
182 inline ntlMatrix4x4<Scalar>::ntlMatrix4x4( void )
183 {
184 #ifdef MATRIX_INIT_ZERO
185         for(int i=0; i<4; i++) {
186                 for(int j=0; j<4; j++) {
187                 value[i][j] = 0.0;
188                 }
189         }
190 #endif
191 }
192
193
194
195 /*************************************************************************
196   Copy-Constructor.
197   */
198 template<class Scalar>
199 inline ntlMatrix4x4<Scalar>::ntlMatrix4x4( const ntlMatrix4x4<Scalar> &v )
200 {
201   value[0][0] = v.value[0][0]; value[0][1] = v.value[0][1]; value[0][2] = v.value[0][2]; value[0][3] = v.value[0][3];
202   value[1][0] = v.value[1][0]; value[1][1] = v.value[1][1]; value[1][2] = v.value[1][2]; value[1][3] = v.value[1][3];
203   value[2][0] = v.value[2][0]; value[2][1] = v.value[2][1]; value[2][2] = v.value[2][2]; value[2][3] = v.value[2][3];
204   value[3][0] = v.value[3][0]; value[3][1] = v.value[3][1]; value[3][2] = v.value[3][2]; value[3][3] = v.value[3][3];
205 }
206
207
208
209 /*************************************************************************
210   Constructor for a vector from a single Scalar. All components of
211   the vector get the same value.
212   \param s The value to set
213   \return The new vector
214   */
215 template<class Scalar>
216 inline ntlMatrix4x4<Scalar>::ntlMatrix4x4(Scalar s )
217 {
218         for(int i=0; i<4; i++) {
219                 for(int j=0; j<4; j++) {
220                 value[i][j] = s;
221                 }
222         }
223 }
224
225
226
227 /*************************************************************************
228   Copy a ntlMatrix4x4 componentwise.
229   \param v vector with values to be copied
230   \return Reference to self
231   */
232 template<class Scalar>
233 inline const ntlMatrix4x4<Scalar>&
234 ntlMatrix4x4<Scalar>::operator=( const ntlMatrix4x4<Scalar> &v )
235 {
236   value[0][0] = v.value[0][0]; value[0][1] = v.value[0][1]; value[0][2] = v.value[0][2]; value[0][3] = v.value[0][3];
237   value[1][0] = v.value[1][0]; value[1][1] = v.value[1][1]; value[1][2] = v.value[1][2]; value[1][3] = v.value[1][3];
238   value[2][0] = v.value[2][0]; value[2][1] = v.value[2][1]; value[2][2] = v.value[2][2]; value[2][3] = v.value[2][3];
239   value[3][0] = v.value[3][0]; value[3][1] = v.value[3][1]; value[3][2] = v.value[3][2]; value[3][3] = v.value[3][3];
240   return *this;
241 }
242
243
244 /*************************************************************************
245   Copy a Scalar to each component.
246   \param s The value to copy
247   \return Reference to self
248   */
249 template<class Scalar>
250 inline const ntlMatrix4x4<Scalar>&
251 ntlMatrix4x4<Scalar>::operator=(Scalar s)
252 {
253         for(int i=0; i<4; i++) {
254                 for(int j=0; j<4; j++) {
255                 value[i][j] = s;
256                 }
257         }
258   return *this;
259 }
260
261
262 /*************************************************************************
263   Add another ntlMatrix4x4 componentwise.
264   \param v vector with values to be added
265   \return Reference to self
266   */
267 template<class Scalar>
268 inline const ntlMatrix4x4<Scalar>&
269 ntlMatrix4x4<Scalar>::operator+=( const ntlMatrix4x4<Scalar> &v )
270 {
271   value[0][0] += v.value[0][0]; value[0][1] += v.value[0][1]; value[0][2] += v.value[0][2]; value[0][3] += v.value[0][3];
272   value[1][0] += v.value[1][0]; value[1][1] += v.value[1][1]; value[1][2] += v.value[1][2]; value[1][3] += v.value[1][3];
273   value[2][0] += v.value[2][0]; value[2][1] += v.value[2][1]; value[2][2] += v.value[2][2]; value[2][3] += v.value[2][3];
274   value[3][0] += v.value[3][0]; value[3][1] += v.value[3][1]; value[3][2] += v.value[3][2]; value[3][3] += v.value[3][3];
275   return *this;
276 }
277
278
279 /*************************************************************************
280   Add a Scalar value to each component.
281   \param s Value to add
282   \return Reference to self
283   */
284 template<class Scalar>
285 inline const ntlMatrix4x4<Scalar>&
286 ntlMatrix4x4<Scalar>::operator+=(Scalar s)
287 {
288         for(int i=0; i<4; i++) {
289                 for(int j=0; j<4; j++) {
290                 value[i][j] += s;
291                 }
292         }
293   return *this;
294 }
295
296
297 /*************************************************************************
298   Subtract another vector componentwise.
299   \param v vector of values to subtract
300   \return Reference to self
301   */
302 template<class Scalar>
303 inline const ntlMatrix4x4<Scalar>&
304 ntlMatrix4x4<Scalar>::operator-=( const ntlMatrix4x4<Scalar> &v )
305 {
306   value[0][0] -= v.value[0][0]; value[0][1] -= v.value[0][1]; value[0][2] -= v.value[0][2]; value[0][3] -= v.value[0][3];
307   value[1][0] -= v.value[1][0]; value[1][1] -= v.value[1][1]; value[1][2] -= v.value[1][2]; value[1][3] -= v.value[1][3];
308   value[2][0] -= v.value[2][0]; value[2][1] -= v.value[2][1]; value[2][2] -= v.value[2][2]; value[2][3] -= v.value[2][3];
309   value[3][0] -= v.value[3][0]; value[3][1] -= v.value[3][1]; value[3][2] -= v.value[3][2]; value[3][3] -= v.value[3][3];
310   return *this;
311 }
312
313
314 /*************************************************************************
315   Subtract a Scalar value from each component.
316   \param s Value to subtract
317   \return Reference to self
318   */
319 template<class Scalar>
320 inline const ntlMatrix4x4<Scalar>&
321 ntlMatrix4x4<Scalar>::operator-=(Scalar s)
322 {
323         for(int i=0; i<4; i++) {
324                 for(int j=0; j<4; j++) {
325                 value[i][j] -= s;
326                 }
327         }
328   return *this;
329 }
330
331
332 /*************************************************************************
333   Multiply with another vector componentwise.
334   \param v vector of values to multiply with
335   \return Reference to self
336   */
337 template<class Scalar>
338 inline const ntlMatrix4x4<Scalar>&
339 ntlMatrix4x4<Scalar>::operator*=( const ntlMatrix4x4<Scalar> &v )
340 {
341         ntlMatrix4x4<Scalar> nv(0.0);
342         for(int i=0; i<4; i++) {
343                 for(int j=0; j<4; j++) {
344
345                         for(int k=0;k<4;k++)
346                                 nv.value[i][j] += (value[i][k] * v.value[k][j]);
347                 }
348         }
349   *this = nv;
350   return *this;
351 }
352
353
354 /*************************************************************************
355   Multiply each component with a Scalar value.
356   \param s Value to multiply with
357   \return Reference to self
358   */
359 template<class Scalar>
360 inline const ntlMatrix4x4<Scalar>&
361 ntlMatrix4x4<Scalar>::operator*=(Scalar s)
362 {
363         for(int i=0; i<4; i++) {
364                 for(int j=0; j<4; j++) {
365                 value[i][j] *= s;
366                 }
367         }
368   return *this;
369 }
370
371
372
373 /*************************************************************************
374   Divide each component by a Scalar value.
375   \param s Value to divide by
376   \return Reference to self
377   */
378 template<class Scalar>
379 inline const ntlMatrix4x4<Scalar>&
380 ntlMatrix4x4<Scalar>::operator/=(Scalar s)
381 {
382         for(int i=0; i<4; i++) {
383                 for(int j=0; j<4; j++) {
384                 value[i][j] /= s;
385                 }
386         }
387   return *this;
388 }
389
390
391 //------------------------------------------------------------------------------
392 // unary operators
393 //------------------------------------------------------------------------------
394
395
396 /*************************************************************************
397   Build componentwise the negative this vector.
398   \return The new (negative) vector
399   */
400 template<class Scalar>
401 inline ntlMatrix4x4<Scalar>
402 ntlMatrix4x4<Scalar>::operator-() const
403 {
404         ntlMatrix4x4<Scalar> nv;
405         for(int i=0; i<4; i++) {
406                 for(int j=0; j<4; j++) {
407                 nv[i][j] = -value[i][j];
408                 }
409         }
410   return nv;
411 }
412
413
414
415 //------------------------------------------------------------------------------
416 // binary operators
417 //------------------------------------------------------------------------------
418
419
420 /*************************************************************************
421   Build a vector with another vector added componentwise.
422   \param v The second vector to add
423   \return The sum vector
424   */
425 template<class Scalar>
426 inline ntlMatrix4x4<Scalar>
427 ntlMatrix4x4<Scalar>::operator+( const ntlMatrix4x4<Scalar> &v ) const
428 {
429         ntlMatrix4x4<Scalar> nv;
430         for(int i=0; i<4; i++) {
431                 for(int j=0; j<4; j++) {
432                 nv[i][j] = value[i][j] + v.value[i][j];
433                 }
434         }
435   return nv;
436 }
437
438
439 /*************************************************************************
440   Build a vector with a Scalar value added to each component.
441   \param s The Scalar value to add
442   \return The sum vector
443   */
444 template<class Scalar>
445 inline ntlMatrix4x4<Scalar>
446 ntlMatrix4x4<Scalar>::operator+(Scalar s) const
447 {
448         ntlMatrix4x4<Scalar> nv;
449         for(int i=0; i<4; i++) {
450                 for(int j=0; j<4; j++) {
451                 nv[i][j] = value[i][j] + s;
452                 }
453         }
454   return nv;
455 }
456
457
458 /*************************************************************************
459   Build a vector with another vector subtracted componentwise.
460   \param v The second vector to subtract
461   \return The difference vector
462   */
463 template<class Scalar>
464 inline ntlMatrix4x4<Scalar>
465 ntlMatrix4x4<Scalar>::operator-( const ntlMatrix4x4<Scalar> &v ) const
466 {
467         ntlMatrix4x4<Scalar> nv;
468         for(int i=0; i<4; i++) {
469                 for(int j=0; j<4; j++) {
470                 nv[i][j] = value[i][j] - v.value[i][j];
471                 }
472         }
473   return nv;
474 }
475
476
477 /*************************************************************************
478   Build a vector with a Scalar value subtracted componentwise.
479   \param s The Scalar value to subtract
480   \return The difference vector
481   */
482 template<class Scalar>
483 inline ntlMatrix4x4<Scalar>
484 ntlMatrix4x4<Scalar>::operator-(Scalar s ) const
485 {
486         ntlMatrix4x4<Scalar> nv;
487         for(int i=0; i<4; i++) {
488                 for(int j=0; j<4; j++) {
489                 nv[i][j] = value[i][j] - s;
490                 }
491         }
492   return nv;
493 }
494
495
496
497 /*************************************************************************
498   Build a ntlMatrix4x4 with a Scalar value multiplied to each component.
499   \param s The Scalar value to multiply with
500   \return The product vector
501   */
502 template<class Scalar>
503 inline ntlMatrix4x4<Scalar>
504 ntlMatrix4x4<Scalar>::operator*(Scalar s) const
505 {
506         ntlMatrix4x4<Scalar> nv;
507         for(int i=0; i<4; i++) {
508                 for(int j=0; j<4; j++) {
509                 nv[i][j] = value[i][j] * s;
510                 }
511         }
512   return nv;
513 }
514
515
516
517
518 /*************************************************************************
519   Build a vector divided componentwise by a Scalar value.
520   \param s The Scalar value to divide by
521   \return The ratio vector
522   */
523 template<class Scalar>
524 inline ntlMatrix4x4<Scalar>
525 ntlMatrix4x4<Scalar>::operator/(Scalar s) const
526 {
527         ntlMatrix4x4<Scalar> nv;
528         for(int i=0; i<4; i++) {
529                 for(int j=0; j<4; j++) {
530                 nv[i][j] = value[i][j] / s;
531                 }
532         }
533   return nv;
534 }
535
536
537
538
539
540 /*************************************************************************
541   Build a vector with another vector multiplied by componentwise.
542   \param v The second vector to muliply with
543   \return The product vector
544   */
545 template<class Scalar>
546 inline ntlMatrix4x4<Scalar>
547 ntlMatrix4x4<Scalar>::operator*( const ntlMatrix4x4<Scalar>& v) const
548 {
549         ntlMatrix4x4<Scalar> nv(0.0);
550         for(int i=0; i<4; i++) {
551                 for(int j=0; j<4; j++) {
552
553                         for(int k=0;k<4;k++)
554                                 nv.value[i][j] += (value[i][k] * v.value[k][j]);
555                 }
556         }
557   return nv;
558 }
559
560
561 template<class Scalar>
562 inline ntlVector3Dim<Scalar>
563 ntlMatrix4x4<Scalar>::operator*( const ntlVector3Dim<Scalar>& v) const
564 {
565         ntlVector3Dim<Scalar> nvec(0.0);
566         for(int i=0; i<3; i++) {
567                 for(int j=0; j<3; j++) {
568                         nvec[i] += (v[j] * value[i][j]);
569                 }
570         }
571         // assume normalized w coord
572         for(int i=0; i<3; i++) {
573                 nvec[i] += (1.0 * value[i][3]);
574         }
575   return nvec;
576 }
577
578
579
580 //------------------------------------------------------------------------------
581 // Other helper functions
582 //------------------------------------------------------------------------------
583
584 //! init identity matrix
585 template<class Scalar>
586 inline void ntlMatrix4x4<Scalar>::initId()
587 {
588         (*this) = (Scalar)(0.0);
589         value[0][0] = 
590         value[1][1] = 
591         value[2][2] = 
592         value[3][3] = (Scalar)(1.0);
593 }
594
595 //! init rotation matrix
596 template<class Scalar>
597 inline void ntlMatrix4x4<Scalar>::initTranslation(Scalar x, Scalar y, Scalar z)
598 {
599         //(*this) = (Scalar)(0.0);
600         this->initId();
601         value[0][3] = x;
602         value[1][3] = y;
603         value[2][3] = z;
604 }
605
606 //! init rotation matrix
607 template<class Scalar>
608 inline void 
609 ntlMatrix4x4<Scalar>::initRotationX(Scalar rot)
610 {
611         double drot = (double)(rot/360.0*2.0*M_PI);
612         //? while(drot < 0.0) drot += (M_PI*2.0);
613
614         this->initId();
615         value[1][1] = (Scalar)  cos(drot);
616         value[1][2] = (Scalar)  sin(drot);
617         value[2][1] = (Scalar)(-sin(drot));
618         value[2][2] = (Scalar)  cos(drot);
619 }
620 template<class Scalar>
621 inline void 
622 ntlMatrix4x4<Scalar>::initRotationY(Scalar rot)
623 {
624         double drot = (double)(rot/360.0*2.0*M_PI);
625         //? while(drot < 0.0) drot += (M_PI*2.0);
626
627         this->initId();
628         value[0][0] = (Scalar)  cos(drot);
629         value[0][2] = (Scalar)(-sin(drot));
630         value[2][0] = (Scalar)  sin(drot);
631         value[2][2] = (Scalar)  cos(drot);
632 }
633 template<class Scalar>
634 inline void 
635 ntlMatrix4x4<Scalar>::initRotationZ(Scalar rot)
636 {
637         double drot = (double)(rot/360.0*2.0*M_PI);
638         //? while(drot < 0.0) drot += (M_PI*2.0);
639
640         this->initId();
641         value[0][0] = (Scalar)  cos(drot);
642         value[0][1] = (Scalar)  sin(drot);
643         value[1][0] = (Scalar)(-sin(drot));
644         value[1][1] = (Scalar)  cos(drot);
645 }
646 template<class Scalar>
647 inline void 
648 ntlMatrix4x4<Scalar>::initRotationXYZ( Scalar rotx, Scalar roty, Scalar rotz)
649 {
650         ntlMatrix4x4<Scalar> val;
651         ntlMatrix4x4<Scalar> rot;
652         this->initId();
653
654         // org
655         /*rot.initRotationX(rotx);
656         (*this) *= rot;
657         rot.initRotationY(roty);
658         (*this) *= rot;
659         rot.initRotationZ(rotz);
660         (*this) *= rot;
661         // org */
662
663         // blender
664         rot.initRotationZ(rotz);
665         (*this) *= rot;
666         rot.initRotationY(roty);
667         (*this) *= rot;
668         rot.initRotationX(rotx);
669         (*this) *= rot;
670         // blender */
671 }
672
673 //! init scaling matrix
674 template<class Scalar>
675 inline void 
676 ntlMatrix4x4<Scalar>::initScaling(Scalar scale)
677 {
678         this->initId();
679         value[0][0] = scale;
680         value[1][1] = scale;
681         value[2][2] = scale;
682 }
683 //! init scaling matrix
684 template<class Scalar>
685 inline void 
686 ntlMatrix4x4<Scalar>::initScaling(Scalar x, Scalar y, Scalar z)
687 {
688         this->initId();
689         value[0][0] = x;
690         value[1][1] = y;
691         value[2][2] = z;
692 }
693
694
695 //! from 16 value array (init id if all 0)
696 template<class Scalar>
697 inline void 
698 ntlMatrix4x4<Scalar>::initArrayCheck(Scalar *array)
699 {
700         bool allZero = true;
701         for(int i=0; i<4; i++) {
702                 for(int j=0; j<4; j++) {
703                         value[i][j] = array[i*4+j];
704                         if(array[i*4+j]!=0.0) allZero=false;
705                 }
706         }
707         if(allZero) this->initId();
708 }
709
710 //! decompose matrix
711 template<class Scalar>
712 void 
713 ntlMatrix4x4<Scalar>::decompose(ntlVector3Dim<Scalar> &trans,ntlVector3Dim<Scalar> &scale,ntlVector3Dim<Scalar> &rot,ntlVector3Dim<Scalar> &shear) {
714         ntlVec3Gfx row[3],temp;
715
716         for(int i = 0; i < 3; i++) {
717                 trans[i] = this->value[3][i];
718         }
719
720         for(int i = 0; i < 3; i++) {
721                 row[i][0] = this->value[i][0];
722                 row[i][1] = this->value[i][1];
723                 row[i][2] = this->value[i][2];
724         }
725
726         scale[0] = norm(row[0]);
727         normalize (row[0]);
728
729         shear[0] = dot(row[0], row[1]);
730         row[1][0] = row[1][0] - shear[0]*row[0][0];
731         row[1][1] = row[1][1] - shear[0]*row[0][1];
732         row[1][2] = row[1][2] - shear[0]*row[0][2];
733
734         scale[1] = norm(row[1]);
735         normalize (row[1]);
736
737         if(scale[1] != 0.0)
738                 shear[0] /= scale[1];
739
740         shear[1] = dot(row[0], row[2]);
741         row[2][0] = row[2][0] - shear[1]*row[0][0];
742         row[2][1] = row[2][1] - shear[1]*row[0][1];
743         row[2][2] = row[2][2] - shear[1]*row[0][2];
744
745         shear[2] = dot(row[1], row[2]);
746         row[2][0] = row[2][0] - shear[2]*row[1][0];
747         row[2][1] = row[2][1] - shear[2]*row[1][1];
748         row[2][2] = row[2][2] - shear[2]*row[1][2];
749
750         scale[2] = norm(row[2]);
751         normalize (row[2]);
752
753         if(scale[2] != 0.0) {
754                 shear[1] /= scale[2];
755                 shear[2] /= scale[2];
756         }
757
758         temp = cross(row[1], row[2]);
759         if(dot(row[0], temp) < 0.0) {
760                 for(int i = 0; i < 3; i++) {
761                         scale[i]  *= -1.0;
762                         row[i][0] *= -1.0;
763                         row[i][1] *= -1.0;
764                         row[i][2] *= -1.0;
765                 }
766         }
767
768         if(row[0][2] < -1.0) row[0][2] = -1.0;
769         if(row[0][2] > +1.0) row[0][2] = +1.0;
770
771         rot[1] = asin(-row[0][2]);
772
773         if(fabs(cos(rot[1])) > VECTOR_EPSILON) {
774                 rot[0] = atan2 (row[1][2], row[2][2]);
775                 rot[2] = atan2 (row[0][1], row[0][0]);
776         }
777         else {
778                 rot[0] = atan2 (row[1][0], row[1][1]);
779                 rot[2] = 0.0;
780         }
781
782         rot[0] = (180.0/M_PI)*rot[0];
783         rot[1] = (180.0/M_PI)*rot[1];
784         rot[2] = (180.0/M_PI)*rot[2];
785
786
787 #define NTL_MATRICES_H
788 #endif
789