source: trunk/VUT/GtpVisibilityPreprocessor/src/VspKdTree.cpp @ 410

Revision 410, 34.5 KB checked in by mattausch, 19 years ago (diff)
Line 
1
2// ================================================================
3// $Id: lsds_kdtree.cpp,v 1.18 2005/04/16 09:34:21 bittner Exp $
4// ****************************************************************
5/**
6   The KD tree based LSDS
7*/
8// Initial coding by
9/**
10   @author Jiri Bittner
11*/
12
13// Standard headers
14#include <stack>
15#include <queue>
16#include <algorithm>
17#include <fstream>
18#include <string>
19
20#include "VspKdTree.h"
21
22#include "Environment.h"
23#include "VssRay.h"
24#include "Intersectable.h"
25#include "Ray.h"
26
27
28#define DEBUG_DIR_SPLIT 0
29
30
31
32// Static variables
33int
34VspKdTreeLeaf::mailID = 0;
35
36inline void AddObject2Pvs(Intersectable *object,
37                                                  const int side,
38                                                  int &pvsBack,
39                                                  int &pvsFront)
40{
41        if (!object)
42                return;
43               
44        if (side <= 0)
45        {
46                if (!object->Mailed() && !object->Mailed(2))
47                {
48                        ++ pvsBack;
49
50                        if (object->Mailed(1))
51                                object->Mail(2);
52                        else
53                                object->Mail();
54                }
55        }
56       
57        if (side >= 0)
58        {
59                if (!object->Mailed(1) && !object->Mailed(2))
60                {
61                        ++ pvsFront;
62                        if (object->Mailed())
63                                object->Mail(2);
64                        else
65                                object->Mail(1);
66                }
67        }
68}
69
70
71// Constructor
72VspKdTree::VspKdTree()
73{
74        environment->GetIntValue("VspKdTree.maxDepth", termMaxDepth);
75        environment->GetIntValue("VspKdTree.minPvs", termMinPvs);
76        environment->GetIntValue("VspKdTree.minRays", termMinRays);
77        environment->GetFloatValue("VspKdTree.maxRayContribution", termMaxRayContribution);
78        environment->GetFloatValue("VspKdTree.maxCostRatio", termMaxCostRatio);
79
80        environment->GetFloatValue("VspKdTree.minSize", termMinSize);
81        termMinSize = sqr(termMinSize);
82       
83        environment->GetFloatValue("VspKdTree.refDirBoxMaxSize", refDirBoxMaxSize);
84        refDirBoxMaxSize = sqr(refDirBoxMaxSize);
85 
86        environment->GetFloatValue("VspKdTree.epsilon", epsilon);
87        environment->GetFloatValue("VspKdTree.ct_div_ci", ct_div_ci);
88       
89        environment->GetFloatValue("VspKdTree.maxTotalMemory", maxTotalMemory);
90        environment->GetFloatValue("VspKdTree.maxStaticMemory", maxStaticMemory);
91 
92 
93        float refDirAngle;
94        environment->GetFloatValue("VspKdTree.refDirAngle", refDirAngle);
95
96        environment->GetIntValue("VspKdTree.accessTimeThreshold", accessTimeThreshold);
97        //= 1000;
98        environment->GetIntValue("VspKdTree.minCollapseDepth", minCollapseDepth);
99 
100
101        // split type
102        char sname[128];
103        environment->GetStringValue("VspKdTree.splitType", sname);
104        string name(sname);
105       
106        if (name.compare("regular") == 0)
107                splitType = ESplitRegular;
108        else
109        {
110                if (name.compare("heuristic") == 0)
111                        splitType = ESplitHeuristic;
112                else
113                {
114                        cerr << "Invalid VspKdTree split type " << name << endl;
115                        exit(1);
116                }
117        }
118
119        environment->GetBoolValue("VspKdTree.randomize", randomize);
120
121        root = NULL;
122
123        splitCandidates = new vector<SortableEntry>;
124}
125
126
127VspKdTree::~VspKdTree()
128{
129        if (root)
130                delete root;
131}
132
133
134
135
136void
137VspKdStatistics::Print(ostream &app) const
138{
139  app << "===== VspKdTree statistics ===============\n";
140
141  app << "#N_RAYS ( Number of rays )\n"
142      << rays <<endl;
143 
144  app << "#N_INITPVS ( Initial PVS size )\n"
145      << initialPvsSize <<endl;
146 
147  app << "#N_NODES ( Number of nodes )\n" << nodes << "\n";
148
149  app << "#N_LEAVES ( Number of leaves )\n" << Leaves() << "\n";
150
151  app << "#N_SPLITS ( Number of splits in axes x y z dx dy dz \n";
152  for (int i=0; i<7; i++)
153    app << splits[i] <<" ";
154  app <<endl;
155
156  app << "#N_RAYREFS ( Number of rayRefs )\n" <<
157    rayRefs << "\n";
158
159  app << "#N_RAYRAYREFS  ( Number of rayRefs / ray )\n" <<
160    rayRefs/(double)rays << "\n";
161
162  app << "#N_LEAFRAYREFS  ( Number of rayRefs / leaf )\n" <<
163    rayRefs/(double)Leaves() << "\n";
164
165  app << "#N_MAXRAYREFS  ( Max number of rayRefs / leaf )\n" <<
166    maxRayRefs << "\n";
167
168
169  //  app << setprecision(4);
170
171  app << "#N_PMAXDEPTHLEAVES ( Percentage of leaves at maxdepth )\n"<<
172    maxDepthNodes*100/(double)Leaves()<<endl;
173
174  app << "#N_PMINPVSLEAVES  ( Percentage of leaves with minCost )\n"<<
175    minPvsNodes*100/(double)Leaves()<<endl;
176
177  app << "#N_PMINRAYSLEAVES  ( Percentage of leaves with minCost )\n"<<
178    minRaysNodes*100/(double)Leaves()<<endl;
179       
180        app << "#N_PMINSIZELEAVES  ( Percentage of leaves with minSize )\n"<<
181    minSizeNodes*100/(double)Leaves()<<endl;
182
183        app << "#N_PMAXRAYCONTRIBLEAVES  ( Percentage of leaves with maximal ray contribution )\n"<<
184    maxRayContribNodes*100/(double)Leaves()<<endl;
185
186  app << "#N_ADDED_RAYREFS  (Number of dynamically added ray references )\n"<<
187    addedRayRefs<<endl;
188
189  app << "#N_REMOVED_RAYREFS  (Number of dynamically removed ray references )\n"<<
190    removedRayRefs<<endl;
191
192  //  app << setprecision(4);
193
194  app << "#N_CTIME  ( Construction time [s] )\n"
195      << Time() << " \n";
196
197  app << "===== END OF VspKdTree statistics ==========\n";
198
199}
200
201
202void
203VspKdTreeLeaf::UpdatePvsSize()
204{
205        if (!mValidPvs)
206        {
207                Intersectable::NewMail();
208                int pvsSize = 0;
209                for(VspKdTreeNode::RayInfoContainer::iterator ri = rays.begin();
210                        ri != rays.end(); ++ ri)
211                {
212                        if ((*ri).mRay->IsActive())
213                        {
214                                Intersectable *object;
215#if BIDIRECTIONAL_RAY
216                                object = (*ri).mRay->mOriginObject;
217                       
218                                if (object && !object->Mailed())
219                                {
220                                        ++ pvsSize;
221                                        object->Mail();
222                                }
223#endif
224                                object = (*ri).mRay->mTerminationObject;
225                                if (object && !object->Mailed())
226                                {
227                                        ++ pvsSize;
228                                        object->Mail();
229                                }
230                        }
231                }
232
233                mPvsSize = pvsSize;
234                mValidPvs = true;
235        }
236}
237
238
239void
240VspKdTree::Construct(VssRayContainer &rays,
241                                         AxisAlignedBox3 *forcedBoundingBox)
242{
243        stat.Start();
244 
245        maxMemory = maxStaticMemory;
246
247        if (root)
248                delete root;
249
250        root = new VspKdTreeLeaf(NULL, rays.size());
251
252        // first construct a leaf that will get subdivide
253        VspKdTreeLeaf *leaf = (VspKdTreeLeaf *) root;
254
255       
256        stat.nodes = 1;
257       
258        bbox.Initialize();
259        dirBBox.Initialize();
260 
261       
262        for (VssRayContainer::const_iterator ri = rays.begin(); ri != rays.end(); ++ ri)
263        {
264                leaf->AddRay(VspKdTreeNode::RayInfo(*ri));
265
266                bbox.Include((*ri)->GetOrigin());
267                bbox.Include((*ri)->GetTermination());   
268               
269                dirBBox.Include(Vector3((*ri)->GetDirParametrization(0), (*ri)->GetDirParametrization(1), 0));
270        }
271       
272        if (forcedBoundingBox)
273                bbox = *forcedBoundingBox;
274
275        cout<<"Bbox = "<<bbox<<endl;
276        cout<<"Dirr Bbox = "<<dirBBox<<endl;
277
278        stat.rays = leaf->rays.size();
279        leaf->UpdatePvsSize();
280       
281        stat.initialPvsSize = leaf->GetPvsSize();
282       
283        // Subdivide();
284        root = Subdivide(TraversalData(leaf, bbox, 0));
285
286        if (splitCandidates)
287        {
288                // force realease of this vector
289                delete splitCandidates;
290                splitCandidates = new vector<SortableEntry>;
291        }
292 
293        stat.Stop();
294
295        stat.Print(cout);
296        cout<<"#Total memory=" << GetMemUsage() << endl;
297}
298
299
300
301VspKdTreeNode *VspKdTree::Subdivide(const TraversalData &tdata)
302{
303        VspKdTreeNode *result = NULL;
304
305        priority_queue<TraversalData> tStack;
306        //  stack<TraversalData> tStack;
307 
308        tStack.push(tdata);
309
310        AxisAlignedBox3 backBox;
311        AxisAlignedBox3 frontBox;
312 
313        int lastMem = 0;
314        while (!tStack.empty())
315        {
316                float mem = GetMemUsage();
317               
318                if ( lastMem/10 != ((int)mem)/10)
319                {
320                        cout << mem << " MB" << endl;
321                }
322                lastMem = (int)mem;
323               
324                if (  mem > maxMemory )
325                {
326                        // count statistics on unprocessed leafs
327                        while (!tStack.empty())
328                        {
329                                EvaluateLeafStats(tStack.top());
330                                tStack.pop();
331                        }
332                        break;
333                }
334   
335                TraversalData data = tStack.top();
336                tStack.pop();
337   
338               
339                VspKdTreeNode *node = SubdivideNode((VspKdTreeLeaf *) data.node,
340                                                                                        data.bbox, backBox,     frontBox);
341                if (result == NULL)
342                        result = node;
343   
344                if (!node->IsLeaf())
345                {
346                        VspKdTreeInterior *interior = (VspKdTreeInterior *) node;
347
348                        // push the children on the stack
349                        tStack.push(TraversalData(interior->back, backBox, data.depth+1));
350                        tStack.push(TraversalData(interior->front, frontBox, data.depth+1));
351                }
352                else
353                {
354                        EvaluateLeafStats(data);
355                }
356        }
357
358        return result;
359}
360
361
362// returns selected plane for subdivision
363int
364VspKdTree::SelectPlane(VspKdTreeLeaf *leaf,
365                                           const AxisAlignedBox3 &box,
366                                           float &position,
367                                           int &raysBack,
368                                           int &raysFront,
369                                           int &pvsBack,
370                                           int &pvsFront)
371{
372        int minDirDepth = 6;
373        int axis;
374        float costRatio;
375       
376        if (splitType == ESplitRegular) {
377                costRatio = BestCostRatioRegular(leaf,
378                                                                                 axis,
379                                                                                 position,
380                                                                                 raysBack,
381                                                                                 raysFront,
382                                                                                 pvsBack,
383                                                                                 pvsFront);
384
385        }
386        else
387        {
388                if (splitType == ESplitHeuristic)
389                        costRatio = BestCostRatioHeuristic(leaf,
390                                                                                           axis,       
391                                                                                           position,
392                                                                                           raysBack,
393                                                                                           raysFront,
394                                                                                           pvsBack,
395                                                                                           pvsFront);
396                else {
397                        cerr << "VspKdTree: Unknown split heuristics\n";
398                        exit(1);
399                }
400  }
401
402        if (costRatio > termMaxCostRatio)
403        {
404                //              cout<<"Too big cost ratio "<<costRatio<<endl;
405                return -1;
406        }
407
408#if 0   
409        cout<<
410                "pvs="<<leaf->mPvsSize<<
411                " rays="<<leaf->rays.size()<<
412                " rc="<<leaf->GetAvgRayContribution()<<
413                " axis="<<axis<<endl;
414#endif
415       
416        return axis;
417}
418
419
420                                                       
421
422float
423VspKdTree::EvalCostRatio(VspKdTreeLeaf *leaf,
424                                                 const int axis,
425                                                 const float position,
426                                                 int &raysBack,
427                                                 int &raysFront,
428                                                 int &pvsBack,
429                                                 int &pvsFront)
430{
431        raysBack = 0;
432        raysFront = 0;
433        pvsFront = 0;
434        pvsBack = 0;
435
436        float newCost;
437
438        Intersectable::NewMail(3);
439       
440        // eval pvs size
441        int pvsSize = leaf->GetPvsSize();
442
443        Intersectable::NewMail(3);
444
445        if (axis <= VspKdTreeNode::SPLIT_Z) {
446    // this is the main ray classification loop!
447    for(VspKdTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
448                ri != leaf->rays.end(); ri++)
449                if ((*ri).mRay->IsActive())
450                {
451                        // determine the side of this ray with respect to the plane
452                        int side = (*ri).ComputeRayIntersection(axis, position, (*ri).mRay->mT);
453                        //                              (*ri).mRay->mSide = side;
454                       
455                        if (side <= 0)
456                                raysBack ++;
457                               
458                        if (side >= 0)
459                                raysFront++;
460                               
461                        AddObject2Pvs((*ri).mRay->mTerminationObject, side, pvsBack, pvsFront);
462                }
463               
464                AxisAlignedBox3 box = GetBBox(leaf);
465               
466                float minBox = box.Min(axis);
467                float maxBox = box.Max(axis);
468                float sizeBox = maxBox - minBox;
469               
470                //              float sum = raysBack*(position - minBox) + raysFront*(maxBox - position);
471                float sum = pvsBack*(position - minBox) + pvsFront*(maxBox - position);
472
473                newCost = ct_div_ci + sum/sizeBox;
474        }
475        else
476        {
477                // directional split
478                for(VspKdTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
479                        ri != leaf->rays.end(); ++ ri)
480                        if ((*ri).mRay->IsActive())
481                        {
482                                // determine the side of this ray with respect to the plane
483                                int side;
484                                if ((*ri).mRay->GetDirParametrization(axis - 3) > position)
485                                        side = 1;
486                                else
487                                        side = -1;
488
489                                if (side <= 0)
490                                        raysBack++;
491                               
492                                if (side >= 0)
493                                        raysFront++;
494
495                                //                              (*ri).mRay->mSide = side;
496                                AddObject2Pvs((*ri).mRay->mTerminationObject, side, pvsBack, pvsFront);
497                        }
498       
499                        AxisAlignedBox3 box = GetDirBBox(leaf);
500               
501                        float minBox = box.Min(axis - 3);
502                        float maxBox = box.Max(axis - 3);
503                        float sizeBox = maxBox - minBox;
504               
505                        // float sum = raysBack*(position - minBox) + raysFront*(maxBox - position);
506                        float sum =     pvsBack*(position - minBox) + pvsFront*(maxBox - position);
507                //              float sum = pvsBack + pvsFront;
508                newCost = ct_div_ci + sum/sizeBox;
509        }
510
511        //      cout<<axis<<" "<<pvsSize<<" "<<pvsBack<<" "<<pvsFront<<endl;
512        //  float oldCost = leaf->rays.size();
513        float oldCost = pvsSize;
514       
515        float ratio = newCost / oldCost;
516        //      cout<<"ratio="<<ratio<<endl;
517        return ratio;
518}
519
520float
521VspKdTree::BestCostRatioRegular(VspKdTreeLeaf *leaf,
522                                                                int &axis,
523                                                                float &position,
524                                                                int &raysBack,
525                                                                int &raysFront,
526                                                                int &pvsBack,
527                                                                int &pvsFront)
528{
529        int nRaysBack[6], nRaysFront[6];
530        int nPvsBack[6], nPvsFront[6];
531        float nPosition[6];
532        float nCostRatio[6];
533        int bestAxis = -1;
534       
535        AxisAlignedBox3 sBox = GetBBox(leaf);
536        AxisAlignedBox3 dBox = GetDirBBox(leaf);
537
538        // int sAxis = box.Size().DrivingAxis();
539        int sAxis = sBox.Size().DrivingAxis();
540        int dAxis = dBox.Size().DrivingAxis() + 3;
541
542        bool onlyDrivingAxis = true;
543
544        for (axis = 0; axis < 6; ++ axis)
545        {
546                if (!onlyDrivingAxis || axis == sAxis || axis == dAxis)
547                {
548                        if (axis < 3)
549                                nPosition[axis] = (sBox.Min()[axis] + sBox.Max()[axis])*0.5f;
550                        else
551                                nPosition[axis] = (dBox.Min()[axis-3] + dBox.Max()[axis-3])*0.5f;
552                       
553                        nCostRatio[axis] = EvalCostRatio(leaf,
554                                                                                         axis,
555                                                                                         nPosition[axis],
556                                                                                         nRaysBack[axis],
557                                                                                         nRaysFront[axis],
558                                                                                         nPvsBack[axis],
559                                                                                         nPvsFront[axis]);
560                       
561                        if (bestAxis == -1)
562                                bestAxis = axis;
563                        else
564                                if (nCostRatio[axis] < nCostRatio[bestAxis])
565                                        bestAxis = axis;
566                }
567        }
568
569        axis = bestAxis;
570        position = nPosition[bestAxis];
571
572        raysBack = nRaysBack[bestAxis];
573        raysFront = nRaysFront[bestAxis];
574
575        pvsBack = nPvsBack[bestAxis];
576        pvsFront = nPvsFront[bestAxis];
577       
578        return nCostRatio[bestAxis];
579}
580
581float VspKdTree::BestCostRatioHeuristic(VspKdTreeLeaf *leaf,
582                                                                                int &axis,
583                                                                                float &position,
584                                                                                int &raysBack,
585                                                                                int &raysFront,
586                                                                                int &pvsBack,
587                                                                                int &pvsFront)
588{
589        AxisAlignedBox3 box = GetBBox(leaf);
590        //      AxisAlignedBox3 dirBox = GetDirBBox(node);
591       
592        axis = box.Size().DrivingAxis();
593               
594        SortSplitCandidates(leaf, axis);
595 
596        // go through the lists, count the number of objects left and right
597        // and evaluate the following cost funcion:
598        // C = ct_div_ci  + (ql*rl + qr*rr)/queries
599       
600        int rl=0, rr = leaf->rays.size();
601        int pl=0, pr = leaf->GetPvsSize();
602        float minBox = box.Min(axis);
603        float maxBox = box.Max(axis);
604        float sizeBox = maxBox - minBox;
605       
606        float minBand = minBox + 0.1*(maxBox - minBox);
607        float maxBand = minBox + 0.9*(maxBox - minBox);
608       
609        float sum = rr*sizeBox;
610        float minSum = 1e20;
611       
612        Intersectable::NewMail();
613        // set all object as belonging to the fron pvs
614        for(VspKdTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
615                ri != leaf->rays.end(); ++ ri)
616        {
617                if ((*ri).mRay->IsActive())
618                {
619                        Intersectable *object = (*ri).mRay->mTerminationObject;
620                       
621                        if (object)
622                                if (!object->Mailed())
623                                {
624                                        object->Mail();
625                                        object->mCounter = 1;
626                                }
627                                else
628                                        ++ object->mCounter;
629                }
630        }
631       
632        Intersectable::NewMail();
633       
634        for (vector<SortableEntry>::const_iterator ci = splitCandidates->begin();
635                ci < splitCandidates->end(); ++ ci)
636        {
637                VssRay *ray;
638                       
639                switch ((*ci).type)
640                {
641                        case SortableEntry::ERayMin:
642                                {
643                                        ++ rl;
644                                        ray = (VssRay *) (*ci).data;
645                                        Intersectable *object = ray->mTerminationObject;
646                                        if (object && !object->Mailed())
647                                        {
648                                                object->Mail();
649                                                ++ pl;
650                                        }
651                                        break;
652                                }
653                        case SortableEntry::ERayMax:
654                                {
655                                        -- rr;
656                                       
657                                        ray = (VssRay *) (*ci).data;
658                                        Intersectable *object = ray->mTerminationObject;
659                                       
660                                        if (object)
661                                        {
662                                                if (-- object->mCounter == 0)
663                                                -- pr;
664                                        }
665                                               
666                                        break;
667                                }
668                }
669                       
670                if ((*ci).value > minBand && (*ci).value < maxBand)
671                {
672                        sum = pl*((*ci).value - minBox) + pr*(maxBox - (*ci).value);
673               
674                        //  cout<<"pos="<<(*ci).value<<"\t q=("<<ql<<","<<qr<<")\t r=("<<rl<<","<<rr<<")"<<endl;
675                        // cout<<"cost= "<<sum<<endl;
676                       
677                        if (sum < minSum)
678                        {
679                                minSum = sum;
680                                position = (*ci).value;
681                       
682                                raysBack = rl;
683                                raysFront = rr;
684                       
685                                pvsBack = pl;
686                                pvsFront = pr;
687                               
688                        }
689                }
690        }
691
692        float oldCost = leaf->GetPvsSize();
693        float newCost = ct_div_ci + minSum / sizeBox;
694        float ratio = newCost / oldCost;
695 
696        //  cout<<"===================="<<endl;
697        //  cout<<"costRatio="<<ratio<<" pos="<<position<<" t="<<(position - minBox)/(maxBox - minBox)
698        //      <<"\t q=("<<queriesBack<<","<<queriesFront<<")\t r=("<<raysBack<<","<<raysFront<<")"<<endl;
699       
700        return ratio;
701}
702
703void VspKdTree::SortSplitCandidates(VspKdTreeLeaf *node,
704                                                                        const int axis)
705{
706        splitCandidates->clear();
707 
708        int requestedSize = 2*(node->rays.size());
709        // creates a sorted split candidates array
710        if (splitCandidates->capacity() > 500000 &&
711                requestedSize < (int)(splitCandidates->capacity()/10) )
712        {
713        delete splitCandidates;
714                splitCandidates = new vector<SortableEntry>;
715        }
716 
717        splitCandidates->reserve(requestedSize);
718
719        // insert all queries
720        for(VspKdTreeNode::RayInfoContainer::const_iterator ri = node->rays.begin();
721                ri < node->rays.end(); ++ ri)
722        {
723                bool positive = (*ri).mRay->HasPosDir(axis);
724                splitCandidates->push_back(SortableEntry(positive ? SortableEntry::ERayMin : SortableEntry::ERayMax,
725                                                                                                 (*ri).ExtrapOrigin(axis), (void *)&*ri));
726               
727                splitCandidates->push_back(SortableEntry(positive ? SortableEntry::ERayMax : SortableEntry::ERayMin,
728                                                                                                 (*ri).ExtrapTermination(axis), (void *)&*ri));
729        }
730       
731        stable_sort(splitCandidates->begin(), splitCandidates->end());
732}
733
734
735void VspKdTree::EvaluateLeafStats(const TraversalData &data)
736{
737        // the node became a leaf -> evaluate stats for leafs
738        VspKdTreeLeaf *leaf = (VspKdTreeLeaf *)data.node;
739
740        if (data.depth >= termMaxDepth)
741                ++ stat.maxDepthNodes;
742 
743        //  if ( (int)(leaf->rays.size()) < termMinCost)
744        //    stat.minCostNodes++;
745        if (leaf->GetPvsSize() < termMinPvs)
746                ++ stat.minPvsNodes;
747
748        if (leaf->GetPvsSize() < termMinRays)
749                ++ stat.minRaysNodes;
750
751        if (0 && leaf->GetAvgRayContribution() > termMaxRayContribution)
752                ++ stat.maxRayContribNodes;
753       
754        if (SqrMagnitude(data.bbox.Size()) <= termMinSize)
755                ++ stat.minSizeNodes;
756
757        if ((int)(leaf->rays.size()) > stat.maxRayRefs)
758                stat.maxRayRefs = leaf->rays.size();
759}
760
761
762
763VspKdTreeNode *VspKdTree::SubdivideNode(VspKdTreeLeaf *leaf,
764                                                                                const AxisAlignedBox3 &box,
765                                                                                AxisAlignedBox3 &backBBox,
766                                                                                AxisAlignedBox3 &frontBBox)
767{
768        if ( (leaf->GetPvsSize() < termMinPvs) || (leaf->rays.size() < termMinRays) ||
769        // (leaf->GetAvgRayContribution() > termMaxRayContribution ) ||
770       (leaf->depth >= termMaxDepth) || SqrMagnitude(box.Size()) <= termMinSize )
771        {
772
773#if 0
774                if (leaf->depth >= termMaxDepth) {
775                        cout<<"Warning: max depth reached depth="<<(int)leaf->depth<<" rays="<<leaf->rays.size()<<endl;
776                        cout<<"Bbox: "<<GetBBox(leaf)<<" dirbbox:"<<GetDirBBox(leaf)<<endl;
777                }
778#endif
779               
780                return leaf;
781        }
782       
783        float position;
784        // first count ray sides
785        int raysBack;
786        int raysFront;
787        int pvsBack;
788        int pvsFront;
789       
790        // select subdivision axis
791        int axis = SelectPlane( leaf, box, position, raysBack, raysFront, pvsBack, pvsFront);
792
793        //      cout<<"rays back="<<raysBack<<" rays front="<<raysFront<<" pvs back="<<pvsBack<<" pvs front="<<
794        //              pvsFront<<endl;
795
796        if (axis == -1)
797        {
798                return leaf;
799        }
800 
801        stat.nodes += 2;
802        ++ stat.splits[axis];
803
804        // add the new nodes to the tree
805        VspKdTreeInterior *node = new VspKdTreeInterior(leaf->parent);
806
807        node->axis = axis;
808        node->position = position;
809        node->bbox = box;
810        node->dirBBox = GetDirBBox(leaf);
811 
812        backBBox = box;
813        frontBBox = box;
814
815        VspKdTreeLeaf *back = new VspKdTreeLeaf(node, raysBack);
816        back->SetPvsSize(pvsBack);
817        VspKdTreeLeaf *front = new VspKdTreeLeaf(node, raysFront);
818        front->SetPvsSize(pvsFront);
819
820        // replace a link from node's parent
821        if (leaf->parent)
822                leaf->parent->ReplaceChildLink(leaf, node);
823        // and setup child links
824        node->SetupChildLinks(back, front);
825       
826        if (axis <= VspKdTreeNode::SPLIT_Z)
827        {
828                backBBox.SetMax(axis, position);
829                frontBBox.SetMin(axis, position);
830               
831                for(VspKdTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
832                                ri != leaf->rays.end(); ri++)
833                {
834                        if ((*ri).mRay->IsActive())
835                        {
836                                // first unref ray from the former leaf
837                                (*ri).mRay->Unref();
838
839                                // determine the side of this ray with respect to the plane
840                                int side = node->ComputeRayIntersection(*ri, (*ri).mRay->mT);
841
842                                if (side == 0)
843                                {
844                                        if ((*ri).mRay->HasPosDir(axis))
845                                        {
846                                                back->AddRay(VspKdTreeNode::RayInfo((*ri).mRay,
847                                                                                                                        (*ri).mMinT,
848                                                                                                                        (*ri).mRay->mT));
849                                                front->AddRay(VspKdTreeNode::RayInfo((*ri).mRay,
850                                                                                                                         (*ri).mRay->mT,
851                                                                                                                         (*ri).mMaxT));
852                                        }
853                                        else
854                                        {
855                                                back->AddRay(VspKdTreeNode::RayInfo((*ri).mRay,
856                                                                                                                        (*ri).mRay->mT,
857                                                                                                                        (*ri).mMaxT));
858                                                front->AddRay(VspKdTreeNode::RayInfo((*ri).mRay,
859                                                                                                                         (*ri).mMinT,
860                                                                                                                         (*ri).mRay->mT));
861                                        }
862                                }
863                                else
864                                        if (side == 1)
865                                                front->AddRay(*ri);
866                                        else
867                                                back->AddRay(*ri);
868                        } else
869                                (*ri).mRay->Unref();
870                }
871        }
872        else
873        {
874                // rays front/back
875   
876        for(VspKdTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
877                        ri != leaf->rays.end(); ++ ri)
878                {
879                        if ((*ri).mRay->IsActive())
880                        {
881                                // first unref ray from the former leaf
882                                (*ri).mRay->Unref();
883
884                                int side;
885                                if ((*ri).mRay->GetDirParametrization(axis - 3) > position)
886                                        side = 1;
887                                else
888                                        side = -1;
889                               
890                                if (side == 1)
891                                        front->AddRay(*ri);
892                                else
893                                        back->AddRay(*ri);             
894                        }
895                        else
896                                (*ri).mRay->Unref();
897                }
898        }
899       
900        // update stats
901        stat.rayRefs -= leaf->rays.size();
902        stat.rayRefs += raysBack + raysFront;
903
904        delete leaf;
905        return node;
906}
907
908int VspKdTree::ReleaseMemory(const int time)
909{
910        stack<VspKdTreeNode *> tstack;
911
912        // find a node in the tree which subtree will be collapsed
913        int maxAccessTime = time - accessTimeThreshold;
914        int released;
915
916        tstack.push(root);
917
918        while (!tstack.empty())
919        {
920                VspKdTreeNode *node = tstack.top();
921                tstack.pop();
922   
923                if (!node->IsLeaf())
924                {
925                        VspKdTreeInterior *in = (VspKdTreeInterior *)node;
926                        // cout<<"depth="<<(int)in->depth<<" time="<<in->lastAccessTime<<endl;
927                        if (in->depth >= minCollapseDepth && in->lastAccessTime <= maxAccessTime)
928                        {
929                                released = CollapseSubtree(node, time);
930                                break;
931                        }
932                        if (in->back->GetAccessTime() <  in->front->GetAccessTime())
933                        {
934                                tstack.push(in->front);
935                                tstack.push(in->back);
936                        }
937                        else
938                        {
939                                tstack.push(in->back);
940                                tstack.push(in->front);
941                        }
942                }
943        }
944
945        while (tstack.empty())
946        {
947                // could find node to collaps...
948                // cout<<"Could not find a node to release "<<endl;
949                break;
950        }
951 
952        return released;
953}
954
955
956
957
958VspKdTreeNode *VspKdTree::SubdivideLeaf(VspKdTreeLeaf *leaf,
959                                                                                const float sizeThreshold)
960{
961        VspKdTreeNode *node = leaf;
962
963        AxisAlignedBox3 leafBBox = GetBBox(leaf);
964
965        static int pass = 0;
966        ++ pass;
967       
968        // check if we should perform a dynamic subdivision of the leaf
969        if (// leaf->rays.size() > (unsigned)termMinCost &&
970                (leaf->GetPvsSize() >= termMinPvs) &&
971                (SqrMagnitude(leafBBox.Size()) > sizeThreshold) )
972        {
973        // memory check and realese...
974                if (GetMemUsage() > maxTotalMemory)
975                        ReleaseMemory(pass);
976               
977                AxisAlignedBox3 backBBox, frontBBox;
978
979                // subdivide the node
980                node = SubdivideNode(leaf, leafBBox, backBBox, frontBBox);
981        }
982        return node;
983}
984
985
986
987void
988VspKdTree::UpdateRays(VssRayContainer &remove,
989                                          VssRayContainer &add)
990{
991        VspKdTreeLeaf::NewMail();
992
993        // schedule rays for removal
994        for(VssRayContainer::const_iterator ri = remove.begin();
995                ri != remove.end(); ++ ri)
996        {
997                (*ri)->ScheduleForRemoval(); 
998        }
999
1000        int inactive = 0;
1001
1002        for (VssRayContainer::const_iterator ri = remove.begin(); ri != remove.end(); ++ ri)
1003        {
1004                if ((*ri)->ScheduledForRemoval())
1005                        //      RemoveRay(*ri, NULL, false);
1006                        // !!! BUG - with true it does not work correctly - aggreated delete
1007                        RemoveRay(*ri, NULL, true);
1008                else
1009                        ++ inactive;
1010        }
1011
1012        //  cout<<"all/inactive"<<remove.size()<<"/"<<inactive<<endl;
1013 
1014        for(VssRayContainer::const_iterator ri = add.begin(); ri != add.end(); ++ ri)
1015        {
1016                AddRay(*ri);
1017        }
1018}
1019
1020
1021void VspKdTree::RemoveRay(VssRay *ray,
1022                                                  vector<VspKdTreeLeaf *> *affectedLeaves,
1023                                                  const bool removeAllScheduledRays)
1024{
1025        stack<RayTraversalData> tstack;
1026
1027        tstack.push(RayTraversalData(root, VspKdTreeNode::RayInfo(ray)));
1028 
1029        RayTraversalData data;
1030
1031        // cout<<"Number of ray refs = "<<ray->RefCount()<<endl;
1032        while (!tstack.empty())
1033        {
1034                data = tstack.top();
1035                tstack.pop();
1036
1037                if (!data.node->IsLeaf())
1038                {
1039                        // split the set of rays in two groups intersecting the
1040                        // two subtrees
1041
1042                        TraverseInternalNode(data, tstack);
1043        }
1044                else
1045                {
1046                        // remove the ray from the leaf
1047                        // find the ray in the leaf and swap it with the last ray...
1048                        VspKdTreeLeaf *leaf = (VspKdTreeLeaf *)data.node;
1049     
1050                        if (!leaf->Mailed())
1051                        {
1052                                leaf->Mail();
1053                               
1054                                if (affectedLeaves)
1055                                        affectedLeaves->push_back(leaf);
1056       
1057                                if (removeAllScheduledRays)
1058                                {
1059                                        int tail = leaf->rays.size() - 1;
1060
1061                                        for (int i=0; i < (int)(leaf->rays.size()); ++ i)
1062                                        {
1063                                                if (leaf->rays[i].mRay->ScheduledForRemoval())
1064                                                {
1065                                                        // find a ray to replace it with
1066                                                        while (tail >= i && leaf->rays[tail].mRay->ScheduledForRemoval())
1067                                                        {
1068                                                                ++ stat.removedRayRefs;
1069                                                                leaf->rays[tail].mRay->Unref();
1070                                                                leaf->rays.pop_back();
1071                                                               
1072                                                                -- tail;
1073                                                        }
1074                                                       
1075                                                        if (tail < i)
1076                                                                break;
1077             
1078                                                        ++ stat.removedRayRefs;
1079             
1080                                                        leaf->rays[i].mRay->Unref();
1081                                                        leaf->rays[i] = leaf->rays[tail];
1082                                                        leaf->rays.pop_back();
1083                                                        -- tail;
1084                                                }
1085                                        }
1086                                }
1087                        }
1088     
1089                        if (!removeAllScheduledRays)
1090                                for (int i=0; i < (int)leaf->rays.size(); i++)
1091                                {
1092                                        if (leaf->rays[i].mRay == ray)
1093                                        {
1094                                                ++ stat.removedRayRefs;
1095                                                ray->Unref();
1096                                                leaf->rays[i] = leaf->rays[leaf->rays.size() - 1];
1097                                                leaf->rays.pop_back();
1098                                            // check this ray again
1099                                                break;
1100                                        }
1101                                }
1102                }
1103        }
1104
1105        if (ray->RefCount() != 0)
1106        {
1107                cerr << "Error: Number of remaining refs = " << ray->RefCount() << endl;
1108                exit(1);
1109        }
1110
1111}
1112
1113void VspKdTree::AddRay(VssRay *ray)
1114{
1115        stack<RayTraversalData> tstack;
1116 
1117        tstack.push(RayTraversalData(root, VspKdTreeNode::RayInfo(ray)));
1118 
1119        RayTraversalData data;
1120 
1121        while (!tstack.empty())
1122        {
1123                data = tstack.top();
1124                tstack.pop();
1125
1126                if (!data.node->IsLeaf())
1127                {
1128                        TraverseInternalNode(data, tstack);
1129                }
1130                else
1131                {
1132                        // remove the ray from the leaf
1133                        // find the ray in the leaf and swap it with the last ray...
1134                        VspKdTreeLeaf *leaf = (VspKdTreeLeaf *)data.node;
1135                        leaf->AddRay(data.rayData);
1136                        ++ stat.addedRayRefs;
1137                }
1138        }
1139}
1140
1141void VspKdTree::TraverseInternalNode(RayTraversalData &data,
1142                                                                         stack<RayTraversalData> &tstack)
1143{
1144        VspKdTreeInterior *in =  (VspKdTreeInterior *) data.node;
1145 
1146        if (in->axis <= VspKdTreeNode::SPLIT_Z)
1147        {
1148            // determine the side of this ray with respect to the plane
1149                int side = in->ComputeRayIntersection(data.rayData, data.rayData.mRay->mT);
1150   
1151      if (side == 0)
1152          {
1153                if (data.rayData.mRay->HasPosDir(in->axis))
1154                {
1155                        tstack.push(RayTraversalData(in->back,
1156                                                VspKdTreeNode::RayInfo(data.rayData.mRay, data.rayData.mMinT,
1157                                                                                           data.rayData.mRay->mT)));
1158                               
1159                        tstack.push(RayTraversalData(in->front,
1160                                                VspKdTreeNode::RayInfo(data.rayData.mRay, data.rayData.mRay->mT,
1161                                                                                           data.rayData.mMaxT)));
1162       
1163                }
1164                else
1165                {
1166                        tstack.push(RayTraversalData(in->back,
1167                                                                                 VspKdTreeNode::RayInfo(data.rayData.mRay,
1168                                                                                                                                data.rayData.mRay->mT,
1169                                                                                                                                data.rayData.mMaxT)));
1170                        tstack.push(RayTraversalData(in->front,
1171                                                                                 VspKdTreeNode::RayInfo(data.rayData.mRay,
1172                                                                                                                                data.rayData.mMinT,
1173                                                                                                                                data.rayData.mRay->mT)));
1174                }
1175          }
1176          else
1177          if (side == 1)
1178                          tstack.push(RayTraversalData(in->front, data.rayData));
1179                  else
1180                          tstack.push(RayTraversalData(in->back, data.rayData));
1181        }
1182        else
1183        {
1184                // directional split
1185                if (data.rayData.mRay->GetDirParametrization(in->axis - 3) > in->position)
1186                        tstack.push(RayTraversalData(in->front, data.rayData));
1187                else
1188                        tstack.push(RayTraversalData(in->back, data.rayData));
1189        }
1190}
1191
1192
1193int VspKdTree::CollapseSubtree(VspKdTreeNode *sroot, const int time)
1194{
1195        // first count all rays in the subtree
1196        // use mail 1 for this purpose
1197        stack<VspKdTreeNode *> tstack;
1198       
1199        int rayCount = 0;
1200        int totalRayCount = 0;
1201        int collapsedNodes = 0;
1202
1203#if DEBUG_COLLAPSE
1204        cout<<"Collapsing subtree"<<endl;
1205        cout<<"acessTime="<<sroot->GetAccessTime()<<endl;
1206        cout<<"depth="<<(int)sroot->depth<<endl;
1207#endif
1208
1209        // tstat.collapsedSubtrees++;
1210        // tstat.collapseDepths += (int)sroot->depth;
1211        // tstat.collapseAccessTimes += time - sroot->GetAccessTime();
1212        tstack.push(sroot);
1213        VssRay::NewMail();
1214       
1215        while (!tstack.empty())
1216        {
1217                ++ collapsedNodes;
1218               
1219                VspKdTreeNode *node = tstack.top();
1220                tstack.pop();
1221                if (node->IsLeaf())
1222                {
1223                        VspKdTreeLeaf *leaf = (VspKdTreeLeaf *) node;
1224                       
1225                        for(VspKdTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
1226                                        ri != leaf->rays.end(); ++ ri)
1227                        {
1228                                ++ totalRayCount;
1229                               
1230                                if ((*ri).mRay->IsActive() && !(*ri).mRay->Mailed())
1231                                {
1232                                        (*ri).mRay->Mail();
1233                                        ++ rayCount;
1234                                }
1235                        }
1236                }
1237                else
1238                {
1239                        tstack.push(((VspKdTreeInterior *)node)->back);
1240                        tstack.push(((VspKdTreeInterior *)node)->front);
1241                }
1242        }
1243 
1244        VssRay::NewMail();
1245
1246        // create a new node that will hold the rays
1247        VspKdTreeLeaf *newLeaf = new VspKdTreeLeaf(sroot->parent, rayCount);
1248       
1249        if (newLeaf->parent)
1250                newLeaf->parent->ReplaceChildLink(sroot, newLeaf);
1251
1252        tstack.push(sroot);
1253
1254        while (!tstack.empty())
1255        {
1256                VspKdTreeNode *node = tstack.top();
1257                tstack.pop();
1258
1259                if (node->IsLeaf())
1260                {
1261                        VspKdTreeLeaf *leaf = (VspKdTreeLeaf *) node;
1262
1263                        for(VspKdTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
1264                                        ri != leaf->rays.end(); ++ ri)
1265                        {
1266                                // unref this ray from the old node             
1267                                if ((*ri).mRay->IsActive())
1268                                {
1269                                        (*ri).mRay->Unref();
1270                                        if (!(*ri).mRay->Mailed())
1271                                        {
1272                                                (*ri).mRay->Mail();
1273                                                newLeaf->AddRay(*ri);
1274                                        }
1275                                }
1276                                else
1277                                        (*ri).mRay->Unref();
1278                        }
1279                }
1280                else
1281                {
1282                        tstack.push(((VspKdTreeInterior *)node)->back);
1283                        tstack.push(((VspKdTreeInterior *)node)->front);
1284                }
1285        }
1286 
1287        // delete the node and all its children
1288        delete sroot;
1289 
1290        // for(VspKdTreeNode::SRayContainer::iterator ri = newLeaf->rays.begin();
1291    //      ri != newLeaf->rays.end(); ++ ri)
1292        //     (*ri).ray->UnMail(2);
1293
1294#if DEBUG_COLLAPSE
1295        cout<<"Total memory before="<<GetMemUsage()<<endl;
1296#endif
1297
1298        stat.nodes -= collapsedNodes - 1;
1299        stat.rayRefs -= totalRayCount - rayCount;
1300 
1301#if DEBUG_COLLAPSE
1302        cout << "collapsed nodes" << collapsedNodes << endl;
1303        cout << "collapsed rays" << totalRayCount - rayCount << endl;
1304        cout << "Total memory after=" << GetMemUsage() << endl;
1305        cout << "================================" << endl;
1306#endif
1307
1308        //  tstat.collapsedNodes += collapsedNodes;
1309        //  tstat.collapsedRays += totalRayCount - rayCount;
1310    return totalRayCount - rayCount;
1311}
1312
1313
1314int VspKdTree::GetPvsSize(VspKdTreeNode *node, const AxisAlignedBox3 &box) const
1315{
1316        stack<VspKdTreeNode *> tstack;
1317        tstack.push(root);
1318
1319        Intersectable::NewMail();
1320        int pvsSize = 0;
1321       
1322        while (!tstack.empty())
1323        {
1324                VspKdTreeNode *node = tstack.top();
1325                tstack.pop();
1326   
1327          if (node->IsLeaf())
1328                {
1329                        VspKdTreeLeaf *leaf = (VspKdTreeLeaf *)node;
1330                        for (VspKdTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin();
1331                                ri != leaf->rays.end(); ++ ri)
1332                        {
1333                                if ((*ri).mRay->IsActive())
1334                                {
1335                                        Intersectable *object;
1336#if BIDIRECTIONAL_RAY
1337                                        object = (*ri).mRay->mOriginObject;
1338                                       
1339                                        if (object && !object->Mailed())
1340                                        {
1341                                                ++ pvsSize;
1342                                                object->Mail();
1343                                        }
1344#endif
1345                                        object = (*ri).mRay->mTerminationObject;
1346                                        if (object && !object->Mailed())
1347                                        {
1348                                                ++ pvsSize;
1349                                                object->Mail();
1350                                        }
1351                                }
1352                        }
1353                }
1354                else
1355                {
1356                        VspKdTreeInterior *in = (VspKdTreeInterior *)node;
1357       
1358                        if (in->axis < 3)
1359                        {
1360                                if (box.Max(in->axis) >= in->position)
1361                                        tstack.push(in->front);
1362                               
1363                                if (box.Min(in->axis) <= in->position)
1364                                        tstack.push(in->back);
1365                        }
1366                        else
1367                        {
1368                                // both nodes for directional splits
1369                                tstack.push(in->front);
1370                                tstack.push(in->back);
1371                        }
1372                }
1373        }
1374
1375        return pvsSize;
1376}
1377
1378void VspKdTree::GetRayContributionStatistics(float &minRayContribution,
1379                                                                                         float &maxRayContribution,
1380                                                                                         float &avgRayContribution)
1381{
1382        stack<VspKdTreeNode *> tstack;
1383        tstack.push(root);
1384       
1385        minRayContribution = 1.0f;
1386        maxRayContribution = 0.0f;
1387
1388        float sumRayContribution = 0.0f;
1389        int leaves = 0;
1390
1391        while (!tstack.empty())
1392        {
1393                VspKdTreeNode *node = tstack.top();
1394                tstack.pop();
1395                if (node->IsLeaf())
1396                {
1397                        leaves++;
1398                        VspKdTreeLeaf *leaf = (VspKdTreeLeaf *)node;
1399                        float c = leaf->GetAvgRayContribution();
1400
1401                        if (c > maxRayContribution)
1402                                maxRayContribution = c;
1403                        if (c < minRayContribution)
1404                                minRayContribution = c;
1405                        sumRayContribution += c;
1406                       
1407                }
1408                else {
1409                        VspKdTreeInterior *in = (VspKdTreeInterior *)node;
1410                        // both nodes for directional splits
1411                        tstack.push(in->front);
1412                        tstack.push(in->back);
1413                }
1414        }
1415       
1416        cout<<"sum="<<sumRayContribution<<endl;
1417        cout<<"leaves="<<leaves<<endl;
1418        avgRayContribution = sumRayContribution/(float)leaves;
1419}
1420
1421
1422int VspKdTree::GenerateRays(const float ratioPerLeaf,
1423                                                        SimpleRayContainer &rays)
1424{
1425        stack<VspKdTreeNode *> tstack;
1426        tstack.push(root);
1427       
1428        while (!tstack.empty())
1429        {
1430                VspKdTreeNode *node = tstack.top();
1431                tstack.pop();
1432               
1433                if (node->IsLeaf())
1434                {
1435                        VspKdTreeLeaf *leaf = (VspKdTreeLeaf *)node;
1436                        float c = leaf->GetAvgRayContribution();
1437                        int num = (c*ratioPerLeaf + 0.5);
1438                        //                      cout<<num<<" ";
1439
1440                        for (int i=0; i < num; i++)
1441                        {
1442                                Vector3 origin = GetBBox(leaf).GetRandomPoint();
1443                                Vector3 dirVector = GetDirBBox(leaf).GetRandomPoint();
1444                                Vector3 direction = Vector3(sin(dirVector.x), sin(dirVector.y), cos(dirVector.x));
1445                                //cout<<"dir vector.x="<<dirVector.x<<"direction'.x="<<atan2(direction.x, direction.y)<<endl;
1446                                rays.push_back(SimpleRay(origin, direction));
1447                        }
1448                       
1449                }
1450                else
1451                {
1452                        VspKdTreeInterior *in = (VspKdTreeInterior *)node;
1453                        // both nodes for directional splits
1454                        tstack.push(in->front);
1455                        tstack.push(in->back);
1456                }
1457        }
1458
1459        return rays.size();
1460}
1461
1462
1463float VspKdTree::GetAvgPvsSize()
1464{
1465        stack<VspKdTreeNode *> tstack;
1466        tstack.push(root);
1467
1468        int sumPvs = 0;
1469        int leaves = 0;
1470       
1471        while (!tstack.empty())
1472        {
1473                VspKdTreeNode *node = tstack.top();
1474                tstack.pop();
1475               
1476                if (node->IsLeaf())
1477                {
1478                        VspKdTreeLeaf *leaf = (VspKdTreeLeaf *)node;
1479                        // update pvs size
1480                        leaf->UpdatePvsSize();
1481                        sumPvs += leaf->GetPvsSize();
1482                        leaves++;
1483                }
1484                else
1485                {
1486                        VspKdTreeInterior *in = (VspKdTreeInterior *)node;
1487                        // both nodes for directional splits
1488                        tstack.push(in->front);
1489                        tstack.push(in->back);
1490                }
1491        }
1492
1493        return sumPvs / (float)leaves;
1494}
Note: See TracBrowser for help on using the repository browser.