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