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

Revision 904, 24.5 KB checked in by bittner, 18 years ago (diff)

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