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

Revision 1414, 31.5 KB checked in by mattausch, 18 years ago (diff)

fixed kd tree loading / exporting

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