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

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