Cleanup: remove redundant doxygen \file argument
[blender.git] / intern / smoke / intern / VEC3.h
1 /** \file \ingroup smoke
2  */
3 /******************************************************************************
4  * Copyright 2007 Nils Thuerey
5  * Basic vector class 
6  *****************************************************************************/
7 #ifndef BASICVECTOR_H
8 #define BASICVECTOR_H
9
10 #include <math.h>
11 #include <stdlib.h>
12 #include <stdio.h>
13 #include <iostream>
14 #include <sstream>
15
16 // use which fp-precision? 1=float, 2=double
17 #ifndef FLOATINGPOINT_PRECISION
18 #if DDF_DEBUG==1
19 #define FLOATINGPOINT_PRECISION 2
20 #else // DDF_DEBUG==1
21 #define FLOATINGPOINT_PRECISION 1
22 #endif // DDF_DEBUG==1
23 #endif
24
25 // VECTOR_EPSILON is the minimal vector length
26 // In order to be able to discriminate floating point values near zero, and
27 // to be sure not to fail a comparison because of roundoff errors, use this
28 // value as a threshold.  
29
30 #if FLOATINGPOINT_PRECISION==1
31 typedef float Real;
32 #define FP_REAL_MAX __FLT_MAX__
33 #define VECTOR_EPSILON (1e-5f)
34 #else
35 typedef double Real;
36 #define FP_REAL_MAX __DBL_MAX__
37 #define VECTOR_EPSILON (1e-10)
38 #endif
39
40
41 // hardcoded limits for now...
42 // for e.g. MSVC compiler...
43 // some of these defines can be needed
44 // for linux systems as well (e.g. FLT_MAX)
45 #ifndef __FLT_MAX__
46 #       ifdef FLT_MAX  // try to use it instead
47 #               define __FLT_MAX__ FLT_MAX
48 #       else // FLT_MAX
49 #               define __FLT_MAX__ 3.402823466e+38f
50 #       endif // FLT_MAX
51 #endif // __FLT_MAX__
52 #ifndef __DBL_MAX__
53 #       ifdef DBL_MAX // try to use it instead
54 #               define __DBL_MAX__ DBL_MAX
55 #       else // DBL_MAX
56 #               define __DBL_MAX__ 1.7976931348623158e+308
57 #       endif // DBL_MAX
58 #endif // __DBL_MAX__
59
60 #ifndef FLT_MAX
61 #define FLT_MAX __FLT_MAX__
62 #endif
63
64 #ifndef DBL_MAX
65 #define DBL_MAX __DBL_MAX__
66 #endif
67
68 #ifndef M_PI
69 #       define M_PI 3.1415926536
70 #       define M_E  2.7182818284
71 #endif
72
73
74
75 namespace BasicVector {
76
77
78 // basic inlined vector class
79 template<class Scalar>
80 class Vector3Dim
81 {
82 public:
83   // Constructor
84   inline Vector3Dim();
85   // Copy-Constructor
86   inline Vector3Dim(const Vector3Dim<Scalar> &v );
87   inline Vector3Dim(const float *);
88   inline Vector3Dim(const double *);
89   // construct a vector from one Scalar
90   inline Vector3Dim(Scalar);
91   // construct a vector from three Scalars
92   inline Vector3Dim(Scalar, Scalar, Scalar);
93
94         // get address of array for OpenGL
95         Scalar *getAddress() { return value; }
96
97   // Assignment operator
98   inline const Vector3Dim<Scalar>& operator=  (const Vector3Dim<Scalar>& v);
99   // Assignment operator
100   inline const Vector3Dim<Scalar>& operator=  (Scalar s);
101   // Assign and add operator
102   inline const Vector3Dim<Scalar>& operator+= (const Vector3Dim<Scalar>& v);
103   // Assign and add operator
104   inline const Vector3Dim<Scalar>& operator+= (Scalar s);
105   // Assign and sub operator
106   inline const Vector3Dim<Scalar>& operator-= (const Vector3Dim<Scalar>& v);
107   // Assign and sub operator
108   inline const Vector3Dim<Scalar>& operator-= (Scalar s);
109   // Assign and mult operator
110   inline const Vector3Dim<Scalar>& operator*= (const Vector3Dim<Scalar>& v);
111   // Assign and mult operator
112   inline const Vector3Dim<Scalar>& operator*= (Scalar s);
113   // Assign and div operator
114   inline const Vector3Dim<Scalar>& operator/= (const Vector3Dim<Scalar>& v);
115   // Assign and div operator
116   inline const Vector3Dim<Scalar>& operator/= (Scalar s);
117
118
119   // unary operator
120   inline Vector3Dim<Scalar> operator- () const;
121
122   // binary operator add
123   inline Vector3Dim<Scalar> operator+ (const Vector3Dim<Scalar>&) const;
124   // binary operator add
125   inline Vector3Dim<Scalar> operator+ (Scalar) const;
126   // binary operator sub
127   inline Vector3Dim<Scalar> operator- (const Vector3Dim<Scalar>&) const;
128   // binary operator sub
129   inline Vector3Dim<Scalar> operator- (Scalar) const;
130   // binary operator mult
131   inline Vector3Dim<Scalar> operator* (const Vector3Dim<Scalar>&) const;
132   // binary operator mult
133   inline Vector3Dim<Scalar> operator* (Scalar) const;
134   // binary operator div
135   inline Vector3Dim<Scalar> operator/ (const Vector3Dim<Scalar>&) const;
136   // binary operator div
137   inline Vector3Dim<Scalar> operator/ (Scalar) const;
138
139   // Projection normal to a vector
140   inline Vector3Dim<Scalar>       getOrthogonalntlVector3Dim() const;
141   // Project into a plane
142   inline const Vector3Dim<Scalar>& projectNormalTo(const Vector3Dim<Scalar> &v);
143   
144   // minimize
145   inline const Vector3Dim<Scalar> &minimize(const Vector3Dim<Scalar> &);
146   // maximize
147   inline const Vector3Dim<Scalar> &maximize(const Vector3Dim<Scalar> &);
148   
149   // access operator
150   inline Scalar& operator[](unsigned int i);
151   // access operator
152   inline const Scalar& operator[](unsigned int i) const;
153
154         //! actual values
155         union {
156                 struct {
157                 Scalar value[3];  
158                 };
159                 struct {
160                 Scalar x;
161                 Scalar y;
162                 Scalar z;
163                 };
164                 struct {
165                 Scalar X;
166                 Scalar Y;
167                 Scalar Z;
168                 };
169         };
170 protected:
171   
172 };
173
174
175
176
177
178 //------------------------------------------------------------------------------
179 // VECTOR inline FUNCTIONS
180 //------------------------------------------------------------------------------
181
182
183
184 /*************************************************************************
185   Constructor.
186   */
187 template<class Scalar>
188 inline Vector3Dim<Scalar>::Vector3Dim( void )
189 {
190   value[0] = value[1] = value[2] = 0;
191 }
192
193
194
195 /*************************************************************************
196   Copy-Constructor.
197   */
198 template<class Scalar>
199 inline Vector3Dim<Scalar>::Vector3Dim( const Vector3Dim<Scalar> &v )
200 {
201   value[0] = v.value[0];
202   value[1] = v.value[1];
203   value[2] = v.value[2];
204 }
205 template<class Scalar>
206 inline Vector3Dim<Scalar>::Vector3Dim( const float *fvalue)
207 {
208   value[0] = (Scalar)fvalue[0];
209   value[1] = (Scalar)fvalue[1];
210   value[2] = (Scalar)fvalue[2];
211 }
212 template<class Scalar>
213 inline Vector3Dim<Scalar>::Vector3Dim( const double *fvalue)
214 {
215   value[0] = (Scalar)fvalue[0];
216   value[1] = (Scalar)fvalue[1];
217   value[2] = (Scalar)fvalue[2];
218 }
219
220
221
222 /*************************************************************************
223   Constructor for a vector from a single Scalar. All components of
224   the vector get the same value.
225   \param s The value to set
226   \return The new vector
227   */
228 template<class Scalar>
229 inline Vector3Dim<Scalar>::Vector3Dim(Scalar s )
230 {
231   value[0]= s;
232   value[1]= s;
233   value[2]= s;
234 }
235
236
237 /*************************************************************************
238   Constructor for a vector from three Scalars.
239   \param s1 The value for the first vector component
240   \param s2 The value for the second vector component
241   \param s3 The value for the third vector component
242   \return The new vector
243   */
244 template<class Scalar>
245 inline Vector3Dim<Scalar>::Vector3Dim(Scalar s1, Scalar s2, Scalar s3)
246 {
247   value[0]= s1;
248   value[1]= s2;
249   value[2]= s3;
250 }
251
252
253
254 /*************************************************************************
255   Copy a Vector3Dim componentwise.
256   \param v vector with values to be copied
257   \return Reference to self
258   */
259 template<class Scalar>
260 inline const Vector3Dim<Scalar>&
261 Vector3Dim<Scalar>::operator=( const Vector3Dim<Scalar> &v )
262 {
263   value[0] = v.value[0];
264   value[1] = v.value[1];
265   value[2] = v.value[2];  
266   return *this;
267 }
268
269
270 /*************************************************************************
271   Copy a Scalar to each component.
272   \param s The value to copy
273   \return Reference to self
274   */
275 template<class Scalar>
276 inline const Vector3Dim<Scalar>&
277 Vector3Dim<Scalar>::operator=(Scalar s)
278 {
279   value[0] = s;
280   value[1] = s;
281   value[2] = s;  
282   return *this;
283 }
284
285
286 /*************************************************************************
287   Add another Vector3Dim componentwise.
288   \param v vector with values to be added
289   \return Reference to self
290   */
291 template<class Scalar>
292 inline const Vector3Dim<Scalar>&
293 Vector3Dim<Scalar>::operator+=( const Vector3Dim<Scalar> &v )
294 {
295   value[0] += v.value[0];
296   value[1] += v.value[1];
297   value[2] += v.value[2];  
298   return *this;
299 }
300
301
302 /*************************************************************************
303   Add a Scalar value to each component.
304   \param s Value to add
305   \return Reference to self
306   */
307 template<class Scalar>
308 inline const Vector3Dim<Scalar>&
309 Vector3Dim<Scalar>::operator+=(Scalar s)
310 {
311   value[0] += s;
312   value[1] += s;
313   value[2] += s;  
314   return *this;
315 }
316
317
318 /*************************************************************************
319   Subtract another vector componentwise.
320   \param v vector of values to subtract
321   \return Reference to self
322   */
323 template<class Scalar>
324 inline const Vector3Dim<Scalar>&
325 Vector3Dim<Scalar>::operator-=( const Vector3Dim<Scalar> &v )
326 {
327   value[0] -= v.value[0];
328   value[1] -= v.value[1];
329   value[2] -= v.value[2];  
330   return *this;
331 }
332
333
334 /*************************************************************************
335   Subtract a Scalar value from each component.
336   \param s Value to subtract
337   \return Reference to self
338   */
339 template<class Scalar>
340 inline const Vector3Dim<Scalar>&
341 Vector3Dim<Scalar>::operator-=(Scalar s)
342 {
343   value[0]-= s;
344   value[1]-= s;
345   value[2]-= s;  
346   return *this;
347 }
348
349
350 /*************************************************************************
351   Multiply with another vector componentwise.
352   \param v vector of values to multiply with
353   \return Reference to self
354   */
355 template<class Scalar>
356 inline const Vector3Dim<Scalar>&
357 Vector3Dim<Scalar>::operator*=( const Vector3Dim<Scalar> &v )
358 {
359   value[0] *= v.value[0];
360   value[1] *= v.value[1];
361   value[2] *= v.value[2];  
362   return *this;
363 }
364
365
366 /*************************************************************************
367   Multiply each component with a Scalar value.
368   \param s Value to multiply with
369   \return Reference to self
370   */
371 template<class Scalar>
372 inline const Vector3Dim<Scalar>&
373 Vector3Dim<Scalar>::operator*=(Scalar s)
374 {
375   value[0] *= s;
376   value[1] *= s;
377   value[2] *= s;  
378   return *this;
379 }
380
381
382 /*************************************************************************
383   Divide by another Vector3Dim componentwise.
384   \param v vector of values to divide by
385   \return Reference to self
386   */
387 template<class Scalar>
388 inline const Vector3Dim<Scalar>&
389 Vector3Dim<Scalar>::operator/=( const Vector3Dim<Scalar> &v )
390 {
391   value[0] /= v.value[0];
392   value[1] /= v.value[1];
393   value[2] /= v.value[2];  
394   return *this;
395 }
396
397
398 /*************************************************************************
399   Divide each component by a Scalar value.
400   \param s Value to divide by
401   \return Reference to self
402   */
403 template<class Scalar>
404 inline const Vector3Dim<Scalar>&
405 Vector3Dim<Scalar>::operator/=(Scalar s)
406 {
407   value[0] /= s;
408   value[1] /= s;
409   value[2] /= s;
410   return *this;
411 }
412
413
414 //------------------------------------------------------------------------------
415 // unary operators
416 //------------------------------------------------------------------------------
417
418
419 /*************************************************************************
420   Build componentwise the negative this vector.
421   \return The new (negative) vector
422   */
423 template<class Scalar>
424 inline Vector3Dim<Scalar>
425 Vector3Dim<Scalar>::operator-() const
426 {
427   return Vector3Dim<Scalar>(-value[0], -value[1], -value[2]);
428 }
429
430
431
432 //------------------------------------------------------------------------------
433 // binary operators
434 //------------------------------------------------------------------------------
435
436
437 /*************************************************************************
438   Build a vector with another vector added componentwise.
439   \param v The second vector to add
440   \return The sum vector
441   */
442 template<class Scalar>
443 inline Vector3Dim<Scalar>
444 Vector3Dim<Scalar>::operator+( const Vector3Dim<Scalar> &v ) const
445 {
446   return Vector3Dim<Scalar>(value[0]+v.value[0],
447                         value[1]+v.value[1],
448                         value[2]+v.value[2]);
449 }
450
451
452 /*************************************************************************
453   Build a vector with a Scalar value added to each component.
454   \param s The Scalar value to add
455   \return The sum vector
456   */
457 template<class Scalar>
458 inline Vector3Dim<Scalar>
459 Vector3Dim<Scalar>::operator+(Scalar s) const
460 {
461   return Vector3Dim<Scalar>(value[0]+s,
462                         value[1]+s,
463                         value[2]+s);
464 }
465
466
467 /*************************************************************************
468   Build a vector with another vector subtracted componentwise.
469   \param v The second vector to subtract
470   \return The difference vector
471   */
472 template<class Scalar>
473 inline Vector3Dim<Scalar>
474 Vector3Dim<Scalar>::operator-( const Vector3Dim<Scalar> &v ) const
475 {
476   return Vector3Dim<Scalar>(value[0]-v.value[0],
477                         value[1]-v.value[1],
478                         value[2]-v.value[2]);
479 }
480
481
482 /*************************************************************************
483   Build a vector with a Scalar value subtracted componentwise.
484   \param s The Scalar value to subtract
485   \return The difference vector
486   */
487 template<class Scalar>
488 inline Vector3Dim<Scalar>
489 Vector3Dim<Scalar>::operator-(Scalar s ) const
490 {
491   return Vector3Dim<Scalar>(value[0]-s,
492                         value[1]-s,
493                         value[2]-s);
494 }
495
496
497
498 /*************************************************************************
499   Build a vector with another vector multiplied by componentwise.
500   \param v The second vector to muliply with
501   \return The product vector
502   */
503 template<class Scalar>
504 inline Vector3Dim<Scalar>
505 Vector3Dim<Scalar>::operator*( const Vector3Dim<Scalar>& v) const
506 {
507   return Vector3Dim<Scalar>(value[0]*v.value[0],
508                         value[1]*v.value[1],
509                         value[2]*v.value[2]);
510 }
511
512
513 /*************************************************************************
514   Build a Vector3Dim with a Scalar value multiplied to each component.
515   \param s The Scalar value to multiply with
516   \return The product vector
517   */
518 template<class Scalar>
519 inline Vector3Dim<Scalar>
520 Vector3Dim<Scalar>::operator*(Scalar s) const
521 {
522   return Vector3Dim<Scalar>(value[0]*s, value[1]*s, value[2]*s);
523 }
524
525
526 /*************************************************************************
527   Build a vector divided componentwise by another vector.
528   \param v The second vector to divide by
529   \return The ratio vector
530   */
531 template<class Scalar>
532 inline Vector3Dim<Scalar>
533 Vector3Dim<Scalar>::operator/(const Vector3Dim<Scalar>& v) const
534 {
535   return Vector3Dim<Scalar>(value[0]/v.value[0],
536                         value[1]/v.value[1],
537                         value[2]/v.value[2]);
538 }
539
540
541
542 /*************************************************************************
543   Build a vector divided componentwise by a Scalar value.
544   \param s The Scalar value to divide by
545   \return The ratio vector
546   */
547 template<class Scalar>
548 inline Vector3Dim<Scalar>
549 Vector3Dim<Scalar>::operator/(Scalar s) const
550 {
551   return Vector3Dim<Scalar>(value[0]/s,
552                         value[1]/s,
553                         value[2]/s);
554 }
555
556
557
558
559
560 /*************************************************************************
561   Get a particular component of the vector.
562   \param i Number of Scalar to get
563   \return Reference to the component
564   */
565 template<class Scalar>
566 inline Scalar&
567 Vector3Dim<Scalar>::operator[]( unsigned int i )
568 {
569   return value[i];
570 }
571
572
573 /*************************************************************************
574   Get a particular component of a constant vector.
575   \param i Number of Scalar to get
576   \return Reference to the component
577   */
578 template<class Scalar>
579 inline const Scalar&
580 Vector3Dim<Scalar>::operator[]( unsigned int i ) const
581 {
582   return value[i];
583 }
584
585
586
587 //------------------------------------------------------------------------------
588 // BLITZ compatibility functions
589 //------------------------------------------------------------------------------
590
591
592
593 /*************************************************************************
594   Compute the scalar product with another vector.
595   \param v The second vector to work with
596   \return The value of the scalar product
597   */
598 template<class Scalar>
599 inline Scalar dot(const Vector3Dim<Scalar> &t, const Vector3Dim<Scalar> &v )
600 {
601   //return t.value[0]*v.value[0] + t.value[1]*v.value[1] + t.value[2]*v.value[2];
602   return ((t[0]*v[0]) + (t[1]*v[1]) + (t[2]*v[2]));
603 }
604
605
606 /*************************************************************************
607   Calculate the cross product of this and another vector
608  */
609 template<class Scalar>
610 inline Vector3Dim<Scalar> cross(const Vector3Dim<Scalar> &t, const Vector3Dim<Scalar> &v)
611 {
612   Vector3Dim<Scalar> cp( 
613                         ((t[1]*v[2]) - (t[2]*v[1])),
614                   ((t[2]*v[0]) - (t[0]*v[2])),
615                   ((t[0]*v[1]) - (t[1]*v[0])) );
616   return cp;
617 }
618
619
620
621
622 /*************************************************************************
623   Compute a vector that is orthonormal to self. Nothing else can be assumed
624   for the direction of the new vector.
625   \return The orthonormal vector
626   */
627 template<class Scalar>
628 Vector3Dim<Scalar>
629 Vector3Dim<Scalar>::getOrthogonalntlVector3Dim() const
630 {
631   // Determine the  component with max. absolute value
632   int max= (fabs(value[0]) > fabs(value[1])) ? 0 : 1;
633   max= (fabs(value[max]) > fabs(value[2])) ? max : 2;
634
635   /*************************************************************************
636     Choose another axis than the one with max. component and project
637     orthogonal to self
638     */
639   Vector3Dim<Scalar> vec(0.0);
640   vec[(max+1)%3]= 1;
641   vec.normalize();
642   vec.projectNormalTo(this->getNormalized());
643   return vec;
644 }
645
646
647 /*************************************************************************
648   Projects the vector into a plane normal to the given vector, which must
649   have unit length. Self is modified.
650   \param v The plane normal
651   \return The projected vector
652   */
653 template<class Scalar>
654 inline const Vector3Dim<Scalar>&
655 Vector3Dim<Scalar>::projectNormalTo(const Vector3Dim<Scalar> &v)
656 {
657   Scalar sprod = dot(*this,v);
658   value[0]= value[0] - v.value[0] * sprod;
659   value[1]= value[1] - v.value[1] * sprod;
660   value[2]= value[2] - v.value[2] * sprod;  
661   return *this;
662 }
663
664
665
666 //------------------------------------------------------------------------------
667 // Other helper functions
668 //------------------------------------------------------------------------------
669
670
671
672 /*************************************************************************
673   Minimize the vector, i.e. set each entry of the vector to the minimum
674   of both values.
675   \param pnt The second vector to compare with
676   \return Reference to the modified self
677   */
678 template<class Scalar>
679 inline const Vector3Dim<Scalar> &
680 Vector3Dim<Scalar>::minimize(const Vector3Dim<Scalar> &pnt)
681 {
682   for (unsigned int i = 0; i < 3; i++)
683     value[i] = MIN(value[i],pnt[i]);
684   return *this;
685 }
686
687
688
689 /*************************************************************************
690   Maximize the vector, i.e. set each entry of the vector to the maximum
691   of both values.
692   \param pnt The second vector to compare with
693   \return Reference to the modified self
694   */
695 template<class Scalar>
696 inline const Vector3Dim<Scalar> &
697 Vector3Dim<Scalar>::maximize(const Vector3Dim<Scalar> &pnt)
698 {
699   for (unsigned int i = 0; i < 3; i++)
700     value[i] = MAX(value[i],pnt[i]);
701   return *this;
702 }
703
704
705
706
707
708
709 /************************************************************************/
710 // HELPER FUNCTIONS, independent of implementation
711 /************************************************************************/
712
713 #define VECTOR_TYPE Vector3Dim<Scalar>
714
715
716 /*************************************************************************
717   Compute the length (norm) of the vector.
718   \return The value of the norm
719   */
720 template<class Scalar>
721 inline Scalar norm( const VECTOR_TYPE &v)
722 {
723   Scalar l = v[0]*v[0] + v[1]*v[1] + v[2]*v[2];
724   return (fabs(l-1.) < VECTOR_EPSILON*VECTOR_EPSILON) ? 1. : sqrt(l);
725 }
726
727 // for e.g. min max operator
728 inline Real normHelper(const Vector3Dim<Real> &v) {
729         return norm(v);
730 }       
731 inline Real normHelper(const Real &v) {
732         return (0.0f < v) ? v : -v ;
733 }       
734 inline Real normHelper(const int &v) {
735         return (0 < v) ? (Real)(v) : (Real)(-v) ; 
736 }       
737
738
739 /*************************************************************************
740   Same as getNorm but doesnt sqrt  
741   */
742 template<class Scalar>
743 inline Scalar normNoSqrt( const VECTOR_TYPE &v)
744 {
745   return v[0]*v[0] + v[1]*v[1] + v[2]*v[2];
746 }
747
748
749 /*************************************************************************
750   Compute a normalized vector based on this vector.
751   \return The new normalized vector
752   */
753 template<class Scalar>
754 inline VECTOR_TYPE getNormalized( const VECTOR_TYPE &v)
755 {
756   Scalar l = v[0]*v[0] + v[1]*v[1] + v[2]*v[2];
757   if (fabs(l-1.) < VECTOR_EPSILON*VECTOR_EPSILON)
758     return v; /* normalized "enough"... */
759   else if (l > VECTOR_EPSILON*VECTOR_EPSILON)
760   {
761     Scalar fac = 1./sqrt(l);
762     return VECTOR_TYPE(v[0]*fac, v[1]*fac, v[2]*fac);
763   }
764   else
765     return VECTOR_TYPE((Scalar)0);
766 }
767
768
769 /*************************************************************************
770   Compute the norm of the vector and normalize it.
771   \return The value of the norm
772   */
773 template<class Scalar>
774 inline Scalar normalize( VECTOR_TYPE &v) 
775 {
776   Scalar norm;
777   Scalar l = v[0]*v[0] + v[1]*v[1] + v[2]*v[2];  
778   if (fabs(l-1.) < VECTOR_EPSILON*VECTOR_EPSILON) {
779     norm = 1.;
780         } else if (l > VECTOR_EPSILON*VECTOR_EPSILON) {
781     norm = sqrt(l);
782     Scalar fac = 1./norm;
783     v[0] *= fac;
784     v[1] *= fac;
785     v[2] *= fac; 
786         } else {
787     v[0]= v[1]= v[2]= 0;
788     norm = 0.;
789   }
790   return (Scalar)norm;
791 }
792
793
794 /*************************************************************************
795   Compute a vector, that is self (as an incoming
796   vector) reflected at a surface with a distinct normal vector. Note
797   that the normal is reversed, if the scalar product with it is positive.
798   \param n The surface normal
799   \return The new reflected vector
800   */
801 template<class Scalar>
802 inline VECTOR_TYPE reflectVector(const VECTOR_TYPE &t, const VECTOR_TYPE &n) 
803 {
804   VECTOR_TYPE nn= (dot(t, n) > 0.0) ? (n*-1.0) : n;
805   return ( t - nn * (2.0 * dot(nn, t)) );
806 }
807
808
809
810 /*************************************************************************
811  * My own refraction calculation
812  * Taken from Glassner's book, section 5.2 (Heckberts method)
813  */
814 template<class Scalar>
815 inline VECTOR_TYPE refractVector(const VECTOR_TYPE &t, const VECTOR_TYPE &normal, Scalar nt, Scalar nair, int &refRefl) 
816 {
817         Scalar eta = nair / nt;
818         Scalar n = -dot(t, normal);
819         Scalar tt = 1.0 + eta*eta* (n*n-1.0);
820         if(tt<0.0) {
821                 // we have total reflection!
822                 refRefl = 1;
823         } else {
824                 // normal reflection
825                 tt = eta*n - sqrt(tt);
826                 return( t*eta + normal*tt );
827         }
828         return t;
829 }
830
831
832 /*************************************************************************
833   Test two ntlVector3Dims for equality based on the equality of their
834   values within a small threshold.
835   \param c The second vector to compare
836   \return TRUE if both are equal
837   \sa getEpsilon()
838   */
839 template<class Scalar>
840 inline bool equal(const VECTOR_TYPE &v, const VECTOR_TYPE &c)
841 {
842   return (ABS(v[0]-c[0]) + 
843           ABS(v[1]-c[1]) + 
844           ABS(v[2]-c[2]) < VECTOR_EPSILON);
845 }
846
847
848 /*************************************************************************
849  * Assume this vector is an RGB color, and convert it to HSV
850  */
851 template<class Scalar>
852 inline void rgbToHsv( VECTOR_TYPE &V )
853 {
854         Scalar h=0,s=0,v=0;
855         Scalar maxrgb, minrgb, delta;
856         // convert to hsv...
857         maxrgb = V[0];
858         int maxindex = 1;
859         if(V[2] > maxrgb){ maxrgb = V[2]; maxindex = 2; }
860         if(V[1] > maxrgb){ maxrgb = V[1]; maxindex = 3; }
861         minrgb = V[0];
862         if(V[2] < minrgb) minrgb = V[2];
863         if(V[1] < minrgb) minrgb = V[1];
864
865         v = maxrgb;
866         delta = maxrgb-minrgb;
867
868         if(maxrgb > 0) s = delta/maxrgb;
869         else s = 0;
870
871         h = 0;
872         if(s > 0) {
873                 if(maxindex == 1) {
874                         h = ((V[1]-V[2])/delta)  + 0.0; }
875                 if(maxindex == 2) {
876                         h = ((V[2]-V[0])/delta)  + 2.0; }
877                 if(maxindex == 3) {
878                         h = ((V[0]-V[1])/delta)  + 4.0; }
879                 h *= 60.0;
880                 if(h < 0.0) h += 360.0;
881         }
882
883         V[0] = h;
884         V[1] = s;
885         V[2] = v;
886 }
887
888 /*************************************************************************
889  * Assume this vector is HSV and convert to RGB
890  */
891 template<class Scalar>
892 inline void hsvToRgb( VECTOR_TYPE &V )
893 {
894         Scalar h = V[0], s = V[1], v = V[2];
895         Scalar r=0,g=0,b=0;
896         Scalar p,q,t, fracth;
897         int floorh;
898         // ...and back to rgb
899         if(s == 0) {
900                 r = g = b = v; }
901         else {
902                 h /= 60.0;
903                 floorh = (int)h;
904                 fracth = h - floorh;
905                 p = v * (1.0 - s);
906                 q = v * (1.0 - (s * fracth));
907                 t = v * (1.0 - (s * (1.0 - fracth)));
908                 switch (floorh) {
909                 case 0: r = v; g = t; b = p; break;
910                 case 1: r = q; g = v; b = p; break;
911                 case 2: r = p; g = v; b = t; break;
912                 case 3: r = p; g = q; b = v; break;
913                 case 4: r = t; g = p; b = v; break;
914                 case 5: r = v; g = p; b = q; break;
915                 }
916         }
917
918         V[0] = r;
919         V[1] = g;
920         V[2] = b;
921 }
922
923 //------------------------------------------------------------------------------
924 // STREAM FUNCTIONS
925 //------------------------------------------------------------------------------
926
927
928
929 //! global string for formatting vector output in utilities.cpp
930 //extern const char *globVecFormatStr;
931 #if 0
932 static const char *globVecFormatStr = "[%6.4f,%6.4f,%6.4f]";
933 #endif
934
935 /*************************************************************************
936   Outputs the object in human readable form using the format
937   [x,y,z]
938   */
939 template<class Scalar>
940 std::ostream&
941 operator<<( std::ostream& os, const BasicVector::Vector3Dim<Scalar>& i )
942 {
943 #if 0
944         char buf[256];
945 #  if _WIN32
946         sprintf(buf,globVecFormatStr, (double)i[0],(double)i[1],(double)i[2]);
947 #  else
948         snprintf(buf,256,globVecFormatStr, (double)i[0],(double)i[1],(double)i[2]);
949 #  endif
950         os << std::string(buf); 
951 #else
952         (void)i;  /* Ignored. */
953 #endif
954         return os;
955 }
956
957
958 /*************************************************************************
959   Reads the contents of the object from a stream using the same format
960   as the output operator.
961   */
962 template<class Scalar>
963 std::istream&
964 operator>>( std::istream& is, BasicVector::Vector3Dim<Scalar>& i )
965 {
966   char c;
967   char dummy[3];
968   is >> c >> i[0] >> dummy >> i[1] >> dummy >> i[2] >> c;
969   return is;
970 }
971
972
973 /**************************************************************************/
974 // typedefs!
975 /**************************************************************************/
976
977 /* get minimal vector length value that can be discriminated.  */
978 inline Real getVecEpsilon() { return (Real)VECTOR_EPSILON; }
979
980 // a 3D integer vector
981 typedef Vector3Dim<int> Vec3Int; 
982
983 // a 3D vector 
984 typedef Vector3Dim<Real> Vec3; 
985
986
987 }; // namespace 
988
989
990 #endif /* BASICVECTOR_H */