source: trunk/VUT/GtpVisibilityPreprocessor/src/KdTree.cpp @ 241

Revision 241, 17.3 KB checked in by mattausch, 19 years ago (diff)

added bsp stuff

Line 
1#include <stack>
2#include <algorithm>
3#include <queue>
4#include "Environment.h"
5#include "Mesh.h"
6#include "KdTree.h"
7
8
9int KdNode::mailID = 1;
10
11
12KdNode::KdNode(KdInterior *parent):mParent(parent), mailbox(0)
13{
14  if (parent)
15    mDepth = parent->mDepth+1;
16  else
17    mDepth = 0;
18}
19
20KdTree::KdTree()
21{
22  mRoot = new KdLeaf(NULL, 0);
23  environment->GetIntValue("KdTree.Termination.maxDepth", mTermMaxDepth);
24  environment->GetIntValue("KdTree.Termination.minCost", mTermMinCost);
25  environment->GetFloatValue("KdTree.Termination.maxCostRatio", mMaxCostRatio);
26  environment->GetFloatValue("KdTree.Termination.ct_div_ci", mCt_div_ci);
27  environment->GetFloatValue("KdTree.splitBorder", mSplitBorder);
28
29  environment->GetBoolValue("KdTree.sahUseFaces", mSahUseFaces);
30
31  char splitType[64];
32  environment->GetStringValue("KdTree.splitMethod", splitType);
33 
34  mSplitMethod = SPLIT_SPATIAL_MEDIAN;
35  if (strcmp(splitType, "spatialMedian") == 0)
36    mSplitMethod = SPLIT_SPATIAL_MEDIAN;
37  else
38    if (strcmp(splitType, "objectMedian") == 0)
39      mSplitMethod = SPLIT_OBJECT_MEDIAN;
40    else
41      if (strcmp(splitType, "SAH") == 0)
42        mSplitMethod = SPLIT_SAH;
43      else {
44        cerr<<"Wrong kd split type "<<splitType<<endl;
45        exit(1);
46      }
47  splitCandidates = NULL;
48}
49
50bool
51KdTree::Construct()
52{
53  if (!splitCandidates)
54    splitCandidates = new vector<SortableEntry>;
55
56  // first construct a leaf that will get subdivide
57  KdLeaf *leaf = (KdLeaf *) mRoot;
58
59  mStat.nodes = 1;
60
61  mBox.Initialize();
62 
63  ObjectContainer::const_iterator mi;
64  for ( mi = leaf->mObjects.begin();
65        mi != leaf->mObjects.end();
66        mi++) {
67    mBox.Include((*mi)->GetBox());
68  }
69
70  cout <<"KdTree Root Box:"<< mBox<<endl;
71  mRoot = Subdivide(TraversalData(leaf, mBox, 0));
72
73  // remove the allocated array
74  delete splitCandidates;
75
76  return true;
77}
78
79KdNode *
80KdTree::Subdivide(const TraversalData &tdata)
81{
82
83  KdNode *result = NULL;
84
85  priority_queue<TraversalData> tStack;
86  //  stack<STraversalData> tStack;
87 
88  tStack.push(tdata);
89  AxisAlignedBox3 backBox, frontBox;
90
91 
92  while (!tStack.empty()) {
93
94#if 0
95    if ( GetMemUsage() > maxMemory ) {
96      // count statistics on unprocessed leafs
97      while (!tStack.empty()) {
98        EvaluateLeafStats(tStack.top());
99        tStack.pop();
100      }
101      break;
102    }
103#endif
104
105    TraversalData data = tStack.top();
106    tStack.pop();
107   
108    KdNode *node = SubdivideNode((KdLeaf *) data.mNode,
109                                 data.mBox,
110                                 backBox,
111                                 frontBox
112                                 );
113    if (result == NULL)
114      result = node;
115   
116    if (!node->IsLeaf()) {
117
118      KdInterior *interior = (KdInterior *) node;
119      // push the children on the stack
120      tStack.push(TraversalData(interior->mBack, backBox, data.mDepth+1));
121      tStack.push(TraversalData(interior->mFront, frontBox, data.mDepth+1));
122     
123    } else {
124      EvaluateLeafStats(data);
125    }
126  }
127 
128  return result;
129
130}
131
132
133
134bool
135KdTree::TerminationCriteriaMet(const KdLeaf *leaf)
136{
137  //  cerr<<"\n OBJECTS="<<leaf->mObjects.size()<<endl;
138  return
139    (leaf->mObjects.size() <= mTermMinCost) ||
140    (leaf->mDepth >= mTermMaxDepth);
141 
142}
143
144
145int
146KdTree::SelectPlane(KdLeaf *leaf,
147                    const AxisAlignedBox3 &box,
148                    float &position
149                    )
150{
151  int axis = -1;
152 
153  switch (mSplitMethod)
154    {
155    case SPLIT_SPATIAL_MEDIAN: {
156      axis = box.Size().DrivingAxis();
157      position = (box.Min()[axis] + box.Max()[axis])*0.5f;
158      break;
159    }
160    case SPLIT_SAH: {
161      int objectsBack, objectsFront;
162      float costRatio;
163      bool mOnlyDrivingAxis = false;
164      if (mOnlyDrivingAxis) {
165        axis = box.Size().DrivingAxis();
166        costRatio = BestCostRatio(leaf,
167                                  box,
168                                  axis,
169                                  position,
170                                  objectsBack,
171                                  objectsFront);
172      } else {
173        costRatio = MAX_FLOAT;
174        for (int i=0; i < 3; i++) {
175          float p;
176          float r = BestCostRatio(leaf,
177                                  box,
178                                  i,
179                                  p,
180                                  objectsBack,
181                                  objectsFront);
182          if (r < costRatio) {
183            costRatio = r;
184            axis = i;
185            position = p;
186          }
187        }
188      }
189     
190      if (costRatio > mMaxCostRatio) {
191        //      cout<<"Too big cost ratio "<<costRatio<<endl;
192        axis = -1;
193      }
194      break;
195    }
196   
197    }
198  return axis;
199}
200
201KdNode *
202KdTree::SubdivideNode(
203                      KdLeaf *leaf,
204                      const AxisAlignedBox3 &box,
205                      AxisAlignedBox3 &backBBox,
206                      AxisAlignedBox3 &frontBBox
207                      )
208{
209 
210  if (TerminationCriteriaMet(leaf))
211    return leaf;
212 
213  float position;
214 
215  // select subdivision axis
216  int axis = SelectPlane( leaf, box, position );
217
218  if (axis == -1) {
219    return leaf;
220  }
221 
222  mStat.nodes+=2;
223  mStat.splits[axis]++;
224 
225  // add the new nodes to the tree
226  KdInterior *node = new KdInterior(leaf->mParent);
227
228  node->mAxis = axis;
229  node->mPosition = position;
230  node->mBox = box;
231 
232  backBBox = box;
233  frontBBox = box;
234 
235  // first count ray sides
236  int objectsBack = 0;
237  int objectsFront = 0;
238 
239  backBBox.SetMax(axis, position);
240  frontBBox.SetMin(axis, position);
241
242  ObjectContainer::const_iterator mi;
243
244  for ( mi = leaf->mObjects.begin();
245        mi != leaf->mObjects.end();
246        mi++) {
247    // determine the side of this ray with respect to the plane
248    AxisAlignedBox3 box = (*mi)->GetBox();
249    if (box.Max(axis) > position )
250      objectsFront++;
251   
252    if (box.Min(axis) < position )
253      objectsBack++;
254  }
255
256 
257  KdLeaf *back = new KdLeaf(node, objectsBack);
258  KdLeaf *front = new KdLeaf(node, objectsFront);
259
260  // replace a link from node's parent
261  if (  leaf->mParent )
262    leaf->mParent->ReplaceChildLink(leaf, node);
263
264  // and setup child links
265  node->SetupChildLinks(back, front);
266 
267  for (mi = leaf->mObjects.begin();
268       mi != leaf->mObjects.end();
269       mi++) {
270    // determine the side of this ray with respect to the plane
271    AxisAlignedBox3 box = (*mi)->GetBox();
272
273    if (box.Max(axis) >= position )
274      front->mObjects.push_back(*mi);
275   
276    if (box.Min(axis) < position )
277      back->mObjects.push_back(*mi);
278   
279    mStat.objectRefs -= leaf->mObjects.size();
280    mStat.objectRefs += objectsBack + objectsFront;
281  }
282 
283  delete leaf;
284  return node;
285}
286
287
288
289void
290KdTreeStatistics::Print(ostream &app) const
291{
292  app << "===== KdTree statistics ===============\n";
293
294  app << "#N_RAYS Number of rays )\n"
295      << rays <<endl;
296  app << "#N_DOMAINS  ( Number of query domains )\n"
297      << queryDomains <<endl;
298 
299  app << "#N_NODES ( Number of nodes )\n" << nodes << "\n";
300
301  app << "#N_LEAVES ( Number of leaves )\n" << Leaves() << "\n";
302
303  app << "#N_SPLITS ( Number of splits in axes x y z dx dy dz \n";
304  for (int i=0; i<7; i++)
305    app << splits[i] <<" ";
306  app <<endl;
307
308  app << "#N_RAYREFS ( Number of rayRefs )\n" <<
309    rayRefs << "\n";
310
311  app << "#N_RAYRAYREFS  ( Number of rayRefs / ray )\n" <<
312    rayRefs/(double)rays << "\n";
313
314  app << "#N_LEAFRAYREFS  ( Number of rayRefs / leaf )\n" <<
315    rayRefs/(double)Leaves() << "\n";
316
317  app << "#N_MAXOBJECTREFS  ( Max number of rayRefs / leaf )\n" <<
318    maxObjectRefs << "\n";
319
320  app << "#N_NONEMPTYRAYREFS  ( Number of rayRefs in nonEmpty leaves / non empty leaf )\n" <<
321    rayRefsNonZeroQuery/(double)(Leaves() - zeroQueryNodes) << "\n";
322
323  app << "#N_LEAFDOMAINREFS  ( Number of query domain Refs / leaf )\n" <<
324    objectRefs/(double)Leaves() << "\n";
325
326  //  app << setprecision(4);
327
328  app << "#N_PEMPTYLEAVES  ( Percentage of leaves with zero query domains )\n"<<
329    zeroQueryNodes*100/(double)Leaves()<<endl;
330
331  app << "#N_PMAXDEPTHLEAVES ( Percentage of leaves at maxdepth )\n"<<
332    maxDepthNodes*100/(double)Leaves()<<endl;
333
334  app << "#N_PMINCOSTLEAVES  ( Percentage of leaves with minCost )\n"<<
335    minCostNodes*100/(double)Leaves()<<endl;
336
337  app << "#N_ADDED_RAYREFS  (Number of dynamically added ray references )\n"<<
338    addedRayRefs<<endl;
339
340  app << "#N_REMOVED_RAYREFS  (Number of dynamically removed ray references )\n"<<
341    removedRayRefs<<endl;
342
343  //  app << setprecision(4);
344
345  //  app << "#N_CTIME  ( Construction time [s] )\n"
346  //      << Time() << " \n";
347
348  app << "===== END OF KdTree statistics ==========\n";
349
350}
351
352
353
354void
355KdTree::EvaluateLeafStats(const TraversalData &data)
356{
357
358  // the node became a leaf -> evaluate stats for leafs
359  KdLeaf *leaf = (KdLeaf *)data.mNode;
360
361  if (data.mDepth > mTermMaxDepth)
362    mStat.maxDepthNodes++;
363 
364  if ( (int)(leaf->mObjects.size()) < mTermMinCost)
365    mStat.minCostNodes++;
366 
367 
368  if ( (int)(leaf->mObjects.size()) > mStat.maxObjectRefs)
369    mStat.maxObjectRefs = leaf->mObjects.size();
370 
371}
372
373
374
375void
376KdTree::SortSplitCandidates(
377                            KdLeaf *node,
378                            const int axis
379                            )
380{
381  splitCandidates->clear();
382 
383  int requestedSize = 2*node->mObjects.size();
384  // creates a sorted split candidates array
385  if (splitCandidates->capacity() > 500000 &&
386      requestedSize < (int)(splitCandidates->capacity()/10) ) {
387    delete splitCandidates;
388    splitCandidates = new vector<SortableEntry>;
389  }
390 
391  splitCandidates->reserve(requestedSize);
392 
393  // insert all queries
394  for(ObjectContainer::const_iterator mi = node->mObjects.begin();
395      mi != node->mObjects.end();
396      mi++) {
397    AxisAlignedBox3 box = (*mi)->GetBox();
398
399    splitCandidates->push_back(SortableEntry(SortableEntry::BOX_MIN,
400                                             box.Min(axis),
401                                             *mi)
402                               );
403   
404   
405    splitCandidates->push_back(SortableEntry(SortableEntry::BOX_MAX,
406                                             box.Max(axis),
407                                             *mi)
408                               );
409  }
410 
411  stable_sort(splitCandidates->begin(), splitCandidates->end());
412}
413
414
415float
416KdTree::BestCostRatio(
417                      KdLeaf *node,
418                      const AxisAlignedBox3 &box,
419                      const int axis,
420                      float &position,
421                      int &objectsBack,
422                      int &objectsFront
423                      )
424{
425
426  SortSplitCandidates(node, axis);
427 
428  // go through the lists, count the number of objects left and right
429  // and evaluate the following cost funcion:
430  // C = ct_div_ci  + (ol + or)/queries
431
432  float totalIntersections = 0.0f;
433  vector<SortableEntry>::const_iterator ci;
434
435  for(ci = splitCandidates->begin();
436      ci < splitCandidates->end();
437      ci++)
438    if ((*ci).type == SortableEntry::BOX_MIN) {
439      totalIntersections += (*ci).intersectable->IntersectionComplexity();
440    }
441
442  float intersectionsLeft = 0;
443  float intersectionsRight = totalIntersections;
444
445  int objectsLeft = 0, objectsRight = node->mObjects.size();
446 
447  float minBox = box.Min(axis);
448  float maxBox = box.Max(axis);
449  float boxArea = box.SurfaceArea();
450 
451  float minBand = minBox + mSplitBorder*(maxBox - minBox);
452  float maxBand = minBox + (1.0f - mSplitBorder)*(maxBox - minBox);
453 
454  float minSum = 1e20;
455 
456  for(ci = splitCandidates->begin();
457      ci < splitCandidates->end();
458      ci++) {
459    switch ((*ci).type) {
460    case SortableEntry::BOX_MIN:
461      objectsLeft++;
462      intersectionsLeft += (*ci).intersectable->IntersectionComplexity();
463      break;
464    case SortableEntry::BOX_MAX:
465      objectsRight--;
466      intersectionsRight -= (*ci).intersectable->IntersectionComplexity();
467      break;
468    }
469
470    if ((*ci).value > minBand && (*ci).value < maxBand) {
471      AxisAlignedBox3 lbox = box;
472      AxisAlignedBox3 rbox = box;
473      lbox.SetMax(axis, (*ci).value);
474      rbox.SetMin(axis, (*ci).value);
475
476      float sum;
477      if (mSahUseFaces)
478        sum = intersectionsLeft*lbox.SurfaceArea() + intersectionsRight*rbox.SurfaceArea();
479      else
480        sum = objectsLeft*lbox.SurfaceArea() + objectsRight*rbox.SurfaceArea();
481     
482      //      cout<<"pos="<<(*ci).value<<"\t q=("<<ql<<","<<qr<<")\t r=("<<rl<<","<<rr<<")"<<endl;
483      //      cout<<"cost= "<<sum<<endl;
484     
485      if (sum < minSum) {
486        minSum = sum;
487        position = (*ci).value;
488       
489        objectsBack = objectsLeft;
490        objectsFront = objectsRight;
491      }
492    }
493  }
494 
495  float oldCost = mSahUseFaces ? totalIntersections : node->mObjects.size();
496  float newCost = mCt_div_ci + minSum/boxArea;
497  float ratio = newCost/oldCost;
498 
499#if 0
500  cout<<"===================="<<endl;
501  cout<<"costRatio="<<ratio<<" pos="<<position<<" t="<<(position - minBox)/(maxBox - minBox)
502      <<"\t o=("<<objectsBack<<","<<objectsFront<<")"<<endl;
503#endif
504  return ratio;
505}
506
507int
508KdTree::CastRay(
509                Ray &ray
510                )
511{
512  int hits = 0;
513 
514  stack<RayTraversalData> tStack;
515 
516  float maxt = 1e6;
517  float mint = 0;
518
519  Intersectable::NewMail();
520
521  if (!mBox.GetMinMaxT(ray, &mint, &maxt))
522    return 0;
523
524  if (mint < 0)
525    mint = 0;
526 
527  maxt += Limits::Threshold;
528 
529  Vector3 entp = ray.Extrap(mint);
530  Vector3 extp = ray.Extrap(maxt);
531 
532  KdNode *node = mRoot;
533  KdNode *farChild;
534  float position;
535  int axis;
536 
537  while (1) {
538    if (!node->IsLeaf()) {
539      KdInterior *in = (KdInterior *) node;
540      position = in->mPosition;
541      axis = in->mAxis;
542
543      if (entp[axis] <= position) {
544        if (extp[axis] <= position) {
545          node = in->mBack;
546          // cases N1,N2,N3,P5,Z2,Z3
547          continue;
548        } else {
549          // case N4
550          node = in->mBack;
551          farChild = in->mFront;
552        }
553      }
554      else {
555        if (position <= extp[axis]) {
556          node = in->mFront;
557          // cases P1,P2,P3,N5,Z1
558          continue;
559        } else {
560          node = in->mFront;
561          farChild = in->mBack;
562          // case P4
563        }
564          }
565      // $$ modification 3.5.2004 - hints from Kamil Ghais
566      // case N4 or P4
567      float tdist = (position - ray.GetLoc(axis)) / ray.GetDir(axis);
568      tStack.push(RayTraversalData(farChild, extp, maxt));
569      extp = ray.GetLoc() + ray.GetDir()*tdist;
570      maxt = tdist;
571        } else {
572      // compute intersection with all objects in this leaf
573      KdLeaf *leaf = (KdLeaf *) node;
574      ray.leaves.push_back(leaf);
575     
576      ObjectContainer::const_iterator mi;
577      for ( mi = leaf->mObjects.begin();
578            mi != leaf->mObjects.end();
579            mi++) {
580        Intersectable *object = *mi;
581        if (!object->Mailed() ) {
582          object->Mail();
583          //ray.meshes.push_back(mesh);
584          hits += object->CastRay(ray);
585        }
586      }
587
588      if (hits && ray.GetType() == Ray::LOCAL_RAY)
589        if (ray.intersections[0].mT <= maxt)
590          break;
591     
592      // get the next node from the stack
593      if (tStack.empty())
594        break;
595     
596      entp = extp;
597      mint = maxt;
598      RayTraversalData &s  = tStack.top();
599      node = s.mNode;
600      extp = s.mExitPoint;
601      maxt = s.mMaxT;
602      tStack.pop();
603    }
604  }
605 
606 
607  return hits;
608}
609
610void
611KdTree::CollectObjects(KdNode *n, ObjectContainer &objects)
612{
613  stack<KdNode *> nodeStack;
614
615  nodeStack.push(n);
616
617  while (!nodeStack.empty()) {
618    KdNode *node = nodeStack.top();
619    nodeStack.pop();
620    if (node->IsLeaf()) {
621      KdLeaf *leaf = (KdLeaf *)node;
622      for (int j=0; j < leaf->mObjects.size(); j++) {
623        Intersectable *object = leaf->mObjects[j];
624        if (!object->Mailed()) {
625          object->Mail();
626          objects.push_back(object);
627        }
628      }
629    } else {
630      KdInterior *interior = (KdInterior *)node;
631      nodeStack.push(interior->mFront);
632      nodeStack.push(interior->mBack);
633    }
634  }
635}
636
637// Find random neighbor which was not mailed
638KdNode *
639KdTree::FindRandomNeighbor(KdNode *n,
640                           bool onlyUnmailed
641                           )
642{
643  stack<KdNode *> nodeStack;
644 
645  nodeStack.push(mRoot);
646
647  AxisAlignedBox3 box = GetBox(n);
648  int mask = rand();
649
650  while (!nodeStack.empty()) {
651    KdNode *node = nodeStack.top();
652    nodeStack.pop();
653    if (node->IsLeaf()) {
654      if ( node != n && (!onlyUnmailed || !node->Mailed()) )
655        return node;
656    } else {
657      KdInterior *interior = (KdInterior *)node;
658      if (interior->mPosition > box.Max(interior->mAxis))
659        nodeStack.push(interior->mBack);
660      else
661        if (interior->mPosition < box.Min(interior->mAxis))
662          nodeStack.push(interior->mFront);
663        else {
664          // random decision
665          if (mask&1)
666            nodeStack.push(interior->mBack);
667          else
668            nodeStack.push(interior->mFront);
669          mask = mask>>1;
670        }
671    }
672  }
673 
674  return NULL;
675}
676
677int
678KdTree::FindNeighbors(KdNode *n,
679                      vector<KdNode *> &neighbors,
680                      bool onlyUnmailed
681                      )
682{
683  stack<KdNode *> nodeStack;
684 
685  nodeStack.push(mRoot);
686
687  AxisAlignedBox3 box = GetBox(n);
688
689  while (!nodeStack.empty()) {
690    KdNode *node = nodeStack.top();
691    nodeStack.pop();
692    if (node->IsLeaf()) {
693      if ( node != n && (!onlyUnmailed || !node->Mailed()) )
694        neighbors.push_back(node);
695    } else {
696      KdInterior *interior = (KdInterior *)node;
697      if (interior->mPosition > box.Max(interior->mAxis))
698        nodeStack.push(interior->mBack);
699      else
700        if (interior->mPosition < box.Min(interior->mAxis))
701          nodeStack.push(interior->mFront);
702        else {
703          // random decision
704            nodeStack.push(interior->mBack);
705            nodeStack.push(interior->mFront);
706        }
707    }
708  }
709 
710  return neighbors.size();
711}
712
713void
714KdTree::CollectLeaves(vector<KdLeaf *> &leaves)
715{
716  stack<KdNode *> nodeStack;
717  nodeStack.push(mRoot);
718
719  while (!nodeStack.empty()) {
720    KdNode *node = nodeStack.top();
721    nodeStack.pop();
722    if (node->IsLeaf()) {
723      KdLeaf *leaf = (KdLeaf *)node;
724      leaves.push_back(leaf);
725    } else {
726      KdInterior *interior = (KdInterior *)node;
727      nodeStack.push(interior->mBack);
728      nodeStack.push(interior->mFront);
729    }
730  }
731}
732
733
734int
735KdTree::CollectLeafPvs()
736{
737  int totalPvsSize = 0;
738  stack<KdNode *> nodeStack;
739 
740  nodeStack.push(mRoot);
741 
742  while (!nodeStack.empty()) {
743    KdNode *node = nodeStack.top();
744    nodeStack.pop();
745    if (node->IsLeaf()) {
746      KdLeaf *leaf = (KdLeaf *)node;
747      for (int j=0; j < leaf->mObjects.size(); j++) {
748        Intersectable *object = leaf->mObjects[j];
749        if (!object->Mailed()) {
750          object->Mail();
751          // add this node to pvs of all nodes it can see
752          KdPvsMap::iterator ni = object->mKdPvs.mEntries.begin();
753          for (; ni != object->mKdPvs.mEntries.end(); ni++) {
754            KdNode *node = (*ni).first;
755            // $$ JB TEMPORARY solution -> should add object PVS or explictly computed
756            // kd tree PVS
757            if (leaf->mKdPvs.AddNodeSample(node))
758              totalPvsSize++;
759          }
760        }
761      }
762    } else {
763      KdInterior *interior = (KdInterior *)node;
764      nodeStack.push(interior->mFront);
765      nodeStack.push(interior->mBack);
766    }
767  }
768
769  return totalPvsSize;
770}
Note: See TracBrowser for help on using the repository browser.