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

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