Merge branch 'master' into blender2.8
[blender.git] / intern / elbeem / intern / ntl_bsptree.cpp
1 /** \file elbeem/intern/ntl_bsptree.cpp
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  * Tree container for fast triangle intersects
10  *
11  *****************************************************************************/
12
13
14 #include "ntl_bsptree.h"
15 #include "utilities.h"
16
17 #include <algorithm>
18
19 /*! Static global variable for sorting direction */
20 int globalSortingAxis;
21 /*! Access to points array for sorting */
22 vector<ntlVec3Gfx> *globalSortingPoints;
23
24 #define TREE_DOUBLEI 300
25
26 /* try axis selection? */
27 bool chooseAxis = 0;
28 /* do median search? */
29 int doSort = 0;
30
31
32 //! struct for a single node in the bsp tree
33 class BSPNode {
34         public:
35                 BSPNode() {};
36
37                 ntlVec3Gfx min,max;              /* AABB for node */
38                 vector<ntlTriangle *> *members;  /* stored triangles */
39                 BSPNode *child[2]; /* pointer to children nodes */
40                 char axis;                  /* division axis */
41                 char cloneVec;              /* is this vector a clone? */
42
43                 //! check if node is a leaf
44                 inline bool isLeaf() const { 
45                         return (child[0] == NULL); 
46                 }
47 };
48
49
50 //! an element node stack
51 class BSPStackElement {
52         public:
53                 //! tree node
54                 BSPNode *node;
55                 //! min and maximum distance along axis
56                 gfxReal mindist, maxdist;
57 };
58
59 //! bsp tree stack
60 class BSPStack {
61         public:
62                 //! current stack element
63                 int stackPtr;
64                 //! stack storage
65                 BSPStackElement elem[ BSP_STACK_SIZE ];
66 };
67
68 //! triangle bounding box for quick tree subdivision
69 class TriangleBBox {
70         public:
71                 //! start and end of triangle bounding box
72                 ntlVec3Gfx start, end;
73 };
74
75
76 /******************************************************************************
77  * calculate tree statistics
78  *****************************************************************************/
79 void calcStats(BSPNode *node, int depth, int &noLeafs, gfxReal &avgDepth, gfxReal &triPerLeaf,int &totalTris)
80 {
81         if(node->members != NULL) {
82                 totalTris += node->members->size();
83         }
84         //depth = 15; // DBEUG!
85
86         if( (node->child[0]==NULL) && (node->child[1]==NULL) ) {
87                 // leaf
88                 noLeafs++;
89                 avgDepth += depth;
90                 triPerLeaf += node->members->size();
91         } else {
92                 for(int i=0;i<2;i++) 
93                 calcStats(node->child[i], depth+1, noLeafs, avgDepth, triPerLeaf, totalTris);
94         }
95 }
96
97
98
99 /******************************************************************************
100  * triangle comparison function for median search 
101  *****************************************************************************/
102 bool lessTriangleAverage(const ntlTriangle *x, const ntlTriangle *y)
103 {
104         return x->getAverage(globalSortingAxis) < y->getAverage(globalSortingAxis);
105 }
106
107
108 /******************************************************************************
109  * triangle AABB intersection
110  *****************************************************************************/
111 bool ntlTree::checkAABBTriangle(ntlVec3Gfx &min, ntlVec3Gfx &max, ntlTriangle *tri)
112 {
113         // test only BB of triangle
114         TriangleBBox *bbox = &mpTBB[ tri->getBBoxId() ];
115         if( bbox->end[0]   < min[0] ) return false;
116         if( bbox->start[0] > max[0] ) return false;
117         if( bbox->end[1]   < min[1] ) return false;
118         if( bbox->start[1] > max[1] ) return false;
119         if( bbox->end[2]   < min[2] ) return false;
120         if( bbox->start[2] > max[2] ) return false;
121         return true;
122 }
123
124
125
126
127
128
129
130 /******************************************************************************
131  * Default constructor
132  *****************************************************************************/
133 ntlTree::ntlTree() :
134   mStart(0.0), mEnd(0.0), mMaxDepth( 5 ), mMaxListLength( 5 ), mpRoot( NULL) ,
135   mpNodeStack( NULL), mpVertices( NULL ), mpVertNormals( NULL ), mpTriangles( NULL ),
136   mCurrentDepth(0), mCurrentNodes(0), mTriDoubles(0)
137 {
138   errFatal( "ntlTree","Uninitialized BSP Tree!\n",SIMWORLD_INITERROR );
139   return;
140 }
141
142
143 /******************************************************************************
144  * Constructor with init
145  *****************************************************************************/
146 //ntlTree::ntlTree(int depth, int objnum, vector<ntlVec3Gfx> *vertices, vector<ntlVec3Gfx> *normals, vector<ntlTriangle> *trilist) :
147 ntlTree::ntlTree(int depth, int objnum, ntlScene *scene, int triFlagMask) :
148   mStart(0.0), mEnd(0.0), mMaxDepth( depth ), mMaxListLength( objnum ), mpRoot( NULL) ,
149   mpNodeStack( NULL), mpTBB( NULL ),
150         mTriangleMask( 0xFFFF ),
151   mCurrentDepth(0), mCurrentNodes(0), mTriDoubles(0)
152 {  
153         // init scene data pointers
154         mpVertices = scene->getVertexPointer();
155         mpVertNormals = scene->getVertexNormalPointer();
156         mpTriangles = scene->getTrianglePointer();
157         mTriangleMask = triFlagMask;
158
159   if(mpTriangles == NULL) {
160     errFatal( "ntlTree Cons","no triangle list!\n",SIMWORLD_INITERROR);
161     return;
162   }
163   if(mpTriangles->size() == 0) {
164     warnMsg( "ntlTree::ntlTree","No triangles ("<< mpTriangles->size()  <<")!\n");
165                 mStart = mEnd = ntlVec3Gfx(0,0,0);
166     return;
167   }
168   if(depth>=BSP_STACK_SIZE) {
169     errFatal( "ntlTree::ntlTree","Depth to high ("<< mMaxDepth  <<")!\n", SIMWORLD_INITERROR );
170     return;
171   }
172
173   /* check triangles (a bit inefficient, but we dont know which vertices belong
174      to this tree), and generate bounding boxes */
175         mppTriangles = new vector<ntlTriangle *>;
176         int noOfTriangles = mpTriangles->size();
177         mpTBB = new TriangleBBox[ noOfTriangles ];
178         int bbCount = 0;
179   mStart = mEnd = (*mpVertices)[ mpTriangles->front().getPoints()[0] ];
180         //errMsg("TreeDebug","Start");
181   for (vector<ntlTriangle>::iterator iter = mpTriangles->begin();
182        iter != mpTriangles->end(); 
183        iter++ ) {
184                 //errorOut(" d "<< convertFlags2String((int)(*iter).getFlags()) <<" "<< convertFlags2String( (int)mTriangleMask)<<" add? "<<( ((int)(*iter).getFlags() & (int)mTriangleMask) != 0 ) );
185                 // discard triangles that dont match mask
186                 if( ((int)(*iter).getFlags() & (int)mTriangleMask) == 0 ) {
187                         continue;
188                 }
189
190                 // test? TODO
191                 ntlVec3Gfx tnormal = (*mpVertNormals)[ (*iter).getPoints()[0] ]+
192                         (*mpVertNormals)[ (*iter).getPoints()[1] ]+
193                         (*mpVertNormals)[ (*iter).getPoints()[2] ];
194                 ntlVec3Gfx triangleNormal = (*iter).getNormal();
195                 if( equal(triangleNormal, ntlVec3Gfx(0.0)) ) continue;
196                 if( equal(       tnormal, ntlVec3Gfx(0.0)) ) continue;
197                 // */
198
199                 ntlVec3Gfx bbs, bbe;
200                 //errMsg("TreeDebug","Triangle");
201                 for(int i=0;i<3;i++) {
202                         int index = (*iter).getPoints()[i];
203                         ntlVec3Gfx tp = (*mpVertices)[ index ];
204                         //errMsg("TreeDebug","  Point "<<i<<" = "<<tp<<" ");
205                         if(tp[0] < mStart[0]) mStart[0]= tp[0];
206                         if(tp[0] > mEnd[0])   mEnd[0]= tp[0];
207                         if(tp[1] < mStart[1]) mStart[1]= tp[1];
208                         if(tp[1] > mEnd[1])   mEnd[1]= tp[1];
209                         if(tp[2] < mStart[2]) mStart[2]= tp[2];
210                         if(tp[2] > mEnd[2])   mEnd[2]= tp[2];
211                         if(i==0) {
212                                 bbs = bbe = tp; 
213                         } else {
214                                 if( tp[0] < bbs[0] ) bbs[0] = tp[0];
215                                 if( tp[0] > bbe[0] ) bbe[0] = tp[0];
216                                 if( tp[1] < bbs[1] ) bbs[1] = tp[1];
217                                 if( tp[1] > bbe[1] ) bbe[1] = tp[1];
218                                 if( tp[2] < bbs[2] ) bbs[2] = tp[2];
219                                 if( tp[2] > bbe[2] ) bbe[2] = tp[2];
220                         }
221                 }
222                 mppTriangles->push_back( &(*iter) );
223                 //errMsg("TreeDebug","Triangle "<<(*mpVertices)[(*iter).getPoints()[0]]<<" "<<(*mpVertices)[(*iter).getPoints()[1]]<<" "<<(*mpVertices)[(*iter).getPoints()[2]]<<" ");
224
225                 // add BB
226                 mpTBB[ bbCount ].start = bbs;
227                 mpTBB[ bbCount ].end = bbe;
228                 (*iter).setBBoxId( bbCount );
229                 bbCount++;
230   }
231         
232         
233
234   /* slighlty enlarge bounding tolerance for tree 
235      to avoid problems with triangles paralell to slabs */
236   mStart -= ntlVec3Gfx( getVecEpsilon() );
237   mEnd   += ntlVec3Gfx( getVecEpsilon() );
238
239   /* init root node and stack */
240   mpNodeStack = new BSPStack;
241   mpRoot = new BSPNode;
242   mpRoot->min = mStart;
243   mpRoot->max = mEnd;
244   mpRoot->axis = AXIS_X;
245   mpRoot->members = mppTriangles;
246         mpRoot->child[0] = mpRoot->child[1] = NULL;
247         mpRoot->cloneVec = 0;
248         globalSortingPoints = mpVertices;
249         mpTriDist = new char[ mppTriangles->size() ];
250         mNumNodes = 1;
251         mAbortSubdiv = 0;
252
253   /* create tree */
254   debugOutInter( "Generating BSP Tree...  (Nodes "<< mCurrentNodes <<
255                                                 ", Depth "<<mCurrentDepth<< ") ", 2, 2000 );
256   subdivide(mpRoot, 0, AXIS_X);
257   debMsgStd("ntlTree::ntlTree",DM_MSG,"Generated Tree: Nodes "<< mCurrentNodes <<
258                                                          ", Depth "<<mCurrentDepth<< " with "<<noOfTriangles<<" triangles", 2 );
259
260         delete [] mpTriDist;
261         delete [] mpTBB;
262         mpTriDist = NULL;
263         mpTBB = NULL;
264
265         /* calculate some stats about tree */
266         int noLeafs = 0;
267         gfxReal avgDepth = 0.0;
268         gfxReal triPerLeaf = 0.0;
269         int totalTris = 0;
270         
271         calcStats(mpRoot,0, noLeafs, avgDepth, triPerLeaf, totalTris);
272         avgDepth /= (gfxReal)noLeafs;
273         triPerLeaf /= (gfxReal)noLeafs;
274         debMsgStd("ntlTree::ntlTree",DM_MSG,"Tree ("<<doSort<<","<<chooseAxis<<") Stats: Leafs:"<<noLeafs<<", avgDepth:"<<avgDepth<<
275                         ", triPerLeaf:"<<triPerLeaf<<", triDoubles:"<<mTriDoubles<<", totalTris:"<<totalTris
276                         <<" nodes:"<<mNumNodes
277                         //<<" T"<< (totalTris%3)  // 0=ich, 1=f, 2=a
278                         , 2 );
279
280         if(mAbortSubdiv) {
281                 errMsg("ntlTree::ntlTree","Aborted... "<<mNumNodes);
282         deleteNode(mpRoot);
283                 mpRoot = NULL;
284         }
285 }
286
287 /******************************************************************************
288  * Destructor
289  *****************************************************************************/
290 ntlTree::~ntlTree()
291 {
292   /* delete tree, and all members except for the root node */
293   deleteNode(mpRoot);
294   if(mpNodeStack) delete mpNodeStack;
295 }
296
297
298 /******************************************************************************
299  * subdivide tree
300  *****************************************************************************/
301 void ntlTree::subdivide(BSPNode *node, int depth, int axis)
302 {
303   int nextAxis=0; /* next axis to partition */
304         int allTriDistSet = (1<<0)|(1<<1); // all mpTriDist flags set?
305         //errorOut(" "<<node<<" depth:"<<depth<<" m:"<<node->members->size() <<"  "<<node->min<<" - "<<node->max );
306
307   if(depth>mCurrentDepth) mCurrentDepth = depth;
308   node->child[0] = node->child[1] = NULL;
309         if( ( (int)node->members->size() > mMaxListLength) &&
310                         (depth < mMaxDepth ) 
311                         && (node->cloneVec<10)
312                         && (!mAbortSubdiv)
313                         ) {
314
315                 gfxReal planeDiv = 0.499999;    // position of plane division
316
317                 // determine next subdivision axis
318                 int newaxis = 0;
319                 gfxReal extX = node->max[0]-node->min[0];
320                 gfxReal extY = node->max[1]-node->min[1];
321                 gfxReal extZ = node->max[2]-node->min[2];
322
323                 if( extY>extX  ) {
324                         if( extZ>extY ) {
325                                 newaxis = 2;
326                         } else {
327                                 newaxis = 1;
328                         }
329                 } else {
330                         if( extZ>extX ) {
331                                 newaxis = 2;
332                         } else {
333                                 newaxis = 0;
334                         }
335                 }
336                 axis = node->axis = newaxis;
337
338                 // init child nodes
339                 for( int i=0; i<2; i++) {
340                         /* status output */
341                         mCurrentNodes++;
342                         if(mCurrentNodes % 13973 ==0) {
343                                 debugOutInter( "NTL Generating BSP Tree ("<<doSort<<","<<chooseAxis<<") ...  (Nodes "<< mCurrentNodes <<
344                                                 ", Depth "<<mCurrentDepth<< ") " , 2, 2000);
345                         }
346
347                         /* create new node */
348                         node->child[i] = new BSPNode;
349                         node->child[i]->min = node->min;
350                         node->child[i]->max = node->max;
351                         node->child[i]->max = node->max;
352                         node->child[i]->child[0] = NULL;
353                         node->child[i]->child[1] = NULL;
354                         node->child[i]->members = NULL;
355                         nextAxis = (axis+1)%3;
356                         node->child[i]->axis = nextAxis;
357                         mNumNodes++;
358                         // abort when using 256MB only for tree...
359                         if(mNumNodes*sizeof(BSPNode)> 1024*1024*512) mAbortSubdiv = 1;
360
361                         /* current division plane */
362                         if(!i) {
363                                 node->child[i]->min[axis] = node->min[axis];
364                                 node->child[i]->max[axis] = node->min[axis] + planeDiv*
365                                         (node->max[axis]-node->min[axis]);
366                         } else {
367                                 node->child[i]->min[axis] = node->min[axis] + planeDiv*
368                                         (node->max[axis]-node->min[axis]);
369                                 node->child[i]->max[axis] = node->max[axis];
370                         }
371                 }
372
373
374                 /* process the two children */
375                 int thisTrisFor[2] = {0,0};
376                 int thisTriDoubles[2] = {0,0};
377                 for(int t=0;t<(int)node->members->size();t++) mpTriDist[t] = 0;
378                 for( int i=0; i<2; i++) {
379                         /* distribute triangles */
380                         int t  = 0;
381                         for (vector<ntlTriangle *>::iterator iter = node->members->begin();
382                                         iter != node->members->end(); iter++ ) {
383
384                                 /* add triangle, check bounding box axis */
385                                 TriangleBBox *bbox = &mpTBB[ (*iter)->getBBoxId() ];
386                                 bool isintersect = true;
387                                 if( bbox->end[axis]   < node->child[i]->min[axis] ) isintersect = false;
388                                 else if( bbox->start[axis] > node->child[i]->max[axis] ) isintersect = false;
389                                 if(isintersect) {
390                                         // add flag to vector 
391                                         mpTriDist[t] |= (1<<i);
392                                         // count no. of triangles for vector init
393                                         thisTrisFor[i]++;
394                                 }
395
396                                 if(mpTriDist[t] == allTriDistSet) {
397                                         thisTriDoubles[i]++;
398                                         mTriDoubles++; // TODO check for small geo tree??
399                                 }
400                                 t++;
401                         } /* end of loop over all triangles */
402                 } // i
403
404                 /* distribute triangles */
405                 bool haveCloneVec[2] = {false, false};
406                 for( int i=0; i<2; i++) {
407                         node->child[i]->members = new vector<ntlTriangle *>( thisTrisFor[i] );
408                         node->child[i]->cloneVec = 0;
409                 }
410
411                 int tind0 = 0;
412                 int tind1 = 0;
413                 if( (!haveCloneVec[0]) || (!haveCloneVec[1]) ){
414                         int t  = 0; // triangle index counter
415                         for (vector<ntlTriangle *>::iterator iter = node->members->begin();
416                                         iter != node->members->end(); iter++ ) {
417                                 if(!haveCloneVec[0]) {
418                                         if( (mpTriDist[t] & 1) == 1) {
419                                                 (*node->child[0]->members)[tind0] = (*iter); // dont use push_back for preinited size!
420                                                 tind0++;
421                                         }
422                                 }
423                                 if(!haveCloneVec[1]) {
424                                         if( (mpTriDist[t] & 2) == 2) {
425                                                 (*node->child[1]->members)[tind1] = (*iter); // dont use push_back for preinited size!
426                                                 tind1++;
427                                         }
428                                 }
429                                 t++;
430                         } /* end of loop over all triangles */
431                 }
432
433                 // subdivide children
434                 for( int i=0; i<2; i++) {
435                         /* recurse */
436                         subdivide( node->child[i], depth+1, nextAxis );
437                 }
438
439                 /* if we are here, this are childs, so we dont need the members any more... */
440                 /* delete unecessary members */
441                 if( (!haveCloneVec[0]) && (!haveCloneVec[1]) && (node->cloneVec == 0) ){ 
442                         delete node->members; 
443                 }
444                 node->members = NULL;
445
446         } /* subdivision necessary */
447 }
448
449 /******************************************************************
450  * triangle intersection with triangle pointer,
451  * returns t,u,v by references 
452  */
453 #if GFX_PRECISION==1
454 // float values
455 //! the minimal triangle determinant length
456 #define RAY_TRIANGLE_EPSILON (1e-07)
457
458 #else 
459 // double values
460 //! the minimal triangle determinant length
461 #define RAY_TRIANGLE_EPSILON (1e-15)
462
463 #endif 
464
465
466 /******************************************************************************
467  * intersect ray with BSPtree
468  *****************************************************************************/
469 inline void ntlRay::intersectTriangle(vector<ntlVec3Gfx> *mpV, ntlTriangle *tri, gfxReal &t, gfxReal &u, gfxReal &v) const
470 {
471   /* (cf. moeller&haines, page 305) */
472   t = GFX_REAL_MAX;
473   ntlVec3Gfx  e0 = (*mpV)[ tri->getPoints()[0] ];
474   ntlVec3Gfx  e1 = (*mpV)[ tri->getPoints()[1] ] - e0;
475   ntlVec3Gfx  e2 = (*mpV)[ tri->getPoints()[2] ] - e0;
476   ntlVec3Gfx  p  = cross( mDirection, e2 );
477   gfxReal divisor  = dot(e1, p);        
478   if((divisor > -RAY_TRIANGLE_EPSILON)&&(divisor < RAY_TRIANGLE_EPSILON)) return;
479       
480   gfxReal invDivisor  = 1/divisor;
481   ntlVec3Gfx  s  = mOrigin - e0;
482   u  = invDivisor * dot(s, p);
483   if( (u<0.0-RAY_TRIANGLE_EPSILON) || (u>1.0+RAY_TRIANGLE_EPSILON) ) return;
484       
485   ntlVec3Gfx  q  = cross( s,e1 );
486   v  = invDivisor * dot(mDirection, q);
487   if( (v<0.0-RAY_TRIANGLE_EPSILON) || ((u+v)>1.0+RAY_TRIANGLE_EPSILON) ) return;
488       
489   t = invDivisor * dot(e2, q);      
490 }
491 void ntlTree::intersect(const ntlRay &ray, gfxReal &distance, 
492                 ntlVec3Gfx &normal, 
493                 ntlTriangle *&tri, 
494                 int flags, bool forceNonsmooth) const
495 {
496   gfxReal mint = GFX_REAL_MAX;  /* current minimal t */
497   ntlVec3Gfx  retnormal;       /* intersection (interpolated) normal */
498         gfxReal mintu=0.0, mintv=0.0;    /* u,v for min t intersection */
499
500   BSPNode *curr, *nearChild, *farChild; /* current node and children */
501   gfxReal  planedist, mindist, maxdist;
502   ntlVec3Gfx   pos;
503
504         ntlTriangle *hit = NULL;
505         tri = NULL;
506
507   ray.intersectCompleteAABB(mStart,mEnd,mindist,maxdist);
508
509   if((maxdist < 0.0) ||
510                  (!mpRoot) ||
511      (mindist == GFX_REAL_MAX) ||
512      (maxdist == GFX_REAL_MAX) ) {
513     distance = -1.0;
514     return;
515   }
516   mindist -= getVecEpsilon();
517   maxdist += getVecEpsilon();
518
519   /* stack init */
520   mpNodeStack->elem[0].node = NULL;
521   mpNodeStack->stackPtr = 1;
522
523   curr = mpRoot;  
524   mint = GFX_REAL_MAX;
525   while(curr != NULL) {
526
527     while( !curr->isLeaf() ) {
528       planedist = distanceToPlane(curr, curr->child[0]->max, ray );
529       getChildren(curr, ray.getOrigin(), nearChild, farChild );
530
531                         // check ray direction for small plane distances
532       if( (planedist>-getVecEpsilon() )&&(planedist< getVecEpsilon() ) ) {
533                                 // ray origin on intersection plane
534                                 planedist = 0.0;
535                                 if(ray.getDirection()[curr->axis]>getVecEpsilon() ) {
536                                         // larger coords
537                                         curr = curr->child[1];
538                                 } else if(ray.getDirection()[curr->axis]<-getVecEpsilon() ) {
539                                         // smaller coords
540                                         curr = curr->child[0];
541                                 } else {
542                                         // paralell, order doesnt really matter are min/max/plane ok?
543                                         mpNodeStack->elem[ mpNodeStack->stackPtr ].node    = curr->child[0];
544                                         mpNodeStack->elem[ mpNodeStack->stackPtr ].mindist = planedist;
545                                         mpNodeStack->elem[ mpNodeStack->stackPtr ].maxdist = maxdist;
546                                         (mpNodeStack->stackPtr)++;
547                                         curr    = curr->child[1];
548                                         maxdist = planedist;
549                                 }
550                         } else {
551                                 // normal ray
552                                 if( (planedist>maxdist) || (planedist<0.0-getVecEpsilon() ) ) {
553                                         curr = nearChild;
554                                 } else if(planedist < mindist) {
555                                         curr = farChild;
556                                 } else {
557                                         mpNodeStack->elem[ mpNodeStack->stackPtr ].node    = farChild;
558                                         mpNodeStack->elem[ mpNodeStack->stackPtr ].mindist = planedist;
559                                         mpNodeStack->elem[ mpNodeStack->stackPtr ].maxdist = maxdist;
560                                         (mpNodeStack->stackPtr)++;
561
562                                         curr    = nearChild;
563                                         maxdist = planedist;
564                                 }
565                         } 
566     }
567         
568     
569     /* intersect with current node */
570     for (vector<ntlTriangle *>::iterator iter = curr->members->begin();
571                                  iter != curr->members->end(); iter++ ) {
572
573                         /* check for triangle flags before intersecting */
574                         if((!flags) || ( ((*iter)->getFlags() & flags) > 0 )) {
575
576                                 if( ((*iter)->getLastRay() == ray.getID() )&&((*iter)->getLastRay()>0) ) {
577                                         // was already intersected...
578                                 } else {
579                                         // we still need to intersect this triangle
580                                         gfxReal u=0.0,v=0.0, t=-1.0;
581                                         ray.intersectTriangle( mpVertices, (*iter), t,u,v);
582                                         (*iter)->setLastRay( ray.getID() );
583                                         
584                                         if( (t > 0.0) && (t<mint) )  {
585                                                 mint = t;         
586                                                 hit = (*iter);
587                                                 mintu = u; mintv = v;
588                                         }
589                                 }
590
591                         } // flags check
592     }
593
594     /* check if intersection is valid */
595     if( (mint>0.0) && (mint < GFX_REAL_MAX) ) {
596       pos = ray.getOrigin() + ray.getDirection()*mint;
597
598       if( (pos[0] >= curr->min[0]) && (pos[0] <= curr->max[0]) &&
599                                         (pos[1] >= curr->min[1]) && (pos[1] <= curr->max[1]) &&
600                                         (pos[2] >= curr->min[2]) && (pos[2] <= curr->max[2]) ) 
601                         {
602
603                                 if(forceNonsmooth) {
604                                         // calculate triangle normal
605                                         ntlVec3Gfx e0,e1,e2;
606                                         e0 = (*mpVertices)[ hit->getPoints()[0] ];
607                                         e1 = (*mpVertices)[ hit->getPoints()[1] ];
608                                         e2 = (*mpVertices)[ hit->getPoints()[2] ];
609                                         retnormal = cross( -(e2-e0), (e1-e0) );
610                                 } else {
611                                         // calculate interpolated normal
612                                         retnormal = (*mpVertNormals)[ hit->getPoints()[0] ] * (1.0-mintu-mintv)+
613                                                 (*mpVertNormals)[ hit->getPoints()[1] ]*mintu +
614                                                 (*mpVertNormals)[ hit->getPoints()[2] ]*mintv;
615                                 }
616                                 normalize(retnormal);
617                                 normal = retnormal;
618                                 distance = mint;
619                                 tri = hit;
620                                 return;
621       }
622     }    
623
624     (mpNodeStack->stackPtr)--;
625     curr    = mpNodeStack->elem[ mpNodeStack->stackPtr ].node;
626     mindist = mpNodeStack->elem[ mpNodeStack->stackPtr ].mindist;
627     maxdist = mpNodeStack->elem[ mpNodeStack->stackPtr ].maxdist;
628   } /* traverse tree */
629
630         if(mint == GFX_REAL_MAX) {
631                 distance = -1.0;
632         } else {
633                 // intersection outside the BSP bounding volumes might occur due to roundoff...
634                 //retnormal = (*mpVertNormals)[ hit->getPoints()[0] ] * (1.0-mintu-mintv)+ (*mpVertNormals)[ hit->getPoints()[1] ]*mintu + (*mpVertNormals)[ hit->getPoints()[2] ]*mintv;
635                 if(forceNonsmooth) {
636                         // calculate triangle normal
637                         ntlVec3Gfx e0,e1,e2;
638                         e0 = (*mpVertices)[ hit->getPoints()[0] ];
639                         e1 = (*mpVertices)[ hit->getPoints()[1] ];
640                         e2 = (*mpVertices)[ hit->getPoints()[2] ];
641                         retnormal = cross( -(e2-e0), (e1-e0) );
642                 } else {
643                         // calculate interpolated normal
644                         retnormal = (*mpVertNormals)[ hit->getPoints()[0] ] * (1.0-mintu-mintv)+
645                                 (*mpVertNormals)[ hit->getPoints()[1] ]*mintu +
646                                 (*mpVertNormals)[ hit->getPoints()[2] ]*mintv;
647                 }
648
649                 normalize(retnormal);
650                 normal = retnormal;
651                 distance = mint;
652                 tri = hit;
653         }
654         return;
655 }
656
657 inline void ntlRay::intersectTriangleX(vector<ntlVec3Gfx> *mpV, ntlTriangle *tri, gfxReal &t, gfxReal &u, gfxReal &v) const
658 {
659   /* (cf. moeller&haines, page 305) */
660   t = GFX_REAL_MAX;
661   ntlVec3Gfx  e0 = (*mpV)[ tri->getPoints()[0] ];
662   ntlVec3Gfx  e1 = (*mpV)[ tri->getPoints()[1] ] - e0;
663   ntlVec3Gfx  e2 = (*mpV)[ tri->getPoints()[2] ] - e0;
664
665   //ntlVec3Gfx  p  = cross( mDirection, e2 );
666   //ntlVector3Dim<Scalar> cp( (-), (- (1.0 *v[2])), ((1.0 *v[1]) -) );
667   ntlVec3Gfx  p(0.0, -e2[2], e2[1]);
668
669   gfxReal divisor  = dot(e1, p);        
670   if((divisor > -RAY_TRIANGLE_EPSILON)&&(divisor < RAY_TRIANGLE_EPSILON)) return;
671       
672   gfxReal invDivisor  = 1/divisor;
673   ntlVec3Gfx  s  = mOrigin - e0;
674   u  = invDivisor * dot(s, p);
675   if( (u<0.0-RAY_TRIANGLE_EPSILON) || (u>1.0+RAY_TRIANGLE_EPSILON) ) return;
676       
677   ntlVec3Gfx  q  = cross( s,e1 );
678   //v  = invDivisor * dot(mDirection, q);
679   v  = invDivisor * q[0];
680   if( (v<0.0-RAY_TRIANGLE_EPSILON) || ((u+v)>1.0+RAY_TRIANGLE_EPSILON) ) return;
681       
682   t = invDivisor * dot(e2, q);
683 }
684 void ntlTree::intersectX(const ntlRay &ray, gfxReal &distance, 
685                 ntlVec3Gfx &normal, 
686                 ntlTriangle *&tri, 
687                 int flags, bool forceNonsmooth) const
688 {
689   gfxReal mint = GFX_REAL_MAX;  /* current minimal t */
690   ntlVec3Gfx  retnormal;       /* intersection (interpolated) normal */
691         gfxReal mintu=0.0, mintv=0.0;    /* u,v for min t intersection */
692
693   BSPNode *curr, *nearChild, *farChild; /* current node and children */
694   gfxReal  planedist, mindist, maxdist;
695   ntlVec3Gfx   pos;
696
697         ntlTriangle *hit = NULL;
698         tri = NULL;
699
700   ray.intersectCompleteAABB(mStart,mEnd,mindist,maxdist); // +X
701
702   if((maxdist < 0.0) ||
703                  (!mpRoot) ||
704      (mindist == GFX_REAL_MAX) ||
705      (maxdist == GFX_REAL_MAX) ) {
706     distance = -1.0;
707     return;
708   }
709   mindist -= getVecEpsilon();
710   maxdist += getVecEpsilon();
711
712   /* stack init */
713   mpNodeStack->elem[0].node = NULL;
714   mpNodeStack->stackPtr = 1;
715
716   curr = mpRoot;  
717   mint = GFX_REAL_MAX;
718   while(curr != NULL) { // +X
719
720     while( !curr->isLeaf() ) {
721       planedist = distanceToPlane(curr, curr->child[0]->max, ray );
722       getChildren(curr, ray.getOrigin(), nearChild, farChild );
723
724                         // check ray direction for small plane distances
725       if( (planedist>-getVecEpsilon() )&&(planedist< getVecEpsilon() ) ) {
726                                 // ray origin on intersection plane
727                                 planedist = 0.0;
728                                 if(ray.getDirection()[curr->axis]>getVecEpsilon() ) {
729                                         // larger coords
730                                         curr = curr->child[1];
731                                 } else if(ray.getDirection()[curr->axis]<-getVecEpsilon() ) {
732                                         // smaller coords
733                                         curr = curr->child[0];
734                                 } else {
735                                         // paralell, order doesnt really matter are min/max/plane ok?
736                                         mpNodeStack->elem[ mpNodeStack->stackPtr ].node    = curr->child[0];
737                                         mpNodeStack->elem[ mpNodeStack->stackPtr ].mindist = planedist;
738                                         mpNodeStack->elem[ mpNodeStack->stackPtr ].maxdist = maxdist;
739                                         (mpNodeStack->stackPtr)++;
740                                         curr    = curr->child[1];
741                                         maxdist = planedist;
742                                 }
743                         } else {
744                                 // normal ray
745                                 if( (planedist>maxdist) || (planedist<0.0-getVecEpsilon() ) ) {
746                                         curr = nearChild;
747                                 } else if(planedist < mindist) {
748                                         curr = farChild;
749                                 } else {
750                                         mpNodeStack->elem[ mpNodeStack->stackPtr ].node    = farChild;
751                                         mpNodeStack->elem[ mpNodeStack->stackPtr ].mindist = planedist;
752                                         mpNodeStack->elem[ mpNodeStack->stackPtr ].maxdist = maxdist;
753                                         (mpNodeStack->stackPtr)++;
754
755                                         curr    = nearChild;
756                                         maxdist = planedist;
757                                 }
758                         } 
759     } // +X
760         
761     
762     /* intersect with current node */
763     for (vector<ntlTriangle *>::iterator iter = curr->members->begin();
764                                  iter != curr->members->end(); iter++ ) {
765
766                         /* check for triangle flags before intersecting */
767                         if((!flags) || ( ((*iter)->getFlags() & flags) > 0 )) {
768
769                                 if( ((*iter)->getLastRay() == ray.getID() )&&((*iter)->getLastRay()>0) ) {
770                                         // was already intersected...
771                                 } else {
772                                         // we still need to intersect this triangle
773                                         gfxReal u=0.0,v=0.0, t=-1.0;
774                                         ray.intersectTriangleX( mpVertices, (*iter), t,u,v);
775                                         (*iter)->setLastRay( ray.getID() );
776                                         
777                                         if( (t > 0.0) && (t<mint) )  {
778                                                 mint = t;         
779                                                 hit = (*iter);
780                                                 mintu = u; mintv = v;
781                                         }
782                                 }
783
784                         } // flags check
785     } // +X
786
787     /* check if intersection is valid */
788     if( (mint>0.0) && (mint < GFX_REAL_MAX) ) {
789       pos = ray.getOrigin() + ray.getDirection()*mint;
790
791       if( (pos[0] >= curr->min[0]) && (pos[0] <= curr->max[0]) &&
792                                         (pos[1] >= curr->min[1]) && (pos[1] <= curr->max[1]) &&
793                                         (pos[2] >= curr->min[2]) && (pos[2] <= curr->max[2]) ) 
794                         {
795
796                                 if(forceNonsmooth) {
797                                         // calculate triangle normal
798                                         ntlVec3Gfx e0,e1,e2;
799                                         e0 = (*mpVertices)[ hit->getPoints()[0] ];
800                                         e1 = (*mpVertices)[ hit->getPoints()[1] ];
801                                         e2 = (*mpVertices)[ hit->getPoints()[2] ];
802                                         retnormal = cross( -(e2-e0), (e1-e0) );
803                                 } else {
804                                         // calculate interpolated normal
805                                         retnormal = (*mpVertNormals)[ hit->getPoints()[0] ] * (1.0-mintu-mintv)+
806                                                 (*mpVertNormals)[ hit->getPoints()[1] ]*mintu +
807                                                 (*mpVertNormals)[ hit->getPoints()[2] ]*mintv;
808                                 }
809                                 normalize(retnormal);
810                                 normal = retnormal;
811                                 distance = mint;
812                                 tri = hit;
813                                 return;
814       }
815     }     // +X
816
817     (mpNodeStack->stackPtr)--;
818     curr    = mpNodeStack->elem[ mpNodeStack->stackPtr ].node;
819     mindist = mpNodeStack->elem[ mpNodeStack->stackPtr ].mindist;
820     maxdist = mpNodeStack->elem[ mpNodeStack->stackPtr ].maxdist;
821   } /* traverse tree */
822
823         if(mint == GFX_REAL_MAX) {
824                 distance = -1.0;
825         } else {
826
827                 // intersection outside the BSP bounding volumes might occur due to roundoff...
828                 if(forceNonsmooth) {
829                         // calculate triangle normal
830                         ntlVec3Gfx e0,e1,e2;
831                         e0 = (*mpVertices)[ hit->getPoints()[0] ];
832                         e1 = (*mpVertices)[ hit->getPoints()[1] ];
833                         e2 = (*mpVertices)[ hit->getPoints()[2] ];
834                         retnormal = cross( -(e2-e0), (e1-e0) );
835                 } else {
836                         // calculate interpolated normal
837                         retnormal = (*mpVertNormals)[ hit->getPoints()[0] ] * (1.0-mintu-mintv)+
838                                 (*mpVertNormals)[ hit->getPoints()[1] ]*mintu +
839                                 (*mpVertNormals)[ hit->getPoints()[2] ]*mintv;
840                 }
841
842                 normalize(retnormal);
843                 normal = retnormal;
844                 distance = mint;
845                 tri = hit;
846         } // +X
847         return;
848 }
849
850
851
852 /******************************************************************************
853  * distance to plane function for nodes
854  *****************************************************************************/
855 gfxReal ntlTree::distanceToPlane(BSPNode *curr, ntlVec3Gfx plane, ntlRay ray) const
856 {
857   return ( (plane[curr->axis]-ray.getOrigin()[curr->axis]) / ray.getDirection()[curr->axis] );
858 }
859
860
861 /******************************************************************************
862  * return ordering of children nodes relatice to origin point
863  *****************************************************************************/
864 void ntlTree::getChildren(BSPNode *curr, ntlVec3Gfx origin, BSPNode *&node_near, BSPNode *&node_far) const 
865 {
866   if(curr->child[0]->max[ curr->axis ] >= origin[ curr->axis ]) {
867     node_near = curr->child[0];
868     node_far = curr->child[1];
869   } else {
870     node_near = curr->child[1];
871     node_far = curr->child[0];
872   }
873 }
874
875
876 /******************************************************************************
877  * delete a node of the tree with all sub nodes
878  *  dont delete root members 
879  *****************************************************************************/
880 void ntlTree::deleteNode(BSPNode *curr) 
881 {
882         if(!curr) return;
883
884   if(curr->child[0] != NULL)
885     deleteNode(curr->child[0]);
886   if(curr->child[1] != NULL)
887     deleteNode(curr->child[1]);
888
889   if(curr->members != NULL) delete curr->members;
890   delete curr;
891 }
892
893
894
895 /******************************************************************
896  * intersect only front or backsides
897  * currently unused
898  */
899 inline void ntlRay::intersectTriangleFront(vector<ntlVec3Gfx> *mpV, ntlTriangle *tri, gfxReal &t, gfxReal &u, gfxReal &v) const
900 {
901   t = GFX_REAL_MAX;
902   ntlVec3Gfx  e0 = (*mpV)[ tri->getPoints()[0] ];
903   ntlVec3Gfx  e1 = (*mpV)[ tri->getPoints()[1] ] - e0;
904   ntlVec3Gfx  e2 = (*mpV)[ tri->getPoints()[2] ] - e0;
905   ntlVec3Gfx  p  = cross( mDirection, e2 );
906   gfxReal a  = dot(e1, p);      
907   //if((a > -RAY_TRIANGLE_EPSILON)&&(a < RAY_TRIANGLE_EPSILON)) return;
908   if(a < RAY_TRIANGLE_EPSILON) return; // cull backsides
909       
910   gfxReal f  = 1/a;
911   ntlVec3Gfx  s  = mOrigin - e0;
912   u  = f * dot(s, p);
913   if( (u<0.0-RAY_TRIANGLE_EPSILON) || (u>1.0+RAY_TRIANGLE_EPSILON) ) return;
914       
915   ntlVec3Gfx  q  = cross( s,e1 );
916   v  = f * dot(mDirection, q);
917   if( (v<0.0-RAY_TRIANGLE_EPSILON) || ((u+v)>1.0+RAY_TRIANGLE_EPSILON) ) return;
918       
919   t = f * dot(e2, q);      
920 }
921 inline void ntlRay::intersectTriangleBack(vector<ntlVec3Gfx> *mpV, ntlTriangle *tri, gfxReal &t, gfxReal &u, gfxReal &v) const
922 {
923   t = GFX_REAL_MAX;
924   ntlVec3Gfx  e0 = (*mpV)[ tri->getPoints()[0] ];
925   ntlVec3Gfx  e1 = (*mpV)[ tri->getPoints()[1] ] - e0;
926   ntlVec3Gfx  e2 = (*mpV)[ tri->getPoints()[2] ] - e0;
927   ntlVec3Gfx  p  = cross( mDirection, e2 );
928   gfxReal a  = dot(e1, p);      
929   //if((a > -RAY_TRIANGLE_EPSILON)&&(a < RAY_TRIANGLE_EPSILON)) return;
930   if(a > -RAY_TRIANGLE_EPSILON) return; // cull frontsides
931       
932   gfxReal f  = 1/a;
933   ntlVec3Gfx  s  = mOrigin - e0;
934   u  = f * dot(s, p);
935   if( (u<0.0-RAY_TRIANGLE_EPSILON) || (u>1.0+RAY_TRIANGLE_EPSILON) ) return;
936       
937   ntlVec3Gfx  q  = cross( s,e1 );
938   v  = f * dot(mDirection, q);
939   if( (v<0.0-RAY_TRIANGLE_EPSILON) || ((u+v)>1.0+RAY_TRIANGLE_EPSILON) ) return;
940       
941   t = f * dot(e2, q);      
942 }
943
944
945