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