[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 | 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 |
|
---|
[3262] | 175 | int BvhConstructor::SortTriangles(BvhLeaf *leaf,
|
---|
[3267] | 176 | int axis,
|
---|
| 177 | float position
|
---|
| 178 | )
|
---|
[3262] | 179 | {
|
---|
| 180 | int i = leaf->mFirst;
|
---|
| 181 | int j = leaf->mLast;
|
---|
| 182 |
|
---|
| 183 | while (1)
|
---|
| 184 | {
|
---|
[3265] | 185 | //while (mEntities[i]->GetWorldCenter()[axis] < position) ++ i;
|
---|
| 186 | //while (position < mEntities[j]->GetWorldCenter()[axis]) -- j;
|
---|
[3262] | 187 |
|
---|
[3265] | 188 | while ((i < leaf->mLast) && (mEntities[i]->GetWorldCenter()[axis] < position)) ++ i;
|
---|
| 189 | while ((j > leaf->mFirst) && (position < mEntities[j]->GetWorldCenter()[axis])) -- j;
|
---|
| 190 |
|
---|
[3262] | 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 |
|
---|
[3267] | 205 | void 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 |
|
---|
[3262] | 229 | int 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 |
|
---|
[3269] | 238 | int 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 |
|
---|
[3280] | 285 | int 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 |
|
---|
[3262] | 313 | BvhNode *BvhConstructor::SubdivideLeaf(BvhLeaf *leaf,
|
---|
| 314 | int parentAxis
|
---|
| 315 | )
|
---|
| 316 | {
|
---|
| 317 | if (TerminationCriteriaMet(leaf))
|
---|
| 318 | {
|
---|
[3267] | 319 | //if (leaf->CountPrimitives() <= 0)
|
---|
| 320 | //cout << "leaf constructed:" << leaf->mBox << " " << leaf->mFirst << " " << leaf->mLast << endl;
|
---|
[3262] | 321 | return leaf;
|
---|
| 322 | }
|
---|
| 323 |
|
---|
| 324 | //const int axis = (parentAxis + 1) % 3;
|
---|
[3269] | 325 | int axis = leaf->mBox.MajorAxis();
|
---|
[3262] | 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;
|
---|
[3269] | 332 |
|
---|
[3262] | 333 | // Spatial median subdivision
|
---|
[3269] | 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;
|
---|
[3280] | 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;
|
---|
[3269] | 386 | default:
|
---|
| 387 | cerr << "should not come here" << endl;
|
---|
| 388 | break;
|
---|
| 389 | }
|
---|
[3262] | 390 |
|
---|
| 391 | if (split == leaf->mLast)
|
---|
| 392 | {
|
---|
| 393 | // no split could be achieved => just halve number of objects
|
---|
[3265] | 394 | split = (leaf->mLast + leaf->mFirst) / 2;
|
---|
[3268] | 395 | //cerr << "no reduction " << leaf->CountPrimitives() << " " << leaf->mFirst << " " << leaf->mLast << endl;
|
---|
[3262] | 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 |
|
---|
[3269] | 433 | int 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 |
|
---|
[3262] | 446 | bool BvhConstructor::TerminationCriteriaMet(BvhLeaf *leaf) const
|
---|
| 447 | {
|
---|
| 448 | const bool criteriaMet =
|
---|
| 449 | (leaf->mDepth > mMaxDepth) ||
|
---|
[3269] | 450 | (leaf->CountPrimitives() <= mMaxObjects) ||
|
---|
| 451 | (CountTriangles(leaf) <= mMaxTriangles);
|
---|
[3262] | 452 |
|
---|
| 453 | return criteriaMet;
|
---|
| 454 | }
|
---|
| 455 |
|
---|
| 456 |
|
---|
| 457 | BvhNode *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;
|
---|
[3268] | 466 |
|
---|
[3269] | 467 | cout << "\n================================" << endl;
|
---|
| 468 | cout << "constructing bvh for " << l->mLast - l->mFirst + 1 << " entities" << endl;
|
---|
[3268] | 469 |
|
---|
| 470 | l->mBox = SceneEntity::ComputeBoundingBox(mEntities + mFirst, mLast - mFirst + 1);
|
---|
[3262] | 471 | l->mArea = l->mBox.SurfaceArea();
|
---|
| 472 |
|
---|
| 473 | root = SubdivideLeaf(l, 0);
|
---|
| 474 |
|
---|
[3267] | 475 | UpdateBoundingBoxes(root);
|
---|
| 476 |
|
---|
[3262] | 477 | numNodes = mNumNodes;
|
---|
| 478 |
|
---|
| 479 | return root;
|
---|
| 480 | }
|
---|
| 481 |
|
---|
| 482 |
|
---|
| 483 | } |
---|