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

Revision 2176, 35.8 KB checked in by mattausch, 18 years ago (diff)

removed using namespace std from .h

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