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