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