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

Revision 2332, 34.7 KB checked in by mattausch, 18 years ago (diff)

implemented part of rendering estimation of wimmer et al. for view space / object space subdivision.
warning: not working with undersampling estimation + local visibility based subdivision.

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