source: GTP/trunk/Lib/Vis/Preprocessing/src/KdTree.cpp @ 1004

Revision 1004, 24.8 KB checked in by mattausch, 18 years ago (diff)

environment as a singleton
Line 
1#include <stack>
2#include <algorithm>
3#include <queue>
4#include "Environment.h"
5#include "Mesh.h"
6#include "KdTree.h"
7#include "ViewCell.h"
8#include "Beam.h"
9
10
11
12namespace GtpVisibilityPreprocessor {
13
14int KdNode::mailID = 1;
15
16KdNode::KdNode(KdInterior *parent):mParent(parent), mailbox(0)
17{
18  if (parent)
19    mDepth = parent->mDepth+1;
20  else
21    mDepth = 0;
22}
23
24
25KdInterior::~KdInterior()
26{
27        // recursivly destroy children
28        DEL_PTR(mFront);
29        DEL_PTR(mBack);
30}
31
32
33KdTree::KdTree()
34{
35
36 
37  mRoot = new KdLeaf(NULL, 0);
38  Environment::GetSingleton()->GetIntValue("KdTree.Termination.maxNodes", mTermMaxNodes);
39  Environment::GetSingleton()->GetIntValue("KdTree.Termination.maxDepth", mTermMaxDepth);
40  Environment::GetSingleton()->GetIntValue("KdTree.Termination.minCost", mTermMinCost);
41  Environment::GetSingleton()->GetFloatValue("KdTree.Termination.maxCostRatio", mMaxCostRatio);
42  Environment::GetSingleton()->GetFloatValue("KdTree.Termination.ct_div_ci", mCt_div_ci);
43  Environment::GetSingleton()->GetFloatValue("KdTree.splitBorder", mSplitBorder);
44
45  Environment::GetSingleton()->GetBoolValue("KdTree.sahUseFaces", mSahUseFaces);
46
47  char splitType[64];
48  Environment::GetSingleton()->GetStringValue("KdTree.splitMethod", splitType);
49 
50  mSplitMethod = SPLIT_SPATIAL_MEDIAN;
51  if (strcmp(splitType, "spatialMedian") == 0)
52    mSplitMethod = SPLIT_SPATIAL_MEDIAN;
53  else
54    if (strcmp(splitType, "objectMedian") == 0)
55      mSplitMethod = SPLIT_OBJECT_MEDIAN;
56    else
57      if (strcmp(splitType, "SAH") == 0)
58        mSplitMethod = SPLIT_SAH;
59      else {
60        cerr<<"Wrong kd split type "<<splitType<<endl;
61        exit(1);
62      }
63  splitCandidates = NULL;
64}
65
66
67KdTree::~KdTree()
68{
69        DEL_PTR(mRoot);
70}
71
72
73bool
74KdTree::Construct()
75{
76  if (!splitCandidates)
77    splitCandidates = new vector<SortableEntry>;
78
79  // first construct a leaf that will get subdivide
80  KdLeaf *leaf = (KdLeaf *) mRoot;
81
82  mStat.nodes = 1;
83
84  mBox.Initialize();
85  ObjectContainer::const_iterator mi;
86  for ( mi = leaf->mObjects.begin();
87        mi != leaf->mObjects.end();
88        mi++) {
89        //      cout<<(*mi)->GetBox()<<endl;
90        mBox.Include((*mi)->GetBox());
91  }
92
93  cout <<"KdTree Root Box:"<<mBox<<endl;
94  mRoot = Subdivide(TraversalData(leaf, mBox, 0));
95
96  // remove the allocated array
97  delete splitCandidates;
98
99  return true;
100}
101
102KdNode *
103KdTree::Subdivide(const TraversalData &tdata)
104{
105
106  KdNode *result = NULL;
107
108  priority_queue<TraversalData> tStack;
109  //  stack<STraversalData> tStack;
110 
111  tStack.push(tdata);
112  AxisAlignedBox3 backBox, frontBox;
113
114  while (!tStack.empty()) {
115        //      cout<<mStat.Nodes()<<" "<<mTermMaxNodes<<endl;
116        if (mStat.Nodes() > mTermMaxNodes) {
117          //    if ( GetMemUsage() > maxMemory ) {
118      // count statistics on unprocessed leafs
119      while (!tStack.empty()) {
120                EvaluateLeafStats(tStack.top());
121                tStack.pop();
122      }
123      break;
124    }
125         
126       
127    TraversalData data = tStack.top();
128    tStack.pop();
129   
130    KdNode *node = SubdivideNode((KdLeaf *) data.mNode,
131                                                                 data.mBox,
132                                                                 backBox,
133                                                                 frontBox
134                                                                 );
135
136        if (result == NULL)
137      result = node;
138   
139    if (!node->IsLeaf()) {
140
141      KdInterior *interior = (KdInterior *) node;
142      // push the children on the stack
143      tStack.push(TraversalData(interior->mBack, backBox, data.mDepth+1));
144      tStack.push(TraversalData(interior->mFront, frontBox, data.mDepth+1));
145     
146    } else {
147      EvaluateLeafStats(data);
148    }
149  }
150 
151  return result;
152
153}
154
155
156
157bool
158KdTree::TerminationCriteriaMet(const KdLeaf *leaf)
159{
160  //  cerr<<"\n OBJECTS="<<leaf->mObjects.size()<<endl;
161  return
162    (leaf->mObjects.size() <= mTermMinCost) ||
163    (leaf->mDepth >= mTermMaxDepth);
164 
165}
166
167
168int
169KdTree::SelectPlane(KdLeaf *leaf,
170                                        const AxisAlignedBox3 &box,
171                                        float &position
172                                        )
173{
174  int axis = -1;
175 
176  switch (mSplitMethod)
177    {
178    case SPLIT_SPATIAL_MEDIAN: {
179      axis = box.Size().DrivingAxis();
180      position = (box.Min()[axis] + box.Max()[axis])*0.5f;
181      break;
182    }
183    case SPLIT_SAH: {
184      int objectsBack, objectsFront;
185      float costRatio;
186      bool mOnlyDrivingAxis = false;
187      if (mOnlyDrivingAxis) {
188                axis = box.Size().DrivingAxis();
189                costRatio = BestCostRatio(leaf,
190                                                                  box,
191                                                                  axis,
192                                                                  position,
193                                                                  objectsBack,
194                                                                  objectsFront);
195      } else {
196                costRatio = MAX_FLOAT;
197                for (int i=0; i < 3; i++) {
198                  float p;
199                  float r = BestCostRatio(leaf,
200                                                                  box,
201                                                                  i,
202                                                                  p,
203                                                                  objectsBack,
204                                                                  objectsFront);
205                  if (r < costRatio) {
206                        costRatio = r;
207                        axis = i;
208                        position = p;
209                  }
210                }
211      }
212     
213      if (costRatio > mMaxCostRatio) {
214                //cout<<"Too big cost ratio "<<costRatio<<endl;
215                axis = -1;
216      }
217      break;
218    }
219   
220    }
221  return axis;
222}
223
224KdNode *
225KdTree::SubdivideNode(
226                      KdLeaf *leaf,
227                      const AxisAlignedBox3 &box,
228                      AxisAlignedBox3 &backBBox,
229                      AxisAlignedBox3 &frontBBox
230                      )
231{
232 
233  if (TerminationCriteriaMet(leaf))
234          return leaf;
235
236  float position;
237 
238  // select subdivision axis
239  int axis = SelectPlane( leaf, box, position );
240
241  if (axis == -1) {
242    return leaf;
243  }
244 
245  mStat.nodes+=2;
246  mStat.splits[axis]++;
247 
248  // add the new nodes to the tree
249  KdInterior *node = new KdInterior(leaf->mParent);
250
251  node->mAxis = axis;
252  node->mPosition = position;
253  node->mBox = box;
254 
255  backBBox = box;
256  frontBBox = box;
257 
258  // first count ray sides
259  int objectsBack = 0;
260  int objectsFront = 0;
261 
262  backBBox.SetMax(axis, position);
263  frontBBox.SetMin(axis, position);
264
265  ObjectContainer::const_iterator mi;
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    if (box.Max(axis) > position )
273      objectsFront++;
274   
275    if (box.Min(axis) < position )
276      objectsBack++;
277  }
278
279 
280  KdLeaf *back = new KdLeaf(node, objectsBack);
281  KdLeaf *front = new KdLeaf(node, objectsFront);
282
283  // replace a link from node's parent
284  if (  leaf->mParent )
285    leaf->mParent->ReplaceChildLink(leaf, node);
286
287  // and setup child links
288  node->SetupChildLinks(back, front);
289 
290  for (mi = leaf->mObjects.begin();
291       mi != leaf->mObjects.end();
292       mi++) {
293    // determine the side of this ray with respect to the plane
294    AxisAlignedBox3 box = (*mi)->GetBox();
295
296    if (box.Max(axis) >= position )
297      front->mObjects.push_back(*mi);
298   
299    if (box.Min(axis) < position )
300      back->mObjects.push_back(*mi);
301   
302    mStat.objectRefs -= (int)leaf->mObjects.size();
303    mStat.objectRefs += objectsBack + objectsFront;
304  }
305 
306  delete leaf;
307  return node;
308}
309
310
311
312void
313KdTreeStatistics::Print(ostream &app) const
314{
315  app << "===== KdTree statistics ===============\n";
316
317  app << "#N_NODES ( Number of nodes )\n" << nodes << "\n";
318
319  app << "#N_LEAVES ( Number of leaves )\n" << Leaves() << "\n";
320
321  app << "#N_SPLITS ( Number of splits in axes x y z dx dy dz)\n";
322  for (int i=0; i<7; i++)
323    app << splits[i] <<" ";
324  app <<endl;
325
326  app << "#N_RAYREFS ( Number of rayRefs )\n" <<
327    rayRefs << "\n";
328
329  app << "#N_RAYRAYREFS  ( Number of rayRefs / ray )\n" <<
330    rayRefs/(double)rays << "\n";
331
332  app << "#N_LEAFRAYREFS  ( Number of rayRefs / leaf )\n" <<
333    rayRefs/(double)Leaves() << "\n";
334
335  app << "#N_MAXOBJECTREFS  ( Max number of rayRefs / leaf )\n" <<
336    maxObjectRefs << "\n";
337
338  app << "#N_NONEMPTYRAYREFS  ( Number of rayRefs in nonEmpty leaves / non empty leaf )\n" <<
339    rayRefsNonZeroQuery/(double)(Leaves() - zeroQueryNodes) << "\n";
340
341  app << "#N_LEAFDOMAINREFS  ( Number of query domain Refs / leaf )\n" <<
342    objectRefs/(double)Leaves() << "\n";
343
344  //  app << setprecision(4);
345
346  app << "#N_PEMPTYLEAVES  ( Percentage of leaves with zero query domains )\n"<<
347    zeroQueryNodes*100/(double)Leaves()<<endl;
348
349  app << "#N_PMAXDEPTHLEAVES ( Percentage of leaves at maxdepth )\n"<<
350    maxDepthNodes*100/(double)Leaves()<<endl;
351
352  app << "#N_PMINCOSTLEAVES  ( Percentage of leaves with minCost )\n"<<
353    minCostNodes*100/(double)Leaves()<<endl;
354
355  app << "#N_ADDED_RAYREFS  (Number of dynamically added ray references )\n"<<
356    addedRayRefs<<endl;
357
358  app << "#N_REMOVED_RAYREFS  (Number of dynamically removed ray references )\n"<<
359    removedRayRefs<<endl;
360
361  //  app << setprecision(4);
362
363  //  app << "#N_CTIME  ( Construction time [s] )\n"
364  //      << Time() << " \n";
365
366  app << "===== END OF KdTree statistics ==========\n";
367
368}
369
370
371
372void
373KdTree::EvaluateLeafStats(const TraversalData &data)
374{
375
376  // the node became a leaf -> evaluate stats for leafs
377  KdLeaf *leaf = (KdLeaf *)data.mNode;
378
379  if (data.mDepth > mTermMaxDepth)
380    mStat.maxDepthNodes++;
381 
382  if ( (int)(leaf->mObjects.size()) < mTermMinCost)
383    mStat.minCostNodes++;
384 
385 
386  if ( (int)(leaf->mObjects.size()) > mStat.maxObjectRefs)
387    mStat.maxObjectRefs = (int)leaf->mObjects.size();
388 
389}
390
391
392
393void
394KdTree::SortSplitCandidates(
395                            KdLeaf *node,
396                            const int axis
397                            )
398{
399  splitCandidates->clear();
400 
401  int requestedSize = 2*(int)node->mObjects.size();
402  // creates a sorted split candidates array
403  if (splitCandidates->capacity() > 500000 &&
404      requestedSize < (int)(splitCandidates->capacity()/10) ) {
405    delete splitCandidates;
406    splitCandidates = new vector<SortableEntry>;
407  }
408 
409  splitCandidates->reserve(requestedSize);
410 
411  // insert all queries
412  for(ObjectContainer::const_iterator mi = node->mObjects.begin();
413      mi != node->mObjects.end();
414      mi++) {
415    AxisAlignedBox3 box = (*mi)->GetBox();
416
417    splitCandidates->push_back(SortableEntry(SortableEntry::BOX_MIN,
418                                                                                         box.Min(axis),
419                                                                                         *mi)
420                                                           );
421   
422   
423    splitCandidates->push_back(SortableEntry(SortableEntry::BOX_MAX,
424                                                                                         box.Max(axis),
425                                                                                         *mi)
426                                                           );
427  }
428 
429  stable_sort(splitCandidates->begin(), splitCandidates->end());
430}
431
432
433float
434KdTree::BestCostRatio(
435                                          KdLeaf *node,
436                                          const AxisAlignedBox3 &box,
437                                          const int axis,
438                                          float &position,
439                                          int &objectsBack,
440                                          int &objectsFront
441                                          )
442{
443
444  SortSplitCandidates(node, axis);
445 
446  // go through the lists, count the number of objects left and right
447  // and evaluate the following cost funcion:
448  // C = ct_div_ci  + (ol + or)/queries
449
450  float totalIntersections = 0.0f;
451  vector<SortableEntry>::const_iterator ci;
452
453  for(ci = splitCandidates->begin();
454      ci < splitCandidates->end();
455      ci++)
456    if ((*ci).type == SortableEntry::BOX_MIN) {
457      totalIntersections += (*ci).intersectable->IntersectionComplexity();
458    }
459       
460  float intersectionsLeft = 0;
461  float intersectionsRight = totalIntersections;
462       
463  int objectsLeft = 0, objectsRight = (int)node->mObjects.size();
464 
465  float minBox = box.Min(axis);
466  float maxBox = box.Max(axis);
467  float boxArea = box.SurfaceArea();
468 
469  float minBand = minBox + mSplitBorder*(maxBox - minBox);
470  float maxBand = minBox + (1.0f - mSplitBorder)*(maxBox - minBox);
471 
472  float minSum = 1e20f;
473 
474  for(ci = splitCandidates->begin();
475      ci < splitCandidates->end();
476      ci++) {
477    switch ((*ci).type) {
478    case SortableEntry::BOX_MIN:
479      objectsLeft++;
480      intersectionsLeft += (*ci).intersectable->IntersectionComplexity();
481      break;
482    case SortableEntry::BOX_MAX:
483      objectsRight--;
484      intersectionsRight -= (*ci).intersectable->IntersectionComplexity();
485      break;
486    }
487
488    if ((*ci).value > minBand && (*ci).value < maxBand) {
489      AxisAlignedBox3 lbox = box;
490      AxisAlignedBox3 rbox = box;
491      lbox.SetMax(axis, (*ci).value);
492      rbox.SetMin(axis, (*ci).value);
493
494      float sum;
495      if (mSahUseFaces)
496                sum = intersectionsLeft*lbox.SurfaceArea() + intersectionsRight*rbox.SurfaceArea();
497      else
498                sum = objectsLeft*lbox.SurfaceArea() + objectsRight*rbox.SurfaceArea();
499     
500      //      cout<<"pos="<<(*ci).value<<"\t q=("<<ql<<","<<qr<<")\t r=("<<rl<<","<<rr<<")"<<endl;
501      //      cout<<"cost= "<<sum<<endl;
502     
503      if (sum < minSum) {
504                minSum = sum;
505                position = (*ci).value;
506               
507                objectsBack = objectsLeft;
508                objectsFront = objectsRight;
509      }
510    }
511  }
512 
513  float oldCost = mSahUseFaces ? totalIntersections : node->mObjects.size();
514  float newCost = mCt_div_ci + minSum/boxArea;
515  float ratio = newCost/oldCost;
516 
517#if 0
518  cout<<"===================="<<endl;
519  cout<<"costRatio="<<ratio<<" pos="<<position<<" t="<<(position - minBox)/(maxBox - minBox)
520      <<"\t o=("<<objectsBack<<","<<objectsFront<<")"<<endl;
521#endif
522  return ratio;
523}
524
525int
526KdTree::CastRay(
527                                Ray &ray
528                                )
529{
530  int hits = 0;
531 
532  stack<RayTraversalData> tStack;
533 
534  float maxt = 1e6;
535  float mint = 0;
536 
537  Intersectable::NewMail();
538
539  if (!mBox.GetMinMaxT(ray, &mint, &maxt))
540    return 0;
541
542  if (mint < 0)
543    mint = 0;
544 
545  maxt += Limits::Threshold;
546 
547  Vector3 entp = ray.Extrap(mint);
548  Vector3 extp = ray.Extrap(maxt);
549 
550  KdNode *node = mRoot;
551  KdNode *farChild;
552  float position;
553  int axis;
554
555 
556  while (1) {
557    if (!node->IsLeaf()) {
558      KdInterior *in = (KdInterior *) node;
559      position = in->mPosition;
560      axis = in->mAxis;
561
562      if (entp[axis] <= position) {
563                if (extp[axis] <= position) {
564                  node = in->mBack;
565                  // cases N1,N2,N3,P5,Z2,Z3
566                  continue;
567                } else {
568                  // case N4
569                  node = in->mBack;
570                  farChild = in->mFront;
571                }
572      }
573      else {
574                if (position <= extp[axis]) {
575                  node = in->mFront;
576                  // cases P1,P2,P3,N5,Z1
577                  continue;
578                } else {
579                  node = in->mFront;
580                  farChild = in->mBack;
581                  // case P4
582                }
583          }
584      // $$ modification 3.5.2004 - hints from Kamil Ghais
585      // case N4 or P4
586      float tdist = (position - ray.GetLoc(axis)) / ray.GetDir(axis);
587      tStack.push(RayTraversalData(farChild, extp, maxt));
588      extp = ray.GetLoc() + ray.GetDir()*tdist;
589      maxt = tdist;
590        } else {
591          // compute intersection with all objects in this leaf
592          KdLeaf *leaf = (KdLeaf *) node;
593          if (ray.mFlags & Ray::STORE_KDLEAVES)
594                ray.kdLeaves.push_back(leaf);
595         
596          ObjectContainer::const_iterator mi;
597          for ( mi = leaf->mObjects.begin();
598                        mi != leaf->mObjects.end();
599                        mi++) {
600                Intersectable *object = *mi;
601                if (!object->Mailed() ) {
602                  object->Mail();
603                  if (ray.mFlags & Ray::STORE_TESTED_OBJECTS)
604                        ray.testedObjects.push_back(object);
605
606                  static int oi=1;
607                  if (MeshDebug)
608                        cout<<"Object "<<oi++;
609                 
610                  hits += object->CastRay(ray);
611
612                  if (MeshDebug) {
613                        if (ray.intersections.size())
614                        cout<<"nearest t="<<ray.intersections[0].mT<<endl;
615                  else
616                        cout<<"nearest t=-INF"<<endl;
617                  }
618                 
619                }
620          }
621         
622          if (hits && ray.GetType() == Ray::LOCAL_RAY)
623                if (ray.intersections[0].mT <= maxt)
624                  break;
625         
626          // get the next node from the stack
627          if (tStack.empty())
628                break;
629         
630          entp = extp;
631          mint = maxt;
632          if (ray.GetType() == Ray::LINE_SEGMENT && mint > 1.0f)
633                break;
634         
635          RayTraversalData &s  = tStack.top();
636          node = s.mNode;
637          extp = s.mExitPoint;
638          maxt = s.mMaxT;
639          tStack.pop();
640        }
641  }
642  return hits;
643}
644
645int KdTree::CastLineSegment(const Vector3 &origin,
646                                                        const Vector3 &termination,
647                                                        ViewCellContainer &viewcells)
648{
649        int hits = 0;
650
651        float mint = 0.0f, maxt = 1.0f;
652        const Vector3 dir = termination - origin;
653
654        stack<RayTraversalData> tStack;
655
656        Intersectable::NewMail();
657
658        //maxt += Limits::Threshold;
659
660        Vector3 entp = origin;
661        Vector3 extp = termination;
662
663        KdNode *node = mRoot;
664        KdNode *farChild;
665
666        float position;
667        int axis;
668
669        while (1)
670        {
671                if (!node->IsLeaf())
672                {
673                        KdInterior *in = dynamic_cast<KdInterior *>(node);
674                        position = in->mPosition;
675                        axis = in->mAxis;
676
677                        if (entp[axis] <= position)
678                        {
679                                if (extp[axis] <= position)
680                                {
681                                        node = in->mBack;
682                                        // cases N1,N2,N3,P5,Z2,Z3
683                                        continue;
684                                }
685                                else
686                                {
687                                        // case N4
688                                        node = in->mBack;
689                                        farChild = in->mFront;
690                                }
691                        }
692                        else
693                        {
694                                if (position <= extp[axis])
695                                {
696                                        node = in->mFront;
697                                        // cases P1,P2,P3,N5,Z1
698                                        continue;
699                                }
700                                else
701                                {
702                                        node = in->mFront;
703                                        farChild = in->mBack;
704                                        // case P4
705                                }
706                        }
707
708                        // $$ modification 3.5.2004 - hints from Kamil Ghais
709                        // case N4 or P4
710                        float tdist = (position - origin[axis]) / dir[axis];
711                        //tStack.push(RayTraversalData(farChild, extp, maxt)); //TODO
712                        extp = origin + dir * tdist;
713                        maxt = tdist;
714                }
715                else
716                {
717                        // compute intersection with all objects in this leaf
718                        KdLeaf *leaf = dynamic_cast<KdLeaf *>(node);
719
720                        // add view cell to intersections
721                        ViewCell *vc = leaf->mViewCell;
722
723                        if (!vc->Mailed())
724                        {
725                                vc->Mail();
726                                viewcells.push_back(vc);
727                                ++ hits;
728                        }
729
730                        // get the next node from the stack
731                        if (tStack.empty())
732                                break;
733
734                        entp = extp;
735                        mint = maxt;
736                       
737                        RayTraversalData &s  = tStack.top();
738                        node = s.mNode;
739                        extp = s.mExitPoint;
740                        maxt = s.mMaxT;
741                        tStack.pop();
742                }
743        }
744
745        return hits;
746}
747
748
749void
750KdTree::CollectObjects(const AxisAlignedBox3 &box,
751                                           ObjectContainer &objects)
752{
753  stack<KdNode *> nodeStack;
754
755  nodeStack.push(mRoot);
756
757  while (!nodeStack.empty()) {
758    KdNode *node = nodeStack.top();
759    nodeStack.pop();
760    if (node->IsLeaf()) {
761      KdLeaf *leaf = (KdLeaf *)node;
762      for (int j=0; j < leaf->mObjects.size(); j++) {
763                Intersectable *object = leaf->mObjects[j];
764                if (!object->Mailed() && Overlap(box, object->GetBox())) {
765                  object->Mail();
766                  objects.push_back(object);
767                }
768      }
769    } else {
770      KdInterior *interior = (KdInterior *)node;
771
772          if ( box.Max()[interior->mAxis] > interior->mPosition )
773                nodeStack.push(interior->mFront);
774 
775          if (box.Min()[interior->mAxis] < interior->mPosition)
776                nodeStack.push(interior->mBack);
777    }
778  }
779}
780
781void
782KdTree::CollectObjects(KdNode *n, ObjectContainer &objects)
783{
784  stack<KdNode *> nodeStack;
785
786  nodeStack.push(n);
787
788  while (!nodeStack.empty()) {
789    KdNode *node = nodeStack.top();
790    nodeStack.pop();
791    if (node->IsLeaf()) {
792      KdLeaf *leaf = (KdLeaf *)node;
793      for (int j=0; j < leaf->mObjects.size(); j++) {
794                                Intersectable *object = leaf->mObjects[j];
795                                if (!object->Mailed()) {
796                                        object->Mail();
797                                        objects.push_back(object);
798                                }
799      }
800    } else {
801      KdInterior *interior = (KdInterior *)node;
802      nodeStack.push(interior->mFront);
803      nodeStack.push(interior->mBack);
804    }
805  }
806}
807
808// Find random neighbor which was not mailed
809KdNode *
810KdTree::FindRandomNeighbor(KdNode *n,
811                                                                                                         bool onlyUnmailed
812                                                                                                         )
813{
814  stack<KdNode *> nodeStack;
815 
816  nodeStack.push(mRoot);
817
818  AxisAlignedBox3 box = GetBox(n);
819  int mask = rand();
820
821  while (!nodeStack.empty()) {
822    KdNode *node = nodeStack.top();
823    nodeStack.pop();
824    if (node->IsLeaf()) {
825      if ( node != n && (!onlyUnmailed || !node->Mailed()) )
826        return node;
827    } else {
828      KdInterior *interior = (KdInterior *)node;
829      if (interior->mPosition > box.Max(interior->mAxis))
830        nodeStack.push(interior->mBack);
831      else
832        if (interior->mPosition < box.Min(interior->mAxis))
833          nodeStack.push(interior->mFront);
834        else {
835          // random decision
836          if (mask&1)
837            nodeStack.push(interior->mBack);
838          else
839            nodeStack.push(interior->mFront);
840          mask = mask>>1;
841        }
842    }
843  }
844 
845  return NULL;
846}
847
848int
849KdTree::FindNeighbors(KdNode *n,
850                      vector<KdNode *> &neighbors,
851                      bool onlyUnmailed
852                      )
853{
854  stack<KdNode *> nodeStack;
855 
856  nodeStack.push(mRoot);
857
858  AxisAlignedBox3 box = GetBox(n);
859
860  while (!nodeStack.empty()) {
861    KdNode *node = nodeStack.top();
862    nodeStack.pop();
863    if (node->IsLeaf()) {
864      if ( node != n && (!onlyUnmailed || !node->Mailed()) )
865        neighbors.push_back(node);
866    } else {
867      KdInterior *interior = (KdInterior *)node;
868      if (interior->mPosition > box.Max(interior->mAxis))
869                                nodeStack.push(interior->mBack);
870      else
871                                if (interior->mPosition < box.Min(interior->mAxis))
872                                        nodeStack.push(interior->mFront);
873                                else {
874                                        // random decision
875                                        nodeStack.push(interior->mBack);
876                                        nodeStack.push(interior->mFront);
877                                }
878    }
879  }
880 
881  return (int)neighbors.size();
882}
883
884// Find random neighbor which was not mailed
885KdNode *
886KdTree::GetRandomLeaf(const Plane3 &plane)
887{
888  stack<KdNode *> nodeStack;
889 
890  nodeStack.push(mRoot);
891 
892  int mask = rand();
893 
894  while (!nodeStack.empty()) {
895    KdNode *node = nodeStack.top();
896    nodeStack.pop();
897    if (node->IsLeaf()) {
898      return node;
899    } else {
900      KdInterior *interior = (KdInterior *)node;
901      KdNode *next;
902        if (GetBox(interior->mBack).Side(plane) < 0)
903          next = interior->mFront;
904        else
905          if (GetBox(interior->mFront).Side(plane) < 0)
906            next = interior->mBack;
907          else {
908            // random decision
909            if (mask&1)
910              next = interior->mBack;
911            else
912              next = interior->mFront;
913            mask = mask>>1;
914          }
915        nodeStack.push(next);
916    }
917  }
918 
919 
920  return NULL;
921}
922
923void
924KdTree::CollectLeaves(vector<KdLeaf *> &leaves)
925{
926  stack<KdNode *> nodeStack;
927  nodeStack.push(mRoot);
928
929  while (!nodeStack.empty()) {
930    KdNode *node = nodeStack.top();
931    nodeStack.pop();
932    if (node->IsLeaf()) {
933      KdLeaf *leaf = (KdLeaf *)node;
934      leaves.push_back(leaf);
935    } else {
936      KdInterior *interior = (KdInterior *)node;
937      nodeStack.push(interior->mBack);
938      nodeStack.push(interior->mFront);
939    }
940  }
941}
942
943void
944KdTree::CreateAndCollectViewCells(ViewCellContainer &vc) const
945{
946  stack<KdNode *> nodeStack;
947  nodeStack.push(mRoot);
948
949  while (!nodeStack.empty()) {
950    KdNode *node = nodeStack.top();
951    nodeStack.pop();
952    if (node->IsLeaf()) {
953      KdLeaf *leaf = (KdLeaf *)node;
954          // kdtree used as view cell container => create view cell
955          KdViewCell *viewCell = new KdViewCell();
956          leaf->mViewCell = viewCell;
957          // push back pointer to this leaf
958          viewCell->mLeaf = leaf;
959      vc.push_back(viewCell);
960    } else {
961      KdInterior *interior = (KdInterior *)node;
962      nodeStack.push(interior->mBack);
963      nodeStack.push(interior->mFront);
964    }
965  }
966}
967
968int
969KdTree::CollectLeafPvs()
970{
971  int totalPvsSize = 0;
972  stack<KdNode *> nodeStack;
973 
974  nodeStack.push(mRoot);
975 
976  while (!nodeStack.empty()) {
977    KdNode *node = nodeStack.top();
978    nodeStack.pop();
979    if (node->IsLeaf()) {
980      KdLeaf *leaf = (KdLeaf *)node;
981      for (int j=0; j < leaf->mObjects.size(); j++) {
982        Intersectable *object = leaf->mObjects[j];
983        if (!object->Mailed()) {
984          object->Mail();
985          // add this node to pvs of all nodes it can see
986          KdPvsMap::iterator ni = object->mKdPvs.mEntries.begin();
987          for (; ni != object->mKdPvs.mEntries.end(); ni++) {
988            KdNode *node = (*ni).first;
989            // $$ JB TEMPORARY solution -> should add object PVS or explictly computed
990            // kd tree PVS
991                float contribution;
992                if (leaf->mKdPvs.AddSample(node, 1.0f, contribution))
993              totalPvsSize++;
994          }
995        }
996      }
997    } else {
998      KdInterior *interior = (KdInterior *)node;
999      nodeStack.push(interior->mFront);
1000      nodeStack.push(interior->mBack);
1001    }
1002  }
1003
1004  return totalPvsSize;
1005}
1006
1007
1008KdNode *
1009KdTree::GetRandomLeaf(const bool onlyUnmailed)
1010{
1011  stack<KdNode *> nodeStack;
1012  nodeStack.push(mRoot);
1013 
1014  int mask = rand();
1015       
1016  while (!nodeStack.empty()) {
1017    KdNode *node = nodeStack.top();
1018    nodeStack.pop();
1019    if (node->IsLeaf()) {
1020      if ( (!onlyUnmailed || !node->Mailed()) )
1021                                return node;
1022    } else {
1023      KdInterior *interior = (KdInterior *)node;
1024                        // random decision
1025                        if (mask&1)
1026                                nodeStack.push(interior->mBack);
1027                        else
1028                                nodeStack.push(interior->mFront);
1029                        mask = mask>>1;
1030                }
1031        }
1032  return NULL;
1033}
1034
1035
1036int
1037KdTree::CastBeam(
1038                                 Beam &beam
1039                                 )
1040{
1041  stack<KdNode *> nodeStack;
1042  nodeStack.push(mRoot);
1043 
1044  while (!nodeStack.empty()) {
1045    KdNode *node = nodeStack.top();
1046    nodeStack.pop();
1047       
1048        int side = beam.ComputeIntersection(GetBox(node));
1049        switch (side) {
1050        case -1:
1051          beam.mKdNodes.push_back(node);
1052          break;
1053        case 0:
1054          if (node->IsLeaf())
1055                beam.mKdNodes.push_back(node);
1056          else {
1057                KdInterior *interior = (KdInterior *)node;
1058                KdNode *first = interior->mBack;
1059                KdNode *second = interior->mFront;
1060               
1061                if (interior->mAxis < 3) {
1062                  // spatial split -> decide on the order of the nodes
1063                  if (beam.mPlanes[0].mNormal[interior->mAxis] > 0)
1064                        swap(first, second);
1065                }
1066
1067                nodeStack.push(first);
1068                nodeStack.push(second);
1069          }
1070          break;
1071          // default: cull
1072        }
1073  }
1074
1075  if (beam.mFlags & Beam::STORE_OBJECTS)
1076  {
1077          vector<KdNode *>::const_iterator it, it_end = beam.mKdNodes.end();
1078       
1079          Intersectable::NewMail();
1080          for (it = beam.mKdNodes.begin(); it != it_end; ++ it)
1081          {
1082                  CollectObjects(*it, beam.mObjects);
1083          }
1084  }
1085
1086  return (int)beam.mKdNodes.size();
1087}
1088
1089}
Note: See TracBrowser for help on using the repository browser.