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