source: GTP/trunk/App/Demos/Vis/FriendlyCulling/src/BvhConstructor.cpp @ 3295

Revision 3295, 10.0 KB checked in by mattausch, 16 years ago (diff)
Line 
1#include "BvhConstructor.h"
2#include "Triangle3.h"
3#include "SceneEntity.h"
4#include "Geometry.h"
5#include "gzstream.h"
6
7#include <queue>
8#include <stack>
9#include <fstream>
10
11#define MAX_FLOAT 1e20f
12
13
14#ifdef _CRT_SET
15        #define _CRTDBG_MAP_ALLOC
16        #include <stdlib.h>
17        #include <crtdbg.h>
18
19        // redefine new operator
20        #define DEBUG_NEW new(_NORMAL_BLOCK, __FILE__, __LINE__)
21        #define new DEBUG_NEW
22#endif
23
24
25using namespace std;
26
27
28namespace CHCDemoEngine
29{
30
31BvhConstructor::BvhConstructor(SceneEntity **entities,
32                                                           int first,
33                                                           int last):
34mEntities(entities),
35mFirst(first),
36mLast(last),
37mMaxDepth(20),
38mMaxObjects(1),
39mMaxTriangles(1),
40mNumNodes(0),
41//mSplitType(SAH)
42mSplitType(SAH_OR_SIZE)
43{
44}
45
46
47float BvhConstructor::EvaluateSahCost(BvhLeaf *leaf, int axis, float position)
48{
49        // count triangles on the left / right
50        int left = 0;
51        int right = 0;
52        AxisAlignedBox3 leftBox, rightBox;
53
54        leftBox.Initialize();
55        rightBox.Initialize();
56
57        const int candidates = std::max(50, (int)(leaf->CountPrimitives() / 20));
58        const float finc = leaf->CountPrimitives() / (float)candidates;
59
60        int i = leaf->mFirst;
61        float fi = leaf->mFirst + 0.5f;
62
63        AxisAlignedBox3 box;
64
65        for (; i <= leaf->mLast;)
66        {
67                box = mEntities[i]->GetWorldBoundingBox();
68
69                if (box.Center(axis) < position)
70                {
71                        ++ left;
72                        // update the box
73                        leftBox.Include(box);
74                }
75                else
76                {
77                        ++ right;
78                        rightBox.Include(box);
79                }
80                 
81                fi += finc;
82                i = (int)fi;
83        }
84
85        float bW = 1.0f;
86        float leftRatio = left / (float)leaf->CountPrimitives();
87        float rightRatio = right / (float)leaf->CountPrimitives();
88
89        float saLeft = 0.0f;
90        float saRight = 0.0f;
91
92        // not a valid split
93        if (!left || !right)
94                return 1e25f;
95
96        saLeft = leftBox.SurfaceArea();
97        saRight = rightBox.SurfaceArea();
98
99
100        return
101                saLeft  * ((1.0f - bW) + bW * leftRatio) +
102                saRight * ((1.0f - bW) + bW * rightRatio);
103}
104
105
106float BvhConstructor::SelectPlaneSah(BvhLeaf *leaf, int &axis, float &minCost)
107{
108        minCost = MAX_FLOAT;
109        float bestPos = minCost;
110        int bestAxis = 0;
111
112        const int initialPlanes = 3;
113
114        // initiate the costs
115        for (axis = 0; axis < 3; ++ axis)
116        {
117                const float size = leaf->mBox.Size(axis);
118
119                for (int i = 0; i < initialPlanes; ++ i)
120                {
121                        const float pos = leaf->mBox.Min(axis) + (i + 1) * size / (initialPlanes + 1);
122                        const float cost = EvaluateSahCost(leaf, axis, pos);
123
124                        if (cost < minCost)
125                        {
126                                minCost = cost;
127                                bestPos = pos;
128                                bestAxis = axis;
129                        }
130                }
131        }
132
133        axis = bestAxis;
134
135        //  cout<<axis<<endl<<flush;
136        const float shrink = 0.5f;
137
138        // now hierarchically refine the initial interval
139        float size = shrink *
140                (leaf->mBox.Max(axis) - leaf->mBox.Min(axis)) / (initialPlanes + 1);
141
142        const int steps = 4;
143        float cost;
144
145        for (int i = 0; i < steps; ++ i)
146        {
147                const float left = bestPos - size;
148                const float right = bestPos + size;
149
150                cost = EvaluateSahCost(leaf, axis, left);
151
152                if (cost < minCost)
153                {
154                        minCost = cost;
155                        bestPos = left;
156                }
157
158                cost = EvaluateSahCost(leaf, axis, right);
159
160                if (cost < minCost)
161                {
162                        minCost = cost;
163                        bestPos = right;
164                }
165
166                size = shrink * size;
167        }
168
169        //OUT1("best axis: " << axis << " " << bestPos << " " << minCost);
170
171        return bestPos;
172}
173
174
175int BvhConstructor::SortTriangles(BvhLeaf *leaf,
176                                                                  int axis,
177                                                                  float position
178                                                                  )
179{
180        int i = leaf->mFirst;
181        int j = leaf->mLast;
182
183    while (1)
184        {
185                //while (mEntities[i]->GetWorldCenter()[axis] < position) ++ i;
186                //while (position < mEntities[j]->GetWorldCenter()[axis]) -- j;
187
188                while ((i < leaf->mLast) && (mEntities[i]->GetWorldCenter()[axis] < position)) ++ i;
189                while ((j > leaf->mFirst) && (position < mEntities[j]->GetWorldCenter()[axis])) -- j;
190
191                // sorting finished
192                if (i >= j) break;
193
194                // swap entities
195                swap(mEntities[i], mEntities[j]);
196                       
197                ++ i;
198                -- j;
199        }
200
201        return j;
202}
203
204
205void BvhConstructor::UpdateBoundingBoxes(BvhNode *node)
206{
207        if (node->IsLeaf())
208        {
209                int numEntities = node->CountPrimitives();
210                SceneEntity **entities = mEntities + node->mFirst;
211
212                node->mBox = SceneEntity::ComputeBoundingBox(entities, numEntities);
213                //cout << "box: " << node->mBox << endl;
214        }
215        else
216        {
217                BvhNode *f = static_cast<BvhInterior *>(node)->GetFront();
218                BvhNode *b = static_cast<BvhInterior *>(node)->GetBack();
219
220                UpdateBoundingBoxes(f);
221                UpdateBoundingBoxes(b);
222
223                node->mBox = f->mBox;
224                node->mBox.Include(b->mBox);
225        }
226}
227
228
229int BvhConstructor::SortTrianglesSpatialMedian(BvhLeaf *leaf,
230                                                                                           int axis)
231{
232        // spatial median
233        float m = leaf->GetBox().Center()[axis];
234        return SortTriangles(leaf, axis, m);
235}
236
237
238int BvhConstructor::SortTrianglesObjectMedian(BvhLeaf *leaf, int axis, float &pos)
239{
240        // Now distribute the objects into the subnodes
241        // Use QuickMedian sort
242        int l = leaf->mFirst;
243        int r = leaf->mLast;
244        int k = (l + r) / 2;
245
246        while (l < r)
247        {
248                int i = l;
249                int j = r;
250
251                // get some estimation of the pivot
252                pos = mEntities[(l + r) / 2]->GetBoundingBox().Center(axis);
253
254                while (1)
255                {
256                        while ((i <= leaf->mLast) && (mEntities[i]->GetWorldBoundingBox().Center(axis) < pos))
257                                ++ i;
258
259                        while((j >= leaf->mFirst) && (pos < mEntities[j]->GetWorldBoundingBox().Center(axis)))
260                                -- j;
261
262                        if (i <= j)
263                        {
264                                std::swap(mEntities[i], mEntities[j]);
265                                ++ i;
266                                -- j;
267                        }
268                        else
269                        {
270                                break;
271                        }
272                }
273
274                // now check the extents
275                if (i >= k)
276                        r = j;
277                else
278                        l = i;
279        }
280
281        return k;
282}
283
284
285int BvhConstructor::SortTrianglesSurfaceArea(BvhLeaf *leaf, float sa)
286{
287        int i = leaf->mFirst;
288        int j = leaf->mLast;
289
290        while(1)
291        {
292                while ((i <= j) && (mEntities[i]->GetWorldBoundingBox().SurfaceArea() < sa))
293                        ++ i;
294
295                while ((i <= j) && (sa < mEntities[j]->GetWorldBoundingBox().SurfaceArea()))
296                        -- j;
297
298                if (i < j)
299                {
300                        swap(mEntities[i], mEntities[j]);
301                        ++ i;
302                        -- j;
303                }
304                else
305                        break;
306        }
307
308        return j;
309}
310
311
312
313BvhNode *BvhConstructor::SubdivideLeaf(BvhLeaf *leaf,
314                                                                           int parentAxis
315                                                                           )
316{
317        if (TerminationCriteriaMet(leaf))
318        {
319                //if (leaf->CountPrimitives() <= 0)
320                //cout << "leaf constructed:" << leaf->mBox << " " << leaf->mFirst << " " << leaf->mLast << endl;
321                return leaf;
322        }
323
324        //const int axis = (parentAxis + 1) % 3;
325        int axis = leaf->mBox.MajorAxis();
326        const int scale = 20;
327
328        // position of the split in the partailly sorted array of triangles
329        // corresponding to this leaf
330        int split = -1;
331        float pos = -1.0f;
332
333        // Spatial median subdivision
334        switch (mSplitType)
335        {
336        case SPATIAL_MEDIAN:
337                split = SortTrianglesSpatialMedian(leaf, axis);
338                pos = leaf->mBox.Center()[axis];
339                break;
340        case OBJECT_MEDIAN:
341                // Object median subdivision: approximately the same number
342                // of objects on the left of the the splitting point.
343                split = SortTrianglesObjectMedian(leaf, axis, pos);
344                break;
345        case SAH:
346                {
347                        float cost;
348                        pos = SelectPlaneSah(leaf, axis, cost);
349
350                        if (pos != MAX_FLOAT)
351                        {
352                                split = SortTriangles(leaf, axis, pos);
353                        }
354
355                        if ((pos == MAX_FLOAT) || (split == leaf->mLast))
356                        {
357                                split = -1;
358                                split = SortTrianglesObjectMedian(leaf, axis, pos);
359                        }
360                }
361                break;
362        case SAH_OR_SIZE:
363                {
364                        // split by size instead
365                        const float saThreshold = 0.2f * leaf->GetBox().SurfaceArea();
366                        split = SortTrianglesSurfaceArea(leaf, saThreshold);
367
368                        if ((split == leaf->mLast) || (split == leaf->mFirst - 1))
369                        {
370                                // use SAH
371                                float cost;
372                                pos = SelectPlaneSah(leaf, axis, cost);
373
374                                if (pos != MAX_FLOAT)
375                                        split = SortTriangles(leaf, axis, pos);
376                                else
377                                        split = SortTrianglesObjectMedian(leaf, axis, pos);
378                        }
379                        else
380                        {
381                                // note: no position is computed!!
382                                //OUT1("sorted by size");
383                        }
384                }
385                break;
386        default:
387                cerr << "should not come here" << endl;
388                break;
389        }
390       
391        if (split == leaf->mLast)
392        {
393                // no split could be achieved => just halve number of objects
394                split = (leaf->mLast + leaf->mFirst) / 2;
395                //cerr << "no reduction " << leaf->CountPrimitives() << " " << leaf->mFirst << " " << leaf->mLast << endl;
396        }
397
398        // create two more nodes
399        mNumNodes += 2;
400        BvhInterior *parent = new BvhInterior(leaf->GetParent());
401        parent->mFirst = leaf->mFirst;
402        parent->mLast = leaf->mLast;
403       
404        parent->mAxis = axis;
405        parent->mBox = leaf->mBox;
406        parent->mDepth = leaf->mDepth;
407
408        BvhLeaf *front = new BvhLeaf(parent);
409
410        parent->mBack = leaf;
411        parent->mFront = front;
412               
413        // now assign the triangles to the subnodes
414        front->mFirst = split + 1;
415        front->mLast = leaf->mLast;
416        front->mDepth = leaf->mDepth + 1;
417
418        leaf->mLast = split;
419        leaf->mDepth = front->mDepth;
420        leaf->mParent = parent;
421       
422        front->mBox = SceneEntity::ComputeBoundingBox(mEntities + front->mFirst, front->CountPrimitives());
423        leaf->mBox = SceneEntity::ComputeBoundingBox(mEntities + leaf->mFirst, leaf->CountPrimitives());
424
425        // recursively continue subdivision
426        parent->mBack = SubdivideLeaf(static_cast<BvhLeaf *>(parent->mBack), axis);
427        parent->mFront = SubdivideLeaf(static_cast<BvhLeaf *>(parent->mFront), axis);
428       
429        return parent;
430}
431
432
433int BvhConstructor::CountTriangles(BvhNode *node) const
434{
435        int numTriangles = 0;
436
437        for (int i = node->mFirst; i <= node->mLast; ++ i)
438        {
439                numTriangles += mEntities[i]->CountNumTriangles(0);
440        }
441
442        return numTriangles;
443}
444
445
446bool BvhConstructor::TerminationCriteriaMet(BvhLeaf *leaf) const
447{
448        const bool criteriaMet =
449                (leaf->mDepth > mMaxDepth) ||
450                (leaf->CountPrimitives() <= mMaxObjects) ||
451                (CountTriangles(leaf) <= mMaxTriangles);
452
453        return criteriaMet;
454}
455
456
457BvhNode *BvhConstructor::Construct(int &numNodes)
458{
459        BvhNode *root;
460        mNumNodes = 1;
461
462        BvhLeaf *l = new BvhLeaf(NULL);
463       
464        l->mFirst = mFirst;
465        l->mLast = mLast;
466
467        cout << "\n================================" << endl;
468        cout << "constructing bvh for " << l->mLast - l->mFirst + 1 << " entities" << endl;
469
470        l->mBox = SceneEntity::ComputeBoundingBox(mEntities + mFirst, mLast - mFirst + 1);
471        l->mArea = l->mBox.SurfaceArea();
472       
473        root = SubdivideLeaf(l, 0);
474
475        UpdateBoundingBoxes(root);
476
477        numNodes = mNumNodes;
478
479        return root;
480}
481
482
483}
Note: See TracBrowser for help on using the repository browser.