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