KWStyle - itkVersor.txx
 
Matrix View
Description
1 /*=========================================================================
2
3   Program:   Insight Segmentation & Registration Toolkit
4   Module:    $RCSfile: itkVersor.txx.html,v $
5   Language:  C++
6   Date:      $Date: 2006/01/17 19:15:49 $
7   Version:   $Revision: 1.4 $
8
9   Copyright (c) Insight Software Consortium. All rights reserved.
10   See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
11
12      This software is distributed WITHOUT ANY WARRANTY; without even 
13      the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
14      PURPOSE.  See the above copyright notices for more information.
15
16 =========================================================================*/
17 DEF #ifndef _itkVersor_txx
18 DEF #define _itkVersor_txx
19
20 #include "itkVersor.h"
21 #include "itkVector.h" 
22 #include "itkNumericTraits.h" 
23 #include "itkExceptionObject.h"
24
25
26 namespace itk
27 {
28
29
30 /**
31  * Constructor to initialize entire vector to one value.
32  */
33 template<class T>
34 Versor<T>
35 ::Versor()
36 {
37   m_X = NumericTraits<T>::Zero;
38   m_Y = NumericTraits<T>::Zero;
39   m_Z = NumericTraits<T>::Zero;
40   m_W = NumericTraits<T>::One;
41 }
42
43
44 EML
45 /**
46  * Copy Constructor 
47  */
48 template<class T>
49 Versor<T>
50 ::Versor( const Self & v)
51 {
52   m_X = v.m_X;
53   m_Y = v.m_Y;
54   m_Z = v.m_Z;
55   m_W = v.m_W;
56 }
57
58
59 EML
60 /**
61  * Assignment Operator
62  */
63 template<class T>
64 const Versor<T> &
65 Versor<T>
66 ::operator=( const Self & v)
67 {
68   m_X = v.m_X;
69   m_Y = v.m_Y;
70   m_Z = v.m_Z;
71   m_W = v.m_W;
72   return *this;
73 }
74
75
76 /**
77  * Set to an identity transform
78  */
79 template<class T>
80 void
81 Versor<T>
82 ::SetIdentity()
83 {
84   m_X = NumericTraits<T>::Zero;
85   m_Y = NumericTraits<T>::Zero;
86   m_Z = NumericTraits<T>::Zero;
87   m_W = NumericTraits<T>::One;
88 }
89
90
91 EML
92 EML
93 /**
94  * Return a vnl_quaternion
95  */
96 template<class T>
97 vnl_quaternion<T>
98 Versor<T>
99 ::GetVnlQuaternion(void) const
100 {
101   return vnl_quaternion<T>(m_X,m_Y,m_Z,m_W);
102 }
103
104
105 EML
106 /**
107  * Assignment and Composition Operator
108  */
109 template<class T>
110 const Versor<T> &
111 Versor<T>
112 ::operator*=( const Self & v)
113 {
114
115   const double mx =  m_W*v.m_X - m_Z*v.m_Y + m_Y*v.m_Z + m_X*v.m_W;
116   const double my =  m_Z*v.m_X + m_W*v.m_Y - m_X*v.m_Z + m_Y*v.m_W;
117   const double mz = -m_Y*v.m_X + m_X*v.m_Y + m_W*v.m_Z + m_Z*v.m_W;
118   const double mw = -m_X*v.m_X - m_Y*v.m_Y - m_Z*v.m_Z + m_W*v.m_W;
119
120   m_X = mx;
121   m_Y = my;
122   m_Z = mz;
123   m_W = mw;
124
125   return *this;
126 }
127
128
129 EML
130 /**
131  * Composition Operator
132  */
133 template<class T>
134 Versor<T> 
135 Versor<T>
136 ::operator*( const Self & v) const
137 {
138  
139   Self result;
140
141   result.m_X =  m_W*v.m_X - m_Z*v.m_Y + m_Y*v.m_Z + m_X*v.m_W;
142   result.m_Y =  m_Z*v.m_X + m_W*v.m_Y - m_X*v.m_Z + m_Y*v.m_W;
143   result.m_Z = -m_Y*v.m_X + m_X*v.m_Y + m_W*v.m_Z + m_Z*v.m_W;
144   result.m_W = -m_X*v.m_X - m_Y*v.m_Y - m_Z*v.m_Z + m_W*v.m_W;
145
146   return result;
147 }
148
149
150 EML
151 EML
152 /**
153  * Division and Assignment Operator
154  */
155 template<class T>
156 const Versor<T> &
157 Versor<T>
158 ::operator/=( const Self & v)
159 {
160   
161   const double mx = -m_W*v.m_X + m_Z*v.m_Y - m_Y*v.m_Z + m_X*v.m_W;
162   const double my = -m_Z*v.m_X - m_W*v.m_Y + m_X*v.m_Z + m_Y*v.m_W;
163   const double mz =  m_Y*v.m_X - m_X*v.m_Y - m_W*v.m_Z + m_Z*v.m_W;
164   const double mw =  m_X*v.m_X + m_Y*v.m_Y + m_Z*v.m_Z + m_W*v.m_W;
165
166   m_X = mx;
167   m_Y = my;
168   m_Z = mz;
169   m_W = mw;
170
171   return *this;
172
173 }
174
175
176 /**
177  * Division Operator 
178  */
179 template<class T>
180 Versor<T> 
181 Versor<T>
182 ::operator/( const Self & v) const
183 {
184   
185   Self result;
186
187   result.m_X = -m_W*v.m_X + m_Z*v.m_Y - m_Y*v.m_Z + m_X*v.m_W;
188   result.m_Y = -m_Z*v.m_X - m_W*v.m_Y + m_X*v.m_Z + m_Y*v.m_W;
189   result.m_Z =  m_Y*v.m_X - m_X*v.m_Y - m_W*v.m_Z + m_Z*v.m_W;
190   result.m_W =  m_X*v.m_X + m_Y*v.m_Y + m_Z*v.m_Z + m_W*v.m_W;
191
192   return result;
193 }
194
195
196 EML
197 EML
198 /**
199  * Comparision operator
200  */
201 template<class T>
202 bool
203 Versor<T>
204 ::operator!=( const Self & v) const
205 {
206   return !(*this == v);
207 }
208
209
210 EML
211 /**
212  * Comparision operator
213  */
214 template<class T>
215 bool
216 Versor<T>
217 ::operator==( const Self & v) const
218 {
219
220   // Evaluate the quaternion ratio between them
221   Self ratio = *this * v.GetReciprocal();
222   
223   const typename itk::NumericTraits< T >::AccumulateType 
224 IND ****square = ratio.m_W * ratio.m_W;
225   
226   const double epsilon = 1e-300;
227
228   if( fabs( 1.0f - square ) < epsilon )
229     {
230     return true;
231     }
232
233   return false;  
234
235 }
236
237
238 EML
239 /**
240  * Get Conjugate
241  */
242 template<class T>
243 Versor<T> 
244 Versor<T>
245 ::GetConjugate( void ) const
246 {
247   Self result;
248   
249   result.m_X = -m_X;
250   result.m_Y = -m_Y;
251   result.m_Z = -m_Z;
252   result.m_W =  m_W;
253
254   return result;
255 }
256  
257
258
259 EML
260 /**
261  * Get Reciprocal
262  */
263 template<class T>
264 Versor<T> 
265 Versor<T>
266 ::GetReciprocal( void ) const
267 {
268   Self result;
269   
270   result.m_X = -m_X;
271   result.m_Y = -m_Y;
272   result.m_Z = -m_Z;
273   result.m_W =  m_W;
274
275   return result;
276 }
277  
278
279
280 EML
281 /**
282  * Get Tensor part
283  */
284 template<class T>
285 typename Versor<T>::ValueType
286 Versor<T>
287 ::GetTensor( void ) const
288 {
289   
290   const ValueType tensor = 
291 IND ****static_cast< ValueType > (
292       sqrt( m_X*m_X + m_Y*m_Y + m_Z*m_Z + m_W*m_W ) );
293
294   return tensor;
295 }
296  
297
298
299 /**
300  * Normalize
301  */
302 template<class T>
303 void
304 Versor<T>
305 ::Normalize( void ) 
306 {
307   const ValueType tensor = this->GetTensor();
308
309   if( fabs( tensor ) < 1e-20 )
310     {
311     ExceptionObject except;
312 LEN     except.SetDescription("Attempt to normalize a itk::Versor with zero tensor");
313     except.SetLocation(__FILE__);
314     throw except;
315     }
316   m_X /=  tensor;
317   m_Y /=  tensor;
318   m_Z /=  tensor;
319   m_W /=  tensor;
320
321 }
322  
323
324
325 EML
326 EML
327 /**
328  * Get Axis
329  */
330 template<class T>
331 typename Versor<T>::VectorType 
332 Versor<T>
333 ::GetAxis( void ) const
334 {
335
336   VectorType axis;
337   
338   const RealType ax = static_cast<RealType>( m_X );
339   const RealType ay = static_cast<RealType>( m_Y );
340   const RealType az = static_cast<RealType>( m_Z );
341   
342   const RealType vectorNorm = sqrt( ax * ax  +  ay * ay  +  az * az  );
343
344   if( vectorNorm == NumericTraits<RealType>::Zero )
345     {
346     axis[0] = NumericTraits<T>::Zero;
347     axis[1] = NumericTraits<T>::Zero;
348     axis[2] = NumericTraits<T>::Zero;
349     }
350   else 
351     {
352     axis[0] = m_X / vectorNorm;
353     axis[1] = m_Y / vectorNorm;
354     axis[2] = m_Z / vectorNorm;
355     }
356
357   return axis;
358 }
359  
360
361
362 /**
363  * Get Right part
364  */
365 template<class T>
366 typename Versor<T>::VectorType 
367 Versor<T>
368 ::GetRight( void ) const
369 {
370   VectorType axis;
371   
372   axis[0] = m_X;
373   axis[1] = m_Y;
374   axis[2] = m_Z;
375
376   return axis;
377 }
378  
379
380
381 EML
382 /**
383  * Get Scalar part
384  */
385 template<class T>
386 typename Versor<T>::ValueType 
387 Versor<T>
388 ::GetScalar( void ) const
389 {
390   return m_W;
391 }
392  
393
394
395 /**
396  * Get Angle (in radians)
397  */
398 template<class T>
399 typename Versor<T>::ValueType
400 Versor<T>
401 ::GetAngle( void ) const
402 {
403   const RealType ax = static_cast<RealType>( m_X );
404   const RealType ay = static_cast<RealType>( m_Y );
405   const RealType az = static_cast<RealType>( m_Z );
406   
407   const RealType vectorNorm = sqrt( ax * ax  +  ay * ay  +  az * az  );
408
409   const ValueType angle = 2.0 * atan2( vectorNorm, m_W );
410   
411   return angle;
412
413 }
414  
415
416
417 /**
418  * Get the Square root of the unit quaternion
419  */
420 template<class T>
421 Versor<T>
422 Versor<T>
423 ::SquareRoot( void ) const
424 {
425
426   const ValueType newScalar = sqrt( static_cast<double>( 1.0 + m_W ) );
427   const double sqrtOfTwo    = sqrt( 2.0f );
428
429   const double factor = 1.0f / ( newScalar * sqrtOfTwo );
430
431   Self result;
432
433   result.m_X = m_X * factor;
434   result.m_Y = m_Y * factor;
435   result.m_Z = m_Z * factor;
436   result.m_W = newScalar / sqrtOfTwo;
437
438   return result;
439
440 }
441  
442
443 /**
444  *  Compute the Exponential of the quaternion
445  */
446 template<class T>
447 Versor<T>
448 Versor<T>
449 ::Exponential( ValueType exponent ) const
450 {
451
452   Self result;
453   
454   result.Set( this->GetAxis(), 
455 IND **************this->GetAngle() * exponent );
456
457   return result;
458
459 }
460  
461
462
463 /**
464  * Set Axis and Angle (in radians)
465  */
466 template<class T>
467 void
468 Versor<T>
469 ::Set( const VectorType & axis, ValueType angle )
470 {
471
472   const RealType vectorNorm = axis.GetNorm();
473
474   const RealType cosangle2 = cos( angle / 2.0 );
475   const RealType sinangle2 = sin( angle / 2.0 );
476   
477   const RealType factor = sinangle2 / vectorNorm;
478   
479   m_X = axis[0] * factor;
480   m_Y = axis[1] * factor;
481   m_Z = axis[2] * factor;
482   
483   m_W = cosangle2;
484
485 }
486  
487
488 /**
489  * Set using an orthogonal matrix.
490  */
491 template<class T>
492 void
493 Versor<T>
494 ::Set( const MatrixType & mat )
495 {
496   vnl_matrix<double> m( mat.GetVnlMatrix() );
497
498   const double epsilon = 1e-30;
499
500   double trace = m(0,0) + m(1,1) + m(2,2) + 1.0;
501
502   if( trace > epsilon)
503     {
504     const double s = 0.5 / sqrt(trace);
505     m_W = 0.25 / s;
506     m_X = (m(2,1) - m(1,2)) * s;
507     m_Y = (m(0,2) - m(2,0)) * s;
508     m_Z = (m(1,0) - m(0,1)) * s;
509     }
510   else
511     {
512     if( m(0,0) > m(1,1) && m(0,0) > m(2,2) )
513       {
514       const double s = 2.0 * sqrt(1.0 + m(0,0) - m(1,1) - m(2,2));
515       m_X = 0.25 * s;
516       m_Y = (m(0,1) + m(1,0)) / s;
517       m_Z = (m(0,2) + m(2,0)) / s;
518       m_W = (m(1,2) - m(2,1)) / s;
519       }
520     else
521       {
522       if( m(1,1) > m(2,2) )
523         {
524         const double s = 2.0 * sqrt(1.0 + m(1,1) - m(0,0) - m(2,2));
525         m_X = (m(0,1) + m(1,0)) / s;
526         m_Y = 0.25 * s;
527         m_Z = (m(1,2) + m(2,1)) / s;
528         m_W = (m(0,2) - m(2,0)) / s;
529         }
530       else
531         {
532         const double s = 2.0 * sqrt(1.0 + m(2,2) - m(0,0) - m(1,1));
533         m_X = (m(0,2) + m(2,0)) / s;
534         m_Y = (m(1,2) + m(2,1)) / s;
535         m_Z = 0.25 * s;
536         m_W = (m(0,1) - m(1,0)) / s;
537         }
538       }
539     }
540   this->Normalize();
541 }
542
543
544 /**
545  * Set right Part (in radians)
546  */
547 template<class T>
548 void
549 Versor<T>
550 ::Set( const VectorType & axis )
551 {
552
553   const ValueType sinangle2 =  axis.GetNorm();
554   if( sinangle2 > 1.0 )
555     {
556     ExceptionObject exception;
557     exception.SetDescription("Trying to initialize a Versor with " \
558                              "a vector whose magnitude is greater than 1");
559     exception.SetLocation("itk::Versor::Set( const VectorType )");
560     throw exception;
561     }
562   
563   const ValueType cosangle2 =  sqrt( 1.0 - sinangle2 * sinangle2 );
564   
565   m_X = axis[0];
566   m_Y = axis[1];
567   m_Z = axis[2];
568   
569   m_W = cosangle2;
570
571 }
572  
573
574
575 /**
576  * Set the Versor from four components.
577  * After assignment, the quaternion is normalized
578  * in order to get a consistent Versor (unit quaternion).
579  */
580 template<class T>
581 void
582 Versor<T>
583 ::Set( T x, T y, T z, T w )
584 {
585   m_X = x;
586   m_Y = y;
587   m_Z = z;
588   m_W = w;
589   this->Normalize();
590 }
591
592
593 EML
594 /**
595  * Set from a vnl_quaternion
596  * After assignment, the quaternion is normalized
597  * in order to get a consistent Versor (unit quaternion).
598  */
599 template<class T>
600 void
601 Versor<T>
602 ::Set( const VnlQuaternionType & quaternion )
603 {
604   m_X = quaternion.x();
605   m_Y = quaternion.y();
606   m_Z = quaternion.z();
607   m_W = quaternion.r();
608   this->Normalize();
609 }
610
611
612 EML
613 EML
614 /**
615  * Set rotation around X axis 
616  */
617 template<class T>
618 void
619 Versor<T>
620 ::SetRotationAroundX( ValueType angle )
621 {
622
623   const ValueType sinangle2 = sin( angle / 2.0 );
624   const ValueType cosangle2 = cos( angle / 2.0 );
625   
626   m_X = sinangle2;
627   m_Y = NumericTraits< T >::Zero;
628   m_Z = NumericTraits< T >::Zero;
629   m_W = cosangle2;
630
631 }
632  
633
634
635 /**
636  * Set rotation around Y axis 
637  */
638 template<class T>
639 void
640 Versor<T>
641 ::SetRotationAroundY( ValueType angle )
642 {
643
644   const ValueType sinangle2 = sin( angle / 2.0 );
645   const ValueType cosangle2 = cos( angle / 2.0 );
646   
647   m_X = NumericTraits< T >::Zero;
648   m_Y = sinangle2;
649   m_Z = NumericTraits< T >::Zero;
650   m_W = cosangle2;
651
652 }
653  
654
655
656 /**
657  * Set rotation around Z axis 
658  */
659 template<class T>
660 void
661 Versor<T>
662 ::SetRotationAroundZ( ValueType angle )
663 {
664
665   const ValueType sinangle2 = sin( angle / 2.0 );
666   const ValueType cosangle2 = cos( angle / 2.0 );
667   
668   m_X = NumericTraits< T >::Zero;
669   m_Y = NumericTraits< T >::Zero;
670   m_Z = sinangle2;
671   m_W = cosangle2;
672
673 }
674  
675
676
677 /**
678  * Transform a Vector
679  */
680 template<class T>
681 typename Versor<T>::VectorType 
682 Versor<T>
683 ::Transform( const VectorType & v ) const
684 {
685   VectorType result;
686   
687   const ValueType xx = m_X * m_X;
688   const ValueType yy = m_Y * m_Y;
689   const ValueType zz = m_Z * m_Z;
690   const ValueType xy = m_X * m_Y;
691   const ValueType xz = m_X * m_Z;
692   const ValueType xw = m_X * m_W;
693   const ValueType yz = m_Y * m_Z;
694   const ValueType yw = m_Y * m_W;
695   const ValueType zw = m_Z * m_W;
696
697   const ValueType mxx = 1.0 - 2.0 * ( yy + zz );
698   const ValueType myy = 1.0 - 2.0 * ( xx + zz );
699   const ValueType mzz = 1.0 - 2.0 * ( xx + yy );
700   const ValueType mxy = 2.0 * ( xy - zw );
701   const ValueType mxz = 2.0 * ( xz + yw );
702   const ValueType myx = 2.0 * ( xy + zw );
703   const ValueType mzx = 2.0 * ( xz - yw );
704   const ValueType mzy = 2.0 * ( yz + xw );
705   const ValueType myz = 2.0 * ( yz - xw );
706     
707   result[0] = mxx * v[0] + mxy * v[1] + mxz * v[2];
708   result[1] = myx * v[0] + myy * v[1] + myz * v[2];
709   result[2] = mzx * v[0] + mzy * v[1] + mzz * v[2];
710
711   return result;
712 }
713  
714
715
716 EML
717 EML
718 /**
719  * Transform a CovariantVector
720  * given that this is an orthogonal transformation
721  * CovariantVectors are transformed as vectors.
722  */
723 template<class T>
724 typename Versor<T>::CovariantVectorType 
725 Versor<T>
726 ::Transform( const CovariantVectorType & v ) const
727 {
728   CovariantVectorType result;
729   
730   const ValueType xx = m_X * m_X;
731   const ValueType yy = m_Y * m_Y;
732   const ValueType zz = m_Z * m_Z;
733   const ValueType xy = m_X * m_Y;
734   const ValueType xz = m_X * m_Z;
735   const ValueType xw = m_X * m_W;
736   const ValueType yz = m_Y * m_Z;
737   const ValueType yw = m_Y * m_W;
738   const ValueType zw = m_Z * m_W;
739
740   const ValueType mxx = 1.0 - 2.0 * ( yy + zz );
741   const ValueType myy = 1.0 - 2.0 * ( xx + zz );
742   const ValueType mzz = 1.0 - 2.0 * ( xx + yy );
743   const ValueType mxy = 2.0 * ( xy - zw );
744   const ValueType mxz = 2.0 * ( xz + yw );
745   const ValueType myx = 2.0 * ( xy + zw );
746   const ValueType mzx = 2.0 * ( xz - yw );
747   const ValueType mzy = 2.0 * ( yz + xw );
748   const ValueType myz = 2.0 * ( yz - xw );
749  
750   result[0] = mxx * v[0] + mxy * v[1] + mxz * v[2];
751   result[1] = myx * v[0] + myy * v[1] + myz * v[2];
752   result[2] = mzx * v[0] + mzy * v[1] + mzz * v[2];
753
754   return result;
755 }
756
757
758 EML
759 EML
760 /**
761  * Transform a Point
762  */
763 template<class T>
764 typename Versor<T>::PointType 
765 Versor<T>
766 ::Transform( const PointType & v ) const
767 {
768   PointType result;
769   
770   const ValueType xx = m_X * m_X;
771   const ValueType yy = m_Y * m_Y;
772   const ValueType zz = m_Z * m_Z;
773   const ValueType xy = m_X * m_Y;
774   const ValueType xz = m_X * m_Z;
775   const ValueType xw = m_X * m_W;
776   const ValueType yz = m_Y * m_Z;
777   const ValueType yw = m_Y * m_W;
778   const ValueType zw = m_Z * m_W;
779  
780   const ValueType mxx = 1.0 - 2.0 * ( yy + zz );
781   const ValueType myy = 1.0 - 2.0 * ( xx + zz );
782   const ValueType mzz = 1.0 - 2.0 * ( xx + yy );
783   const ValueType mxy = 2.0 * ( xy - zw );
784   const ValueType mxz = 2.0 * ( xz + yw );
785   const ValueType myx = 2.0 * ( xy + zw );
786   const ValueType mzx = 2.0 * ( xz - yw );
787   const ValueType mzy = 2.0 * ( yz + xw );
788   const ValueType myz = 2.0 * ( yz - xw );
789  
790   result[0] = mxx * v[0] + mxy * v[1] + mxz * v[2];
791   result[1] = myx * v[0] + myy * v[1] + myz * v[2];
792   result[2] = mzx * v[0] + mzy * v[1] + mzz * v[2];
793
794   return result;
795 }
796
797
798 EML
799 EML
800 /**
801  * Transform a VnlVector
802  */
803 template<class T>
804 typename Versor<T>::VnlVectorType 
805 Versor<T>
806 ::Transform( const VnlVectorType & v ) const
807 {
808   VnlVectorType result;
809   
810   const ValueType xx = m_X * m_X;
811   const ValueType yy = m_Y * m_Y;
812   const ValueType zz = m_Z * m_Z;
813   const ValueType xy = m_X * m_Y;
814   const ValueType xz = m_X * m_Z;
815   const ValueType xw = m_X * m_W;
816   const ValueType yz = m_Y * m_Z;
817   const ValueType yw = m_Y * m_W;
818   const ValueType zw = m_Z * m_W;
819
820   const ValueType mxx = 1.0 - 2.0 * ( yy + zz );
821   const ValueType myy = 1.0 - 2.0 * ( xx + zz );
822   const ValueType mzz = 1.0 - 2.0 * ( xx + yy );
823   const ValueType mxy = 2.0 * ( xy - zw );
824   const ValueType mxz = 2.0 * ( xz + yw );
825   const ValueType myx = 2.0 * ( xy + zw );
826   const ValueType mzx = 2.0 * ( xz - yw );
827   const ValueType mzy = 2.0 * ( yz + xw );
828   const ValueType myz = 2.0 * ( yz - xw );
829     
830   result[0] = mxx * v[0] + mxy * v[1] + mxz * v[2];
831   result[1] = myx * v[0] + myy * v[1] + myz * v[2];
832   result[2] = mzx * v[0] + mzy * v[1] + mzz * v[2];
833
834   return result;
835
836 }
837  
838
839
840 EML
841 /**
842  * Get Matrix representation
843  */
844 template<class T>
845 Matrix<T,3,3>
846 Versor<T>
847 ::GetMatrix( void ) const
848 {
849   Matrix<T,3,3> matrix;
850   
851   const ValueType xx = m_X * m_X;
852   const ValueType yy = m_Y * m_Y;
853   const ValueType zz = m_Z * m_Z;
854   const ValueType xy = m_X * m_Y;
855   const ValueType xz = m_X * m_Z;
856   const ValueType xw = m_X * m_W;
857   const ValueType yz = m_Y * m_Z;
858   const ValueType yw = m_Y * m_W;
859   const ValueType zw = m_Z * m_W;
860
861   matrix[0][0] = 1.0 - 2.0 * ( yy + zz );
862   matrix[1][1] = 1.0 - 2.0 * ( xx + zz );
863   matrix[2][2] = 1.0 - 2.0 * ( xx + yy );
864   matrix[0][1] = 2.0 * ( xy - zw );
865   matrix[0][2] = 2.0 * ( xz + yw );
866   matrix[1][0] = 2.0 * ( xy + zw );
867   matrix[2][0] = 2.0 * ( xz - yw );
868   matrix[2][1] = 2.0 * ( yz + xw );
869   matrix[1][2] = 2.0 * ( yz - xw );
870     
871   return matrix;
872
873 }
874  
875
876
877 } // end namespace itk
878
879
880 #endif
881
Generated by KWStyle 1.0b on Tuesday January,17 at 02:14:53PM
© Kitware Inc.