source: GTP/trunk/Lib/Vis/Preprocessing/src/TraversalTree.cpp @ 2073

Revision 2073, 17.0 KB checked in by mattausch, 17 years ago (diff)
Line 
1#include <stack>
2#include <algorithm>
3#include <queue>
4#include "Environment.h"
5#include "Mesh.h"
6#include "TraversalTree.h"
7#include "ViewCell.h"
8#include "Beam.h"
9
10
11// $$JB HACK
12#define KD_PVS_AREA (1e-5f)
13
14namespace GtpVisibilityPreprocessor {
15
16int TraversalNode::sMailId = 1;
17int TraversalNode::sReservedMailboxes = 1;
18
19
20inline static bool ilt(Intersectable *obj1, Intersectable *obj2)
21{
22        return obj1->mId < obj2->mId;
23}
24
25
26TraversalNode::TraversalNode(TraversalInterior *parent):mParent(parent), mMailbox(0), mIntersectable(NULL)
27
28{
29  if (parent)
30    mDepth = parent->mDepth+1;
31  else
32    mDepth = 0;
33}
34
35
36TraversalInterior::~TraversalInterior()
37{
38        // recursivly destroy children
39        DEL_PTR(mFront);
40        DEL_PTR(mBack);
41}
42
43
44TraversalLeaf::~TraversalLeaf()
45{
46        DEL_PTR(mViewCell); 
47}
48
49
50TraversalTree::TraversalTree()
51{
52
53 
54  mRoot = new TraversalLeaf(NULL, 0);
55  Environment::GetSingleton()->GetIntValue("TraversalTree.Termination.maxNodes", mTermMaxNodes);
56  Environment::GetSingleton()->GetIntValue("TraversalTree.Termination.maxDepth", mTermMaxDepth);
57  Environment::GetSingleton()->GetIntValue("TraversalTree.Termination.minCost", mTermMinCost);
58  Environment::GetSingleton()->GetFloatValue("TraversalTree.Termination.maxCostRatio", mMaxCostRatio);
59  Environment::GetSingleton()->GetFloatValue("TraversalTree.Termination.ct_div_ci", mCt_div_ci);
60  Environment::GetSingleton()->GetFloatValue("TraversalTree.splitBorder", mSplitBorder);
61
62  Environment::GetSingleton()->GetBoolValue("TraversalTree.sahUseFaces", mSahUseFaces);
63
64  char splitType[64];
65  Environment::GetSingleton()->GetStringValue("TraversalTree.splitMethod", splitType);
66 
67  mSplitMethod = SPLIT_SPATIAL_MEDIAN;
68  if (strcmp(splitType, "spatialMedian") == 0)
69    mSplitMethod = SPLIT_SPATIAL_MEDIAN;
70  else
71    if (strcmp(splitType, "objectMedian") == 0)
72      mSplitMethod = SPLIT_OBJECT_MEDIAN;
73    else
74      if (strcmp(splitType, "SAH") == 0)
75        mSplitMethod = SPLIT_SAH;
76      else {
77        cerr<<"Wrong kd split type "<<splitType<<endl;
78        exit(1);
79      }
80  splitCandidates = NULL;
81}
82
83
84TraversalTree::~TraversalTree()
85{
86    DEL_PTR(mRoot);
87
88        CLEAR_CONTAINER(mKdIntersectables);
89}
90
91
92bool
93TraversalTree::Construct()
94{
95
96  if (!splitCandidates)
97    splitCandidates = new vector<SortableEntry *>;
98
99  // first construct a leaf that will get subdivide
100  TraversalLeaf *leaf = (TraversalLeaf *) mRoot;
101
102  mStat.nodes = 1;
103
104  mBox.Initialize();
105  ObjectContainer::const_iterator mi;
106  for ( mi = leaf->mObjects.begin();
107        mi != leaf->mObjects.end();
108        mi++) {
109        //      cout<<(*mi)->GetBox()<<endl;
110        mBox.Include((*mi)->GetBox());
111  }
112
113  cout <<"TraversalTree Root Box:"<<mBox<<endl;
114  mRoot = Subdivide(TraversalData(leaf, mBox, 0));
115
116  // remove the allocated array
117  CLEAR_CONTAINER(*splitCandidates);
118  delete splitCandidates;
119
120  float area = GetBox().SurfaceArea()*KD_PVS_AREA;
121
122  SetPvsTerminationNodes(area);
123 
124  return true;
125}
126
127TraversalNode *
128TraversalTree::Subdivide(const TraversalData &tdata)
129{
130
131  TraversalNode *result = NULL;
132
133  priority_queue<TraversalData> tStack;
134  //  stack<STraversalData> tStack;
135 
136  tStack.push(tdata);
137  AxisAlignedBox3 backBox, frontBox;
138
139  while (!tStack.empty()) {
140        //      cout<<mStat.Nodes()<<" "<<mTermMaxNodes<<endl;
141        if (mStat.Nodes() > mTermMaxNodes) {
142          //    if ( GetMemUsage() > maxMemory ) {
143      // count statistics on unprocessed leafs
144      while (!tStack.empty()) {
145                EvaluateLeafStats(tStack.top());
146                tStack.pop();
147      }
148      break;
149    }
150         
151       
152    TraversalData data = tStack.top();
153    tStack.pop();
154   
155    TraversalNode *node = SubdivideNode((TraversalLeaf *) data.mNode,
156                                                                 data.mBox,
157                                                                 backBox,
158                                                                 frontBox
159                                                                 );
160
161        if (result == NULL)
162      result = node;
163   
164    if (!node->IsLeaf()) {
165
166      TraversalInterior *interior = (TraversalInterior *) node;
167      // push the children on the stack
168      tStack.push(TraversalData(interior->mBack, backBox, data.mDepth+1));
169      tStack.push(TraversalData(interior->mFront, frontBox, data.mDepth+1));
170     
171    } else {
172      EvaluateLeafStats(data);
173    }
174  }
175 
176  return result;
177
178}
179
180
181bool
182TraversalTree::TerminationCriteriaMet(const TraversalLeaf *leaf)
183{
184        const bool criteriaMet =
185                ((int)leaf->mObjects.size() <= mTermMinCost) ||
186                 (leaf->mDepth >= mTermMaxDepth);
187 
188        if (criteriaMet)
189                cerr<<"\n OBJECTS="<<leaf->mObjects.size()<<endl;
190
191        return criteriaMet;
192}
193
194
195int
196TraversalTree::SelectPlane(TraversalLeaf *leaf,
197                                        const AxisAlignedBox3 &box,
198                                        float &position
199                                        )
200{
201  int axis = -1;
202 
203  switch (mSplitMethod)
204    {
205    case SPLIT_SPATIAL_MEDIAN: {
206      axis = box.Size().DrivingAxis();
207      position = (box.Min()[axis] + box.Max()[axis])*0.5f;
208      break;
209    }
210    case SPLIT_SAH: {
211      int objectsBack, objectsFront;
212      float costRatio;
213      bool mOnlyDrivingAxis = true;
214
215          if (mOnlyDrivingAxis) {
216                axis = box.Size().DrivingAxis();
217                costRatio = BestCostRatio(leaf,
218                                                                  box,
219                                                                  axis,
220                                                                  position,
221                                                                  objectsBack,
222                                                                  objectsFront);
223      } else {
224                costRatio = MAX_FLOAT;
225                for (int i=0; i < 3; i++) {
226                  float p;
227                  float r = BestCostRatio(leaf,
228                                                                  box,
229                                                                  i,
230                                                                  p,
231                                                                  objectsBack,
232                                                                  objectsFront);
233                  if (r < costRatio) {
234                        costRatio = r;
235                        axis = i;
236                        position = p;
237                  }
238                }
239      }
240     
241      if (costRatio > mMaxCostRatio) {
242                //cout<<"Too big cost ratio "<<costRatio<<endl;
243                axis = -1;
244      }
245      break;
246    }
247   
248    }
249  return axis;
250}
251
252TraversalNode *TraversalTree::SubdivideNode(TraversalLeaf *leaf,
253                                                                                        const AxisAlignedBox3 &box,
254                                                                                        AxisAlignedBox3 &backBBox,
255                                                                                        AxisAlignedBox3 &frontBBox
256                                                                                        )
257{
258
259        if (TerminationCriteriaMet(leaf))
260                return leaf;
261
262        float position;
263
264        // select subdivision axis
265        int axis = SelectPlane( leaf, box, position );
266
267        if (axis == -1) {
268                return leaf;
269        }
270
271        mStat.nodes+=2;
272        mStat.splits[axis]++;
273
274        // add the new nodes to the tree
275        TraversalInterior *node = new TraversalInterior(leaf->mParent);
276
277        node->mAxis = axis;
278        node->mPosition = position;
279        node->mBox = box;
280
281        backBBox = box;
282        frontBBox = box;
283
284        // first count ray sides
285        int objectsBack = 0;
286        int objectsFront = 0;
287
288        backBBox.SetMax(axis, position);
289        frontBBox.SetMin(axis, position);
290
291        ObjectContainer::const_iterator mi;
292
293        for ( mi = leaf->mObjects.begin();
294                mi != leaf->mObjects.end();
295                mi++)
296        {
297                // determine the side of this ray with respect to the plane
298                AxisAlignedBox3 box = (*mi)->GetBox();
299                if (box.Max(axis) > position )
300                        objectsFront++;
301
302                if (box.Min(axis) < position )
303                        objectsBack++;
304        }
305
306
307        TraversalLeaf *back = new TraversalLeaf(node, objectsBack);
308        TraversalLeaf *front = new TraversalLeaf(node, objectsFront);
309
310
311        // replace a link from node's parent
312        if (  leaf->mParent )
313                leaf->mParent->ReplaceChildLink(leaf, node);
314
315        // and setup child links
316        node->SetupChildLinks(back, front);
317
318        for (mi = leaf->mObjects.begin(); mi != leaf->mObjects.end(); ++ mi)
319        {
320                // determine the side of this ray with respect to the plane
321                AxisAlignedBox3 box = (*mi)->GetBox();
322
323                // matt: no more ref
324                // for handling multiple objects: keep track of references
325                //if (leaf->IsRoot())
326                //      (*mi)->mReferences = 1; // initialise references at root
327
328                // matt: no more ref
329                //-- (*mi)->mReferences; // remove parent ref
330
331
332                if (box.Max(axis) >= position )
333                {
334                        front->mObjects.push_back(*mi);
335                        //++ (*mi)->mReferences;
336                }
337
338                if (box.Min(axis) < position )
339                {
340                        back->mObjects.push_back(*mi);
341                        // matt: no more ref
342                        //  ++ (*mi)->mReferences;
343                }
344
345                mStat.objectRefs -= (int)leaf->mObjects.size();
346                mStat.objectRefs += objectsBack + objectsFront;
347        }
348
349        // store objects referenced in more than one leaf
350        // for easy access
351        ProcessMultipleRefs(back);
352        ProcessMultipleRefs(front);
353
354        delete leaf;
355        return node;
356}
357
358
359void
360TraversalTreeStatistics::Print(ostream &app) const
361{
362  app << "===== TraversalTree statistics ===============\n";
363
364  app << "#N_NODES ( Number of nodes )\n" << nodes << "\n";
365
366  app << "#N_LEAVES ( Number of leaves )\n" << Leaves() << "\n";
367
368  app << "#N_SPLITS ( Number of splits in axes x y z dx dy dz)\n";
369  for (int i=0; i<7; i++)
370    app << splits[i] <<" ";
371  app <<endl;
372
373  app << "#N_RAYREFS ( Number of rayRefs )\n" <<
374    rayRefs << "\n";
375
376  app << "#N_RAYRAYREFS  ( Number of rayRefs / ray )\n" <<
377    rayRefs/(double)rays << "\n";
378
379  app << "#N_LEAFRAYREFS  ( Number of rayRefs / leaf )\n" <<
380    rayRefs/(double)Leaves() << "\n";
381
382  app << "#N_MAXOBJECTREFS  ( Max number of object refs / leaf )\n" <<
383    maxObjectRefs << "\n";
384
385  app << "#N_NONEMPTYRAYREFS  ( Number of rayRefs in nonEmpty leaves / non empty leaf )\n" <<
386    rayRefsNonZeroQuery/(double)(Leaves() - zeroQueryNodes) << "\n";
387
388  app << "#N_LEAFDOMAINREFS  ( Number of query domain Refs / leaf )\n" <<
389    objectRefs/(double)Leaves() << "\n";
390
391  //  app << setprecision(4);
392
393  app << "#N_PEMPTYLEAVES  ( Percentage of leaves with zero query domains )\n"<<
394    zeroQueryNodes*100/(double)Leaves()<<endl;
395
396  app << "#N_PMAXDEPTHLEAVES ( Percentage of leaves at maxdepth )\n"<<
397    maxDepthNodes*100/(double)Leaves()<<endl;
398
399  app << "#N_PMINCOSTLEAVES  ( Percentage of leaves with minCost )\n"<<
400    minCostNodes*100/(double)Leaves()<<endl;
401
402  app << "#N_ADDED_RAYREFS  (Number of dynamically added ray references )\n"<<
403    addedRayRefs<<endl;
404
405  app << "#N_REMOVED_RAYREFS  (Number of dynamically removed ray references )\n"<<
406    removedRayRefs<<endl;
407
408  //  app << setprecision(4);
409
410  //  app << "#N_CTIME  ( Construction time [s] )\n"
411  //      << Time() << " \n";
412
413  app << "===== END OF TraversalTree statistics ==========\n";
414
415}
416
417
418
419void
420TraversalTree::EvaluateLeafStats(const TraversalData &data)
421{
422
423  // the node became a leaf -> evaluate stats for leafs
424  TraversalLeaf *leaf = (TraversalLeaf *)data.mNode;
425
426  if (data.mDepth > mTermMaxDepth)
427    mStat.maxDepthNodes++;
428 
429  if ( (int)(leaf->mObjects.size()) < mTermMinCost)
430    mStat.minCostNodes++;
431 
432 
433  if ( (int)(leaf->mObjects.size()) > mStat.maxObjectRefs)
434    mStat.maxObjectRefs = (int)leaf->mObjects.size();
435 
436}
437
438
439
440void
441TraversalTree::SortSubdivisionCandidates(
442                            TraversalLeaf *node,
443                            const int axis
444                            )
445{
446        CLEAR_CONTAINER(*splitCandidates);
447        //splitCandidates->clear();
448
449    int requestedSize = 2*(int)node->mObjects.size();
450       
451        // creates a sorted split candidates array
452        if (splitCandidates->capacity() > 500000 &&
453                requestedSize < (int)(splitCandidates->capacity()/10) ) {               
454                        delete splitCandidates;
455                        splitCandidates = new vector<SortableEntry *>;
456  }
457 
458  splitCandidates->reserve(requestedSize);
459 
460  // insert all queries
461  for(ObjectContainer::const_iterator mi = node->mObjects.begin();
462      mi != node->mObjects.end();
463      mi++)
464  {
465          AxisAlignedBox3 box = (*mi)->GetBox();
466
467          splitCandidates->push_back(new SortableEntry(SortableEntry::BOX_MIN,
468                                                                 box.Min(axis),
469                                                                 *mi)
470                                                                 );
471   
472      splitCandidates->push_back(new SortableEntry(SortableEntry::BOX_MAX,
473                                                                 box.Max(axis),
474                                                                 *mi)
475                                                                 );
476  }
477 
478  stable_sort(splitCandidates->begin(), splitCandidates->end(), iltS);
479}
480
481
482float
483TraversalTree::BestCostRatio(
484                                          TraversalLeaf *node,
485                                          const AxisAlignedBox3 &box,
486                                          const int axis,
487                                          float &position,
488                                          int &objectsBack,
489                                          int &objectsFront
490                                          )
491{
492
493#define DEBUG_COST 0
494
495#if DEBUG_COST
496  static int nodeId = -1;
497  char filename[256];
498 
499  static int lastAxis = 100;
500  if (axis <= lastAxis)
501        nodeId++;
502
503  lastAxis = axis;
504 
505  sprintf(filename, "sah-cost%d-%d.log", nodeId, axis);
506  ofstream costStream;
507 
508  if (nodeId < 100)
509        costStream.open(filename);
510
511#endif
512 
513  SortSubdivisionCandidates(node, axis);
514 
515  // go through the lists, count the number of objects left and right
516  // and evaluate the following cost funcion:
517  // C = ct_div_ci  + (ol + or)/queries
518
519  float totalIntersections = 0.0f;
520  vector<SortableEntry *>::const_iterator ci;
521
522  for(ci = splitCandidates->begin();
523      ci < splitCandidates->end();
524      ci++)
525    if ((*ci)->type == SortableEntry::BOX_MIN) {
526      totalIntersections += (*ci)->intersectable->IntersectionComplexity();
527    }
528       
529  float intersectionsLeft = 0;
530  float intersectionsRight = totalIntersections;
531       
532  int objectsLeft = 0, objectsRight = (int)node->mObjects.size();
533 
534  float minBox = box.Min(axis);
535  float maxBox = box.Max(axis);
536  float boxArea = box.SurfaceArea();
537 
538  float minBand = minBox + mSplitBorder*(maxBox - minBox);
539  float maxBand = minBox + (1.0f - mSplitBorder)*(maxBox - minBox);
540 
541  float minSum = 1e20f;
542 
543  for(ci = splitCandidates->begin();
544      ci < splitCandidates->end();
545      ci++) {
546    switch ((*ci)->type) {
547    case SortableEntry::BOX_MIN:
548      objectsLeft++;
549      intersectionsLeft += (*ci)->intersectable->IntersectionComplexity();
550      break;
551    case SortableEntry::BOX_MAX:
552      objectsRight--;
553      intersectionsRight -= (*ci)->intersectable->IntersectionComplexity();
554      break;
555    }
556
557    if ((*ci)->value > minBand && (*ci)->value < maxBand) {
558      AxisAlignedBox3 lbox = box;
559      AxisAlignedBox3 rbox = box;
560      lbox.SetMax(axis, (*ci)->value);
561      rbox.SetMin(axis, (*ci)->value);
562
563      float sum;
564      if (mSahUseFaces)
565                sum = intersectionsLeft*lbox.SurfaceArea() + intersectionsRight*rbox.SurfaceArea();
566      else
567                sum = objectsLeft*lbox.SurfaceArea() + objectsRight*rbox.SurfaceArea();
568     
569      //      cout<<"pos="<<(*ci).value<<"\t q=("<<ql<<","<<qr<<")\t r=("<<rl<<","<<rr<<")"<<endl;
570      //      cout<<"cost= "<<sum<<endl;
571
572#if DEBUG_COST
573  if (nodeId < 100) {
574        float oldCost = mSahUseFaces ? totalIntersections : node->mObjects.size();
575        float newCost = mCt_div_ci + sum/boxArea;
576        float ratio = newCost/oldCost;
577        costStream<<(*ci)->value<<" "<<ratio<<endl;
578  }
579#endif
580         
581      if (sum < minSum) {
582                minSum = sum;
583                position = (*ci)->value;
584               
585                objectsBack = objectsLeft;
586                objectsFront = objectsRight;
587      }
588    }
589  }
590 
591  float oldCost = mSahUseFaces ? totalIntersections : node->mObjects.size();
592  float newCost = mCt_div_ci + minSum/boxArea;
593  float ratio = newCost/oldCost;
594 
595#if 0
596  cout<<"===================="<<endl;
597  cout<<"costRatio="<<ratio<<" pos="<<position<<" t="<<(position - minBox)/(maxBox - minBox)
598      <<"\t o=("<<objectsBack<<","<<objectsFront<<")"<<endl;
599#endif
600  return ratio;
601}
602
603
604int TraversalTree::CastLineSegment(const Vector3 &origin,
605                                                        const Vector3 &termination,
606                                                        ViewCellContainer &viewcells)
607{
608        int hits = 0;
609
610        float mint = 0.0f, maxt = 1.0f;
611        const Vector3 dir = termination - origin;
612
613        stack<RayTraversalData> tStack;
614
615        Intersectable::NewMail();
616
617        //maxt += Limits::Threshold;
618
619        Vector3 entp = origin;
620        Vector3 extp = termination;
621
622        TraversalNode *node = mRoot;
623        TraversalNode *farChild;
624
625        float position;
626        int axis;
627
628        while (1)
629        {
630                if (!node->IsLeaf())
631                {
632                        TraversalInterior *in = static_cast<TraversalInterior *>(node);
633                        position = in->mPosition;
634                        axis = in->mAxis;
635
636                        if (entp[axis] <= position)
637                        {
638                                if (extp[axis] <= position)
639                                {
640                                        node = in->mBack;
641                                        // cases N1,N2,N3,P5,Z2,Z3
642                                        continue;
643                                }
644                                else
645                                {
646                                        // case N4
647                                        node = in->mBack;
648                                        farChild = in->mFront;
649                                }
650                        }
651                        else
652                        {
653                                if (position <= extp[axis])
654                                {
655                                        node = in->mFront;
656                                        // cases P1,P2,P3,N5,Z1
657                                        continue;
658                                }
659                                else
660                                {
661                                        node = in->mFront;
662                                        farChild = in->mBack;
663                                        // case P4
664                                }
665                        }
666
667                        // $$ modification 3.5.2004 - hints from Kamil Ghais
668                        // case N4 or P4
669                        float tdist = (position - origin[axis]) / dir[axis];
670                        //tStack.push(RayTraversalData(farChild, extp, maxt)); //TODO
671                        extp = origin + dir * tdist;
672                        maxt = tdist;
673                }
674                else
675                {
676                        // compute intersection with all objects in this leaf
677                        TraversalLeaf *leaf = static_cast<TraversalLeaf *>(node);
678
679                        // add view cell to intersections
680                        ViewCell *vc = leaf->mViewCell;
681
682                        if (!vc->Mailed())
683                        {
684                                vc->Mail();
685                                viewcells.push_back(vc);
686                                ++ hits;
687                        }
688
689                        // get the next node from the stack
690                        if (tStack.empty())
691                                break;
692
693                        entp = extp;
694                        mint = maxt;
695                       
696                        RayTraversalData &s  = tStack.top();
697                        node = s.mNode;
698                        extp = s.mExitPoint;
699                        maxt = s.mMaxT;
700                        tStack.pop();
701                }
702        }
703
704        return hits;
705}
706
707
708void
709TraversalTree::CollectLeaves(vector<TraversalLeaf *> &leaves)
710{
711  stack<TraversalNode *> nodeStack;
712  nodeStack.push(mRoot);
713
714  while (!nodeStack.empty()) {
715    TraversalNode *node = nodeStack.top();
716    nodeStack.pop();
717    if (node->IsLeaf()) {
718      TraversalLeaf *leaf = (TraversalLeaf *)node;
719      leaves.push_back(leaf);
720    } else {
721      TraversalInterior *interior = (TraversalInterior *)node;
722      nodeStack.push(interior->mBack);
723      nodeStack.push(interior->mFront);
724    }
725  }
726}
727
728
729}
Note: See TracBrowser for help on using the repository browser.