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

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