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