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

Revision 1594, 32.3 KB checked in by bittner, 18 years ago (diff)

support for kd tree based pvs

RevLine 
[372]1#include <stack>
2#include <algorithm>
3#include <queue>
4#include "Environment.h"
5#include "Mesh.h"
6#include "KdTree.h"
[469]7#include "ViewCell.h"
[505]8#include "Beam.h"
[372]9
10
11
[863]12namespace GtpVisibilityPreprocessor {
[860]13
[1176]14int KdNode::sMailId = 1;
15int KdNode::sReservedMailboxes = 1;
[863]16
[1197]17inline static bool ilt(Intersectable *obj1, Intersectable *obj2)
18{
19        return obj1->mId < obj2->mId;
20}
21
22
[1176]23KdNode::KdNode(KdInterior *parent):mParent(parent), mMailbox(0)
[372]24{
25  if (parent)
26    mDepth = parent->mDepth+1;
27  else
28    mDepth = 0;
29}
30
[1002]31
32KdInterior::~KdInterior()
33{
34        // recursivly destroy children
35        DEL_PTR(mFront);
36        DEL_PTR(mBack);
37}
38
39
[372]40KdTree::KdTree()
41{
[878]42
43 
[372]44  mRoot = new KdLeaf(NULL, 0);
[1004]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);
[372]51
[1004]52  Environment::GetSingleton()->GetBoolValue("KdTree.sahUseFaces", mSahUseFaces);
[372]53
54  char splitType[64];
[1004]55  Environment::GetSingleton()->GetStringValue("KdTree.splitMethod", splitType);
[372]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
[1002]73
74KdTree::~KdTree()
75{
76        DEL_PTR(mRoot);
77}
78
79
[372]80bool
81KdTree::Construct()
82{
[1076]83
[372]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++) {
[752]97        //      cout<<(*mi)->GetBox()<<endl;
98        mBox.Include((*mi)->GetBox());
[372]99  }
100
[752]101  cout <<"KdTree Root Box:"<<mBox<<endl;
[372]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()) {
[752]123        //      cout<<mStat.Nodes()<<" "<<mTermMaxNodes<<endl;
124        if (mStat.Nodes() > mTermMaxNodes) {
125          //    if ( GetMemUsage() > maxMemory ) {
[372]126      // count statistics on unprocessed leafs
127      while (!tStack.empty()) {
[752]128                EvaluateLeafStats(tStack.top());
129                tStack.pop();
[372]130      }
131      break;
132    }
[752]133         
134       
[372]135    TraversalData data = tStack.top();
136    tStack.pop();
137   
138    KdNode *node = SubdivideNode((KdLeaf *) data.mNode,
[752]139                                                                 data.mBox,
140                                                                 backBox,
141                                                                 frontBox
142                                                                 );
143
144        if (result == NULL)
[372]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
[1108]170    ((int)leaf->mObjects.size() <= mTermMinCost) ||
[372]171    (leaf->mDepth >= mTermMaxDepth);
172 
173}
174
175
176int
177KdTree::SelectPlane(KdLeaf *leaf,
[752]178                                        const AxisAlignedBox3 &box,
179                                        float &position
180                                        )
[372]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;
[1584]194      bool mOnlyDrivingAxis = true;
[1387]195
196          if (mOnlyDrivingAxis) {
[752]197                axis = box.Size().DrivingAxis();
198                costRatio = BestCostRatio(leaf,
199                                                                  box,
200                                                                  axis,
201                                                                  position,
202                                                                  objectsBack,
203                                                                  objectsFront);
[372]204      } else {
[752]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                }
[372]220      }
221     
222      if (costRatio > mMaxCostRatio) {
[752]223                //cout<<"Too big cost ratio "<<costRatio<<endl;
224                axis = -1;
[372]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))
[469]243          return leaf;
244
[372]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
[1074]292 
[372]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();
[1143]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
[372]311
[1143]312
313        if (box.Max(axis) >= position )
314        {
[372]315      front->mObjects.push_back(*mi);
[1143]316          ++ (*mi)->mReferences;
317        }
318       
[372]319    if (box.Min(axis) < position )
[1143]320        {
[372]321      back->mObjects.push_back(*mi);
[1143]322          ++ (*mi)->mReferences;
323        }
324
[1144]325    mStat.objectRefs -= (int)leaf->mObjects.size();
326    mStat.objectRefs += objectsBack + objectsFront;
327  }
328
[1143]329        // store objects referenced in more than one leaf
330        // for easy access
331        ProcessMultipleRefs(back);
332        ProcessMultipleRefs(front);
333
[1286]334        delete leaf;
335        return node;
[372]336}
337
338
[1143]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
[372]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
[667]365  app << "#N_SPLITS ( Number of splits in axes x y z dx dy dz)\n";
[372]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
[1184]379  app << "#N_MAXOBJECTREFS  ( Max number of object refs / leaf )\n" <<
[372]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)
[469]431    mStat.maxObjectRefs = (int)leaf->mObjects.size();
[372]432 
433}
434
435
436
437void
[1233]438KdTree::SortSubdivisionCandidates(
[372]439                            KdLeaf *node,
440                            const int axis
441                            )
442{
443  splitCandidates->clear();
444 
[469]445  int requestedSize = 2*(int)node->mObjects.size();
[372]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,
[752]462                                                                                         box.Min(axis),
463                                                                                         *mi)
464                                                           );
[372]465   
466   
467    splitCandidates->push_back(SortableEntry(SortableEntry::BOX_MAX,
[752]468                                                                                         box.Max(axis),
469                                                                                         *mi)
470                                                           );
[372]471  }
472 
473  stable_sort(splitCandidates->begin(), splitCandidates->end());
474}
475
476
477float
478KdTree::BestCostRatio(
[752]479                                          KdLeaf *node,
480                                          const AxisAlignedBox3 &box,
481                                          const int axis,
482                                          float &position,
483                                          int &objectsBack,
484                                          int &objectsFront
485                                          )
[372]486{
487
[1387]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 
[1233]508  SortSubdivisionCandidates(node, axis);
[372]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       
[469]527  int objectsLeft = 0, objectsRight = (int)node->mObjects.size();
[372]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 
[469]536  float minSum = 1e20f;
[372]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)
[752]560                sum = intersectionsLeft*lbox.SurfaceArea() + intersectionsRight*rbox.SurfaceArea();
[372]561      else
[752]562                sum = objectsLeft*lbox.SurfaceArea() + objectsRight*rbox.SurfaceArea();
[372]563     
564      //      cout<<"pos="<<(*ci).value<<"\t q=("<<ql<<","<<qr<<")\t r=("<<rl<<","<<rr<<")"<<endl;
565      //      cout<<"cost= "<<sum<<endl;
[1387]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         
[372]576      if (sum < minSum) {
[752]577                minSum = sum;
578                position = (*ci).value;
579               
580                objectsBack = objectsLeft;
581                objectsFront = objectsRight;
[372]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(
[462]600                                Ray &ray
601                                )
[372]602{
603  int hits = 0;
604 
605  stack<RayTraversalData> tStack;
606 
607  float maxt = 1e6;
608  float mint = 0;
[1583]609
[1584]610  //  ray.ComputeInvertedDir();
[1583]611
[372]612  Intersectable::NewMail();
613
614  if (!mBox.GetMinMaxT(ray, &mint, &maxt))
615    return 0;
616
617  if (mint < 0)
618    mint = 0;
[1584]619
[372]620 
621  maxt += Limits::Threshold;
622 
623  Vector3 entp = ray.Extrap(mint);
624  Vector3 extp = ray.Extrap(maxt);
625 
626  KdNode *node = mRoot;
627  KdNode *farChild;
628  float position;
629  int axis;
[752]630
[372]631 
632  while (1) {
633    if (!node->IsLeaf()) {
634      KdInterior *in = (KdInterior *) node;
635      position = in->mPosition;
636      axis = in->mAxis;
637
638      if (entp[axis] <= position) {
[752]639                if (extp[axis] <= position) {
640                  node = in->mBack;
641                  // cases N1,N2,N3,P5,Z2,Z3
642                  continue;
643                } else {
644                  // case N4
645                  node = in->mBack;
646                  farChild = in->mFront;
647                }
[372]648      }
649      else {
[752]650                if (position <= extp[axis]) {
651                  node = in->mFront;
652                  // cases P1,P2,P3,N5,Z1
653                  continue;
654                } else {
655                  node = in->mFront;
656                  farChild = in->mBack;
657                  // case P4
658                }
659          }
[372]660      // $$ modification 3.5.2004 - hints from Kamil Ghais
661      // case N4 or P4
662      float tdist = (position - ray.GetLoc(axis)) / ray.GetDir(axis);
663      tStack.push(RayTraversalData(farChild, extp, maxt));
664      extp = ray.GetLoc() + ray.GetDir()*tdist;
665      maxt = tdist;
[752]666        } else {
667          // compute intersection with all objects in this leaf
668          KdLeaf *leaf = (KdLeaf *) node;
669          if (ray.mFlags & Ray::STORE_KDLEAVES)
670                ray.kdLeaves.push_back(leaf);
671         
672          ObjectContainer::const_iterator mi;
673          for ( mi = leaf->mObjects.begin();
674                        mi != leaf->mObjects.end();
675                        mi++) {
676                Intersectable *object = *mi;
677                if (!object->Mailed() ) {
678                  object->Mail();
679                  if (ray.mFlags & Ray::STORE_TESTED_OBJECTS)
680                        ray.testedObjects.push_back(object);
[878]681
682                  static int oi=1;
683                  if (MeshDebug)
684                        cout<<"Object "<<oi++;
685                 
[752]686                  hits += object->CastRay(ray);
[878]687
688                  if (MeshDebug) {
[1344]689                        if (!ray.intersections.empty())
[878]690                        cout<<"nearest t="<<ray.intersections[0].mT<<endl;
691                  else
692                        cout<<"nearest t=-INF"<<endl;
[1344]693                  }       
[376]694                }
[752]695          }
696         
697          if (hits && ray.GetType() == Ray::LOCAL_RAY)
698                if (ray.intersections[0].mT <= maxt)
699                  break;
700         
701          // get the next node from the stack
702          if (tStack.empty())
703                break;
704         
705          entp = extp;
706          mint = maxt;
707          if (ray.GetType() == Ray::LINE_SEGMENT && mint > 1.0f)
708                break;
709         
710          RayTraversalData &s  = tStack.top();
711          node = s.mNode;
712          extp = s.mExitPoint;
713          maxt = s.mMaxT;
714          tStack.pop();
715        }
[372]716  }
717  return hits;
718}
719
[469]720int KdTree::CastLineSegment(const Vector3 &origin,
721                                                        const Vector3 &termination,
[882]722                                                        ViewCellContainer &viewcells)
[469]723{
724        int hits = 0;
725
726        float mint = 0.0f, maxt = 1.0f;
727        const Vector3 dir = termination - origin;
728
[475]729        stack<RayTraversalData> tStack;
[469]730
731        Intersectable::NewMail();
732
733        //maxt += Limits::Threshold;
734
735        Vector3 entp = origin;
736        Vector3 extp = termination;
737
738        KdNode *node = mRoot;
739        KdNode *farChild;
740
741        float position;
742        int axis;
743
744        while (1)
745        {
746                if (!node->IsLeaf())
747                {
748                        KdInterior *in = dynamic_cast<KdInterior *>(node);
749                        position = in->mPosition;
750                        axis = in->mAxis;
751
752                        if (entp[axis] <= position)
753                        {
754                                if (extp[axis] <= position)
755                                {
756                                        node = in->mBack;
757                                        // cases N1,N2,N3,P5,Z2,Z3
758                                        continue;
759                                }
760                                else
761                                {
762                                        // case N4
763                                        node = in->mBack;
764                                        farChild = in->mFront;
765                                }
766                        }
767                        else
768                        {
769                                if (position <= extp[axis])
770                                {
771                                        node = in->mFront;
772                                        // cases P1,P2,P3,N5,Z1
773                                        continue;
774                                }
775                                else
776                                {
777                                        node = in->mFront;
778                                        farChild = in->mBack;
779                                        // case P4
780                                }
781                        }
782
783                        // $$ modification 3.5.2004 - hints from Kamil Ghais
784                        // case N4 or P4
785                        float tdist = (position - origin[axis]) / dir[axis];
786                        //tStack.push(RayTraversalData(farChild, extp, maxt)); //TODO
787                        extp = origin + dir * tdist;
788                        maxt = tdist;
789                }
790                else
791                {
792                        // compute intersection with all objects in this leaf
793                        KdLeaf *leaf = dynamic_cast<KdLeaf *>(node);
794
795                        // add view cell to intersections
796                        ViewCell *vc = leaf->mViewCell;
797
798                        if (!vc->Mailed())
799                        {
800                                vc->Mail();
801                                viewcells.push_back(vc);
802                                ++ hits;
803                        }
804
805                        // get the next node from the stack
806                        if (tStack.empty())
807                                break;
808
809                        entp = extp;
810                        mint = maxt;
811                       
812                        RayTraversalData &s  = tStack.top();
813                        node = s.mNode;
814                        extp = s.mExitPoint;
815                        maxt = s.mMaxT;
816                        tStack.pop();
817                }
818        }
819
820        return hits;
821}
822
823
[372]824void
[859]825KdTree::CollectObjects(const AxisAlignedBox3 &box,
826                                           ObjectContainer &objects)
827{
828  stack<KdNode *> nodeStack;
829
830  nodeStack.push(mRoot);
831
832  while (!nodeStack.empty()) {
833    KdNode *node = nodeStack.top();
834    nodeStack.pop();
835    if (node->IsLeaf()) {
836      KdLeaf *leaf = (KdLeaf *)node;
837      for (int j=0; j < leaf->mObjects.size(); j++) {
838                Intersectable *object = leaf->mObjects[j];
[904]839                if (!object->Mailed() && Overlap(box, object->GetBox())) {
[859]840                  object->Mail();
841                  objects.push_back(object);
842                }
843      }
844    } else {
845      KdInterior *interior = (KdInterior *)node;
846
847          if ( box.Max()[interior->mAxis] > interior->mPosition )
848                nodeStack.push(interior->mFront);
849 
850          if (box.Min()[interior->mAxis] < interior->mPosition)
851                nodeStack.push(interior->mBack);
852    }
853  }
854}
855
856void
[372]857KdTree::CollectObjects(KdNode *n, ObjectContainer &objects)
858{
859  stack<KdNode *> nodeStack;
860
861  nodeStack.push(n);
862
863  while (!nodeStack.empty()) {
864    KdNode *node = nodeStack.top();
865    nodeStack.pop();
866    if (node->IsLeaf()) {
867      KdLeaf *leaf = (KdLeaf *)node;
868      for (int j=0; j < leaf->mObjects.size(); j++) {
[374]869                                Intersectable *object = leaf->mObjects[j];
870                                if (!object->Mailed()) {
871                                        object->Mail();
872                                        objects.push_back(object);
873                                }
[372]874      }
875    } else {
876      KdInterior *interior = (KdInterior *)node;
877      nodeStack.push(interior->mFront);
878      nodeStack.push(interior->mBack);
879    }
880  }
881}
882
883// Find random neighbor which was not mailed
884KdNode *
885KdTree::FindRandomNeighbor(KdNode *n,
[1286]886                                                   bool onlyUnmailed
887                                                   )
[372]888{
889  stack<KdNode *> nodeStack;
890 
891  nodeStack.push(mRoot);
892
893  AxisAlignedBox3 box = GetBox(n);
894  int mask = rand();
895
896  while (!nodeStack.empty()) {
897    KdNode *node = nodeStack.top();
898    nodeStack.pop();
899    if (node->IsLeaf()) {
900      if ( node != n && (!onlyUnmailed || !node->Mailed()) )
901        return node;
902    } else {
903      KdInterior *interior = (KdInterior *)node;
904      if (interior->mPosition > box.Max(interior->mAxis))
905        nodeStack.push(interior->mBack);
906      else
907        if (interior->mPosition < box.Min(interior->mAxis))
908          nodeStack.push(interior->mFront);
909        else {
910          // random decision
911          if (mask&1)
912            nodeStack.push(interior->mBack);
913          else
914            nodeStack.push(interior->mFront);
915          mask = mask>>1;
916        }
917    }
918  }
919 
920  return NULL;
921}
922
923int
924KdTree::FindNeighbors(KdNode *n,
925                      vector<KdNode *> &neighbors,
926                      bool onlyUnmailed
927                      )
928{
929  stack<KdNode *> nodeStack;
930 
931  nodeStack.push(mRoot);
932
933  AxisAlignedBox3 box = GetBox(n);
934
935  while (!nodeStack.empty()) {
936    KdNode *node = nodeStack.top();
937    nodeStack.pop();
938    if (node->IsLeaf()) {
939      if ( node != n && (!onlyUnmailed || !node->Mailed()) )
940        neighbors.push_back(node);
941    } else {
942      KdInterior *interior = (KdInterior *)node;
943      if (interior->mPosition > box.Max(interior->mAxis))
[374]944                                nodeStack.push(interior->mBack);
[372]945      else
[374]946                                if (interior->mPosition < box.Min(interior->mAxis))
947                                        nodeStack.push(interior->mFront);
948                                else {
949                                        // random decision
950                                        nodeStack.push(interior->mBack);
951                                        nodeStack.push(interior->mFront);
952                                }
[372]953    }
954  }
955 
[469]956  return (int)neighbors.size();
[372]957}
958
959// Find random neighbor which was not mailed
960KdNode *
961KdTree::GetRandomLeaf(const Plane3 &plane)
962{
963  stack<KdNode *> nodeStack;
964 
965  nodeStack.push(mRoot);
966 
967  int mask = rand();
968 
969  while (!nodeStack.empty()) {
970    KdNode *node = nodeStack.top();
971    nodeStack.pop();
972    if (node->IsLeaf()) {
973      return node;
974    } else {
975      KdInterior *interior = (KdInterior *)node;
976      KdNode *next;
977        if (GetBox(interior->mBack).Side(plane) < 0)
978          next = interior->mFront;
979        else
980          if (GetBox(interior->mFront).Side(plane) < 0)
981            next = interior->mBack;
982          else {
983            // random decision
984            if (mask&1)
985              next = interior->mBack;
986            else
987              next = interior->mFront;
988            mask = mask>>1;
989          }
990        nodeStack.push(next);
991    }
992  }
993 
994 
995  return NULL;
996}
997
998void
999KdTree::CollectLeaves(vector<KdLeaf *> &leaves)
1000{
1001  stack<KdNode *> nodeStack;
1002  nodeStack.push(mRoot);
1003
1004  while (!nodeStack.empty()) {
1005    KdNode *node = nodeStack.top();
1006    nodeStack.pop();
1007    if (node->IsLeaf()) {
1008      KdLeaf *leaf = (KdLeaf *)node;
1009      leaves.push_back(leaf);
1010    } else {
1011      KdInterior *interior = (KdInterior *)node;
1012      nodeStack.push(interior->mBack);
1013      nodeStack.push(interior->mFront);
1014    }
1015  }
1016}
1017
[469]1018void
1019KdTree::CreateAndCollectViewCells(ViewCellContainer &vc) const
1020{
1021  stack<KdNode *> nodeStack;
1022  nodeStack.push(mRoot);
[372]1023
[469]1024  while (!nodeStack.empty()) {
1025    KdNode *node = nodeStack.top();
1026    nodeStack.pop();
1027    if (node->IsLeaf()) {
1028      KdLeaf *leaf = (KdLeaf *)node;
1029          // kdtree used as view cell container => create view cell
[471]1030          KdViewCell *viewCell = new KdViewCell();
1031          leaf->mViewCell = viewCell;
1032          // push back pointer to this leaf
[1551]1033          viewCell->mLeaves[0] = leaf;
[471]1034      vc.push_back(viewCell);
[469]1035    } else {
1036      KdInterior *interior = (KdInterior *)node;
1037      nodeStack.push(interior->mBack);
1038      nodeStack.push(interior->mFront);
1039    }
1040  }
1041}
1042
[372]1043int
1044KdTree::CollectLeafPvs()
1045{
1046  int totalPvsSize = 0;
1047  stack<KdNode *> nodeStack;
1048 
1049  nodeStack.push(mRoot);
1050 
1051  while (!nodeStack.empty()) {
1052    KdNode *node = nodeStack.top();
1053    nodeStack.pop();
1054    if (node->IsLeaf()) {
1055      KdLeaf *leaf = (KdLeaf *)node;
1056      for (int j=0; j < leaf->mObjects.size(); j++) {
1057        Intersectable *object = leaf->mObjects[j];
1058        if (!object->Mailed()) {
1059          object->Mail();
1060          // add this node to pvs of all nodes it can see
1061          KdPvsMap::iterator ni = object->mKdPvs.mEntries.begin();
1062          for (; ni != object->mKdPvs.mEntries.end(); ni++) {
1063            KdNode *node = (*ni).first;
1064            // $$ JB TEMPORARY solution -> should add object PVS or explictly computed
1065            // kd tree PVS
[466]1066                float contribution;
[556]1067                if (leaf->mKdPvs.AddSample(node, 1.0f, contribution))
[372]1068              totalPvsSize++;
1069          }
1070        }
1071      }
1072    } else {
1073      KdInterior *interior = (KdInterior *)node;
1074      nodeStack.push(interior->mFront);
1075      nodeStack.push(interior->mBack);
1076    }
1077  }
1078
1079  return totalPvsSize;
1080}
1081
1082
1083KdNode *
1084KdTree::GetRandomLeaf(const bool onlyUnmailed)
1085{
[507]1086  stack<KdNode *> nodeStack;
[372]1087  nodeStack.push(mRoot);
[507]1088 
1089  int mask = rand();
[372]1090       
1091  while (!nodeStack.empty()) {
1092    KdNode *node = nodeStack.top();
1093    nodeStack.pop();
1094    if (node->IsLeaf()) {
1095      if ( (!onlyUnmailed || !node->Mailed()) )
1096                                return node;
1097    } else {
1098      KdInterior *interior = (KdInterior *)node;
1099                        // random decision
1100                        if (mask&1)
1101                                nodeStack.push(interior->mBack);
1102                        else
1103                                nodeStack.push(interior->mFront);
1104                        mask = mask>>1;
1105                }
1106        }
1107  return NULL;
1108}
[504]1109
1110
1111int
[507]1112KdTree::CastBeam(
[512]1113                                 Beam &beam
[507]1114                                 )
[504]1115{
[507]1116  stack<KdNode *> nodeStack;
1117  nodeStack.push(mRoot);
[504]1118 
[507]1119  while (!nodeStack.empty()) {
1120    KdNode *node = nodeStack.top();
1121    nodeStack.pop();
1122       
1123        int side = beam.ComputeIntersection(GetBox(node));
1124        switch (side) {
1125        case -1:
1126          beam.mKdNodes.push_back(node);
1127          break;
1128        case 0:
1129          if (node->IsLeaf())
1130                beam.mKdNodes.push_back(node);
1131          else {
1132                KdInterior *interior = (KdInterior *)node;
1133                KdNode *first = interior->mBack;
1134                KdNode *second = interior->mFront;
1135               
1136                if (interior->mAxis < 3) {
1137                  // spatial split -> decide on the order of the nodes
1138                  if (beam.mPlanes[0].mNormal[interior->mAxis] > 0)
1139                        swap(first, second);
1140                }
[504]1141
[507]1142                nodeStack.push(first);
1143                nodeStack.push(second);
1144          }
1145          break;
1146          // default: cull
1147        }
1148  }
1149
[532]1150  if (beam.mFlags & Beam::STORE_OBJECTS)
1151  {
1152          vector<KdNode *>::const_iterator it, it_end = beam.mKdNodes.end();
[535]1153       
[532]1154          Intersectable::NewMail();
1155          for (it = beam.mKdNodes.begin(); it != it_end; ++ it)
1156          {
1157                  CollectObjects(*it, beam.mObjects);
1158          }
1159  }
1160
[837]1161  return (int)beam.mKdNodes.size();
[504]1162}
[860]1163
[1074]1164
[1196]1165#define TYPE_INTERIOR -2
1166#define TYPE_LEAF -3
[1194]1167
1168
[1201]1169void KdTree::ExportBinLeaf(OUT_STREAM &stream, KdLeaf *leaf)
[1194]1170{
1171        ObjectContainer::const_iterator it, it_end = leaf->mObjects.end();
1172       
[1196]1173        int type = TYPE_LEAF;
1174        int size = (int)leaf->mObjects.size();
1175
1176        stream.write(reinterpret_cast<char *>(&type), sizeof(int));
1177        stream.write(reinterpret_cast<char *>(&size), sizeof(int));
1178
[1194]1179        for (it = leaf->mObjects.begin(); it != it_end; ++ it)
1180        {       
1181                Intersectable *obj = *it;
1182                int id = obj->mId;             
[1196]1183       
[1194]1184                //stream.write(reinterpret_cast<char *>(&origin), sizeof(Vector3));
1185                stream.write(reinterpret_cast<char *>(&id), sizeof(int));
1186    }
[878]1187}
[1194]1188
1189
[1201]1190KdLeaf *KdTree::ImportBinLeaf(IN_STREAM &stream,
[1197]1191                                                          KdInterior *parent,
1192                                                          const ObjectContainer &objects)
[1194]1193{
[1196]1194        int leafId = TYPE_LEAF;
[1194]1195        int objId = leafId;
[1196]1196        int size;
1197
1198        stream.read(reinterpret_cast<char *>(&size), sizeof(int));
1199        KdLeaf *leaf = new KdLeaf(parent, size);
1200
[1197]1201        MeshInstance dummyInst(NULL);
1202        //Debug << "objects size: " << size << endl;
1203
[1196]1204        // read object ids
1205        for (int i = 0; i < size; ++ i)
[1194]1206        {       
1207                stream.read(reinterpret_cast<char *>(&objId), sizeof(int));
[1197]1208                dummyInst.SetId(objId);
[1196]1209
[1197]1210                ObjectContainer::const_iterator oit =
1211                        lower_bound(objects.begin(), objects.end(), (Intersectable *)&dummyInst, ilt);
1212                                                               
1213                if ((oit != objects.end()) && ((*oit)->GetId() == objId))
1214                {
1215                        leaf->mObjects.push_back(*oit);
1216                }
1217                else
1218                {
1219                        Debug << "error: object with id " << objId << " does not exist" << endl;
1220                }
1221        }
1222
[1196]1223        return leaf;
[1194]1224}
1225
1226
[1201]1227void KdTree::ExportBinInterior(OUT_STREAM &stream, KdInterior *interior)
[1194]1228{
[1196]1229        int interiorid = TYPE_INTERIOR;
[1194]1230        stream.write(reinterpret_cast<char *>(&interiorid), sizeof(int));
1231
1232        int axis = interior->mAxis;
[1196]1233        float pos = interior->mPosition;
[1194]1234
1235        stream.write(reinterpret_cast<char *>(&axis), sizeof(int));
[1196]1236        stream.write(reinterpret_cast<char *>(&pos), sizeof(float));
[1194]1237}
1238
1239
[1201]1240KdInterior *KdTree::ImportBinInterior(IN_STREAM  &stream, KdInterior *parent)
[1194]1241{
[1196]1242        KdInterior *interior = new KdInterior(parent);
[1194]1243
1244        int axis = interior->mAxis;
[1196]1245        float pos = interior->mPosition;
[1194]1246
1247        stream.read(reinterpret_cast<char *>(&axis), sizeof(int));
[1196]1248        stream.read(reinterpret_cast<char *>(&pos), sizeof(float));
1249
1250        interior->mAxis = axis;
1251        interior->mPosition = pos;
1252
1253        return interior;
[1194]1254}
1255
1256
1257bool KdTree::ExportBinTree(const string &filename)
1258{
[1201]1259        OUT_STREAM stream(filename.c_str(), OUT_BIN_MODE);
[1194]1260       
[1201]1261        //if (!stream.is_open()) return false;
[1194]1262
1263        // export binary version of mesh
[1197]1264        queue<KdNode *> tStack;
[1194]1265        tStack.push(mRoot);
1266
1267        while(!tStack.empty())
1268        {
[1197]1269                KdNode *node = tStack.front();
[1196]1270                tStack.pop();
1271               
1272                if (node->IsLeaf())
1273                {
[1197]1274                        //Debug << "l";
[1196]1275                        ExportBinLeaf(stream, dynamic_cast<KdLeaf *>(node));
1276                }
1277                else
1278                {
[1197]1279                        //Debug << "i";
[1196]1280                        KdInterior *interior = dynamic_cast<KdInterior *>(node);
[1194]1281
[1196]1282                        ExportBinInterior(stream, interior);
1283                       
1284                        tStack.push(interior->mFront);
1285                        tStack.push(interior->mBack);
[1194]1286                }
1287        }
1288
1289        return true;
1290}
1291
1292
[1201]1293KdNode *KdTree::LoadNextNode(IN_STREAM &stream,
[1197]1294                                                         KdInterior *parent,
1295                                                         const ObjectContainer &objects)
[1194]1296{
[1196]1297        int nodeType;
1298        stream.read(reinterpret_cast<char *>(&nodeType), sizeof(int));
1299
1300        if (nodeType == TYPE_LEAF)
1301        {
[1197]1302                return ImportBinLeaf(stream, dynamic_cast<KdInterior *>(parent), objects);
[1196]1303        }
1304
1305        if (nodeType == TYPE_INTERIOR)
1306        {
1307                return ImportBinInterior(stream, dynamic_cast<KdInterior *>(parent));
1308        }
1309
1310        Debug << "error! loading failed!" << endl;
[1194]1311       
[1196]1312        return NULL;
1313}
1314
1315
[1197]1316bool KdTree::LoadBinTree(const string &filename, ObjectContainer &objects)
[1196]1317{
1318        // export binary version of mesh
[1197]1319        queue<TraversalData> tStack;
[1201]1320        IN_STREAM stream(filename.c_str(), IN_BIN_MODE);
[1197]1321
[1286]1322        if (!stream.is_open()) return false;
[1194]1323
[1414]1324        // sort objects by their id
[1197]1325        std::stable_sort(objects.begin(), objects.end(), ilt);
[1194]1326
[1197]1327        mBox.Initialize();
[1196]1328        ObjectContainer::const_iterator oit, oit_end = objects.end();
[1194]1329
[1414]1330
1331        ///////////////////////////
1332        //-- compute bounding box of object space
[1196]1333    for (oit = objects.begin(); oit != oit_end; ++ oit)
1334        {
[1414]1335                const AxisAlignedBox3 box = (*oit)->GetBox();
1336                mBox.Include(box);
[1196]1337        }
1338
[1414]1339        // hack: we make a new root
1340        DEL_PTR(mRoot);
[1197]1341 
1342        KdNode *node = LoadNextNode(stream, NULL, objects);
[1196]1343        mRoot = node;
1344
1345        tStack.push(TraversalData(node, mBox, 0));
1346        mStat.Reset();
1347        mStat.nodes = 1;
1348
[1194]1349        while(!tStack.empty())
1350        {
[1197]1351                TraversalData tData = tStack.front();
[1196]1352                tStack.pop();
[1194]1353
[1196]1354                KdNode *node = tData.mNode;
1355
1356                if (!node->IsLeaf())
[1194]1357                {
[1196]1358                        mStat.nodes += 2;
1359
[1197]1360                        //Debug << "i" ;
[1194]1361                        KdInterior *interior = dynamic_cast<KdInterior *>(node);
[1196]1362                        interior->mBox = tData.mBox;
[1194]1363
[1197]1364            KdNode *front = LoadNextNode(stream, interior, objects);
1365                        KdNode *back = LoadNextNode(stream, interior, objects);
[1196]1366       
1367                        interior->SetupChildLinks(back, front);
1368
1369                        ++ mStat.splits[interior->mAxis];
1370
1371                        // compute new bounding box
1372                        AxisAlignedBox3 frontBox, backBox;
[1194]1373                       
[1196]1374                        tData.mBox.Split(interior->mAxis,
1375                                                         interior->mPosition,
1376                                                         frontBox,
1377                                                         backBox);
1378
1379                        tStack.push(TraversalData(front, frontBox, tData.mDepth + 1));                 
[1197]1380                        tStack.push(TraversalData(back, backBox, tData.mDepth + 1));
[1194]1381                }
[1196]1382                else
1383                {
1384                        EvaluateLeafStats(tData);
[1415]1385                        //cout << "l";
[1196]1386                }
[1194]1387        }
[1415]1388
[1196]1389        Debug << mStat << endl;
1390
[1194]1391        return true;
1392}
1393
[1594]1394KdIntersectable *
1395KdTree::GetOrCreateKdIntersectable(KdNode *node)
1396{
[1194]1397
[1594]1398  if (node == NULL)
1399        return NULL;
1400 
1401  // search nodes
1402  std::map<KdNode *, KdIntersectable *>::
1403        const_iterator it = mKdIntersectables.find(node);
1404 
1405  if (it != mKdIntersectables.end())
1406        {
1407          return (*it).second;
1408        }
1409 
1410  // not in map => create new entry
1411  KdIntersectable *kdObj = new KdIntersectable(node);
1412  mKdIntersectables[node] = kdObj;
1413 
1414  return kdObj;
[1387]1415}
[1594]1416
1417KdNode *
1418KdTree::GetNode(const Vector3 &point,
1419                                const float maxArea) const
1420{
1421  KdNode *node = mRoot;
1422 
1423  while (!node->IsLeaf() && (GetSurfaceArea(node) > maxArea) ) {
1424        KdInterior *inter = (KdInterior *)node;
1425        if (point[inter->mAxis] < inter->mPosition)
1426          node = inter->mBack;
1427        else
1428          node = inter->mFront;
1429  }
1430 
1431  return node;
1432}
1433
1434
1435}
Note: See TracBrowser for help on using the repository browser.