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

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