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

Revision 1867, 33.7 KB checked in by bittner, 18 years ago (diff)

merge, global lines, rss sampling updates

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