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