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

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