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

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