00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 #include "stdafx.h"
00027 #include "accelerators/bvh.h"
00028 #include "probes.h"
00029 #include "paramset.h"
00030
00031
00032 struct BVHPrimitiveInfo {
00033 BVHPrimitiveInfo() { }
00034 BVHPrimitiveInfo(int pn, const BBox &b)
00035 : primitiveNumber(pn), bounds(b) {
00036 centroid = .5f * b.pMin + .5f * b.pMax;
00037 }
00038 int primitiveNumber;
00039 Point centroid;
00040 BBox bounds;
00041 };
00042
00043
00044 struct BVHBuildNode {
00045
00046 BVHBuildNode() { children[0] = children[1] = NULL; }
00047 void InitLeaf(uint32_t first, uint32_t n, const BBox &b) {
00048 firstPrimOffset = first;
00049 nPrimitives = n;
00050 bounds = b;
00051 }
00052 void InitInterior(uint32_t axis, BVHBuildNode *c0, BVHBuildNode *c1) {
00053 children[0] = c0;
00054 children[1] = c1;
00055 bounds = Union(c0->bounds, c1->bounds);
00056 splitAxis = axis;
00057 nPrimitives = 0;
00058 }
00059 BBox bounds;
00060 BVHBuildNode *children[2];
00061 uint32_t splitAxis, firstPrimOffset, nPrimitives;
00062 };
00063
00064
00065 struct CompareToMid {
00066 CompareToMid(int d, float m) { dim = d; mid = m; }
00067 int dim;
00068 float mid;
00069 bool operator()(const BVHPrimitiveInfo &a) const {
00070 return a.centroid[dim] < mid;
00071 }
00072 };
00073
00074
00075 struct ComparePoints {
00076 ComparePoints(int d) { dim = d; }
00077 int dim;
00078 bool operator()(const BVHPrimitiveInfo &a,
00079 const BVHPrimitiveInfo &b) const {
00080 return a.centroid[dim] < b.centroid[dim];
00081 }
00082 };
00083
00084
00085 struct CompareToBucket {
00086 CompareToBucket(int split, int num, int d, const BBox &b)
00087 : centroidBounds(b)
00088 { splitBucket = split; nBuckets = num; dim = d; }
00089 bool operator()(const BVHPrimitiveInfo &p) const;
00090
00091 int splitBucket, nBuckets, dim;
00092 const BBox ¢roidBounds;
00093 };
00094
00095
00096 bool CompareToBucket::operator()(const BVHPrimitiveInfo &p) const {
00097 int b = nBuckets * ((p.centroid[dim] - centroidBounds.pMin[dim]) /
00098 (centroidBounds.pMax[dim] - centroidBounds.pMin[dim]));
00099 if (b == nBuckets) b = nBuckets-1;
00100 Assert(b >= 0 && b < nBuckets);
00101 return b <= splitBucket;
00102 }
00103
00104
00105 struct LinearBVHNode {
00106 BBox bounds;
00107 union {
00108 uint32_t primitivesOffset;
00109 uint32_t secondChildOffset;
00110 };
00111
00112 uint8_t nPrimitives;
00113 uint8_t axis;
00114 uint8_t pad[2];
00115 };
00116
00117
00118 static inline bool IntersectP(const BBox &bounds, const Ray &ray,
00119 const Vector &invDir, const uint32_t dirIsNeg[3]) {
00120
00121 float tmin = (bounds[ dirIsNeg[0]].x - ray.o.x) * invDir.x;
00122 float tmax = (bounds[1-dirIsNeg[0]].x - ray.o.x) * invDir.x;
00123 float tymin = (bounds[ dirIsNeg[1]].y - ray.o.y) * invDir.y;
00124 float tymax = (bounds[1-dirIsNeg[1]].y - ray.o.y) * invDir.y;
00125 if ((tmin > tymax) || (tymin > tmax))
00126 return false;
00127 if (tymin > tmin) tmin = tymin;
00128 if (tymax < tmax) tmax = tymax;
00129
00130
00131 float tzmin = (bounds[ dirIsNeg[2]].z - ray.o.z) * invDir.z;
00132 float tzmax = (bounds[1-dirIsNeg[2]].z - ray.o.z) * invDir.z;
00133 if ((tmin > tzmax) || (tzmin > tmax))
00134 return false;
00135 if (tzmin > tmin)
00136 tmin = tzmin;
00137 if (tzmax < tmax)
00138 tmax = tzmax;
00139 return (tmin < ray.maxt) && (tmax > ray.mint);
00140 }
00141
00142
00143
00144
00145 BVHAccel::BVHAccel(const vector<Reference<Primitive> > &p,
00146 uint32_t mp, const string &sm) {
00147 maxPrimsInNode = min(255u, mp);
00148 for (uint32_t i = 0; i < p.size(); ++i)
00149 p[i]->FullyRefine(primitives);
00150 if (sm == "sah") splitMethod = SPLIT_SAH;
00151 else if (sm == "middle") splitMethod = SPLIT_MIDDLE;
00152 else if (sm == "equal") splitMethod = SPLIT_EQUAL_COUNTS;
00153 else {
00154 Warning("BVH split method \"%s\" unknown. Using \"sah\".",
00155 sm.c_str());
00156 splitMethod = SPLIT_SAH;
00157 }
00158
00159 if (primitives.size() == 0) {
00160 nodes = NULL;
00161 return;
00162 }
00163
00164 PBRT_BVH_STARTED_CONSTRUCTION(this, primitives.size());
00165
00166
00167 vector<BVHPrimitiveInfo> buildData;
00168 buildData.reserve(primitives.size());
00169 for (uint32_t i = 0; i < primitives.size(); ++i) {
00170 BBox bbox = primitives[i]->WorldBound();
00171 buildData.push_back(BVHPrimitiveInfo(i, bbox));
00172 }
00173
00174
00175 MemoryArena buildArena;
00176 uint32_t totalNodes = 0;
00177 vector<Reference<Primitive> > orderedPrims;
00178 orderedPrims.reserve(primitives.size());
00179 BVHBuildNode *root = recursiveBuild(buildArena, buildData, 0,
00180 primitives.size(), &totalNodes,
00181 orderedPrims);
00182 primitives.swap(orderedPrims);
00183 Info("BVH created with %d nodes for %d primitives (%.2f MB)", totalNodes,
00184 (int)primitives.size(), float(totalNodes * sizeof(LinearBVHNode))/(1024.f*1024.f));
00185
00186
00187 nodes = AllocAligned<LinearBVHNode>(totalNodes);
00188 for (uint32_t i = 0; i < totalNodes; ++i)
00189 new (&nodes[i]) LinearBVHNode;
00190 uint32_t offset = 0;
00191 flattenBVHTree(root, &offset);
00192 Assert(offset == totalNodes);
00193 PBRT_BVH_FINISHED_CONSTRUCTION(this);
00194 }
00195
00196
00197 BBox BVHAccel::WorldBound() const {
00198 return nodes ? nodes[0].bounds : BBox();
00199 }
00200
00201
00202 BVHBuildNode *BVHAccel::recursiveBuild(MemoryArena &buildArena,
00203 vector<BVHPrimitiveInfo> &buildData, uint32_t start,
00204 uint32_t end, uint32_t *totalNodes,
00205 vector<Reference<Primitive> > &orderedPrims) {
00206 Assert(start != end);
00207 (*totalNodes)++;
00208 BVHBuildNode *node = buildArena.Alloc<BVHBuildNode>();
00209
00210 BBox bbox;
00211 for (uint32_t i = start; i < end; ++i)
00212 bbox = Union(bbox, buildData[i].bounds);
00213 uint32_t nPrimitives = end - start;
00214 if (nPrimitives == 1) {
00215
00216 uint32_t firstPrimOffset = orderedPrims.size();
00217 for (uint32_t i = start; i < end; ++i) {
00218 uint32_t primNum = buildData[i].primitiveNumber;
00219 orderedPrims.push_back(primitives[primNum]);
00220 }
00221 node->InitLeaf(firstPrimOffset, nPrimitives, bbox);
00222 }
00223 else {
00224
00225 BBox centroidBounds;
00226 for (uint32_t i = start; i < end; ++i)
00227 centroidBounds = Union(centroidBounds, buildData[i].centroid);
00228 int dim = centroidBounds.MaximumExtent();
00229
00230
00231 uint32_t mid = (start + end) / 2;
00232 if (centroidBounds.pMax[dim] == centroidBounds.pMin[dim]) {
00233
00234 uint32_t firstPrimOffset = orderedPrims.size();
00235 for (uint32_t i = start; i < end; ++i) {
00236 uint32_t primNum = buildData[i].primitiveNumber;
00237 orderedPrims.push_back(primitives[primNum]);
00238 }
00239 node->InitLeaf(firstPrimOffset, nPrimitives, bbox);
00240 return node;
00241 }
00242
00243
00244 switch (splitMethod) {
00245 case SPLIT_MIDDLE: {
00246
00247 float pmid = .5f * (centroidBounds.pMin[dim] + centroidBounds.pMax[dim]);
00248 BVHPrimitiveInfo *midPtr = std::partition(&buildData[start],
00249 &buildData[end-1]+1,
00250 CompareToMid(dim, pmid));
00251 mid = midPtr - &buildData[0];
00252 break;
00253 }
00254 case SPLIT_EQUAL_COUNTS: {
00255
00256 mid = (start + end) / 2;
00257 std::nth_element(&buildData[start], &buildData[mid],
00258 &buildData[end-1]+1, ComparePoints(dim));
00259 break;
00260 }
00261 case SPLIT_SAH: default: {
00262
00263 if (nPrimitives <= 4) {
00264
00265 mid = (start + end) / 2;
00266 std::nth_element(&buildData[start], &buildData[mid],
00267 &buildData[end-1]+1, ComparePoints(dim));
00268 }
00269 else {
00270
00271 const int nBuckets = 12;
00272 struct BucketInfo {
00273 BucketInfo() { count = 0; }
00274 int count;
00275 BBox bounds;
00276 };
00277 BucketInfo buckets[nBuckets];
00278
00279
00280 for (uint32_t i = start; i < end; ++i) {
00281 int b = nBuckets *
00282 ((buildData[i].centroid[dim] - centroidBounds.pMin[dim]) /
00283 (centroidBounds.pMax[dim] - centroidBounds.pMin[dim]));
00284 if (b == nBuckets) b = nBuckets-1;
00285 Assert(b >= 0 && b < nBuckets);
00286 buckets[b].count++;
00287 buckets[b].bounds = Union(buckets[b].bounds, buildData[i].bounds);
00288 }
00289
00290
00291 float cost[nBuckets-1];
00292 for (int i = 0; i < nBuckets-1; ++i) {
00293 BBox b0, b1;
00294 int count0 = 0, count1 = 0;
00295 for (int j = 0; j <= i; ++j) {
00296 b0 = Union(b0, buckets[j].bounds);
00297 count0 += buckets[j].count;
00298 }
00299 for (int j = i+1; j < nBuckets; ++j) {
00300 b1 = Union(b1, buckets[j].bounds);
00301 count1 += buckets[j].count;
00302 }
00303 cost[i] = .125f + (count0*b0.SurfaceArea() + count1*b1.SurfaceArea()) /
00304 bbox.SurfaceArea();
00305 }
00306
00307
00308 float minCost = cost[0];
00309 uint32_t minCostSplit = 0;
00310 for (int i = 1; i < nBuckets-1; ++i) {
00311 if (cost[i] < minCost) {
00312 minCost = cost[i];
00313 minCostSplit = i;
00314 }
00315 }
00316
00317
00318 if (nPrimitives > maxPrimsInNode ||
00319 minCost < nPrimitives) {
00320 BVHPrimitiveInfo *pmid = std::partition(&buildData[start],
00321 &buildData[end-1]+1,
00322 CompareToBucket(minCostSplit, nBuckets, dim, centroidBounds));
00323 mid = pmid - &buildData[0];
00324 }
00325
00326 else {
00327
00328 uint32_t firstPrimOffset = orderedPrims.size();
00329 for (uint32_t i = start; i < end; ++i) {
00330 uint32_t primNum = buildData[i].primitiveNumber;
00331 orderedPrims.push_back(primitives[primNum]);
00332 }
00333 node->InitLeaf(firstPrimOffset, nPrimitives, bbox);
00334 }
00335 }
00336 break;
00337 }
00338 }
00339 node->InitInterior(dim,
00340 recursiveBuild(buildArena, buildData, start, mid,
00341 totalNodes, orderedPrims),
00342 recursiveBuild(buildArena, buildData, mid, end,
00343 totalNodes, orderedPrims));
00344 }
00345 return node;
00346 }
00347
00348
00349 uint32_t BVHAccel::flattenBVHTree(BVHBuildNode *node, uint32_t *offset) {
00350 LinearBVHNode *linearNode = &nodes[*offset];
00351 linearNode->bounds = node->bounds;
00352 uint32_t myOffset = (*offset)++;
00353 if (node->nPrimitives > 0) {
00354 Assert(!node->children[0] && !node->children[1]);
00355 linearNode->primitivesOffset = node->firstPrimOffset;
00356 linearNode->nPrimitives = node->nPrimitives;
00357 }
00358 else {
00359
00360 linearNode->axis = node->splitAxis;
00361 linearNode->nPrimitives = 0;
00362 flattenBVHTree(node->children[0], offset);
00363 linearNode->secondChildOffset = flattenBVHTree(node->children[1],
00364 offset);
00365 }
00366 return myOffset;
00367 }
00368
00369
00370 BVHAccel::~BVHAccel() {
00371 FreeAligned(nodes);
00372 }
00373
00374
00375 bool BVHAccel::Intersect(const Ray &ray, Intersection *isect) const {
00376 if (!nodes) return false;
00377 PBRT_BVH_INTERSECTION_STARTED(const_cast<BVHAccel *>(this), const_cast<Ray *>(&ray));
00378 bool hit = false;
00379 Point origin = ray(ray.mint);
00380 Vector invDir(1.f / ray.d.x, 1.f / ray.d.y, 1.f / ray.d.z);
00381 uint32_t dirIsNeg[3] = { invDir.x < 0, invDir.y < 0, invDir.z < 0 };
00382
00383 uint32_t todoOffset = 0, nodeNum = 0;
00384 uint32_t todo[64];
00385 while (true) {
00386 const LinearBVHNode *node = &nodes[nodeNum];
00387
00388 if (::IntersectP(node->bounds, ray, invDir, dirIsNeg)) {
00389 if (node->nPrimitives > 0) {
00390
00391 PBRT_BVH_INTERSECTION_TRAVERSED_LEAF_NODE(const_cast<LinearBVHNode *>(node));
00392 for (uint32_t i = 0; i < node->nPrimitives; ++i)
00393 {
00394 PBRT_BVH_INTERSECTION_PRIMITIVE_TEST(const_cast<Primitive *>(primitives[node->primitivesOffset+i].GetPtr()));
00395 if (primitives[node->primitivesOffset+i]->Intersect(ray, isect))
00396 {
00397 PBRT_BVH_INTERSECTION_PRIMITIVE_HIT(const_cast<Primitive *>(primitives[node->primitivesOffset+i].GetPtr()));
00398 hit = true;
00399 }
00400 else {
00401 PBRT_BVH_INTERSECTION_PRIMITIVE_MISSED(const_cast<Primitive *>(primitives[node->primitivesOffset+i].GetPtr()));
00402 }
00403 }
00404 if (todoOffset == 0) break;
00405 nodeNum = todo[--todoOffset];
00406 }
00407 else {
00408
00409 PBRT_BVH_INTERSECTION_TRAVERSED_INTERIOR_NODE(const_cast<LinearBVHNode *>(node));
00410 if (dirIsNeg[node->axis]) {
00411 todo[todoOffset++] = nodeNum + 1;
00412 nodeNum = node->secondChildOffset;
00413 }
00414 else {
00415 todo[todoOffset++] = node->secondChildOffset;
00416 nodeNum = nodeNum + 1;
00417 }
00418 }
00419 }
00420 else {
00421 if (todoOffset == 0) break;
00422 nodeNum = todo[--todoOffset];
00423 }
00424 }
00425 PBRT_BVH_INTERSECTION_FINISHED();
00426 return hit;
00427 }
00428
00429
00430 bool BVHAccel::IntersectP(const Ray &ray) const {
00431 if (!nodes) return false;
00432 PBRT_BVH_INTERSECTIONP_STARTED(const_cast<BVHAccel *>(this), const_cast<Ray *>(&ray));
00433 Vector invDir(1.f / ray.d.x, 1.f / ray.d.y, 1.f / ray.d.z);
00434 uint32_t dirIsNeg[3] = { invDir.x < 0, invDir.y < 0, invDir.z < 0 };
00435 uint32_t todo[64];
00436 uint32_t todoOffset = 0, nodeNum = 0;
00437 while (true) {
00438 const LinearBVHNode *node = &nodes[nodeNum];
00439 if (::IntersectP(node->bounds, ray, invDir, dirIsNeg)) {
00440
00441 if (node->nPrimitives > 0) {
00442 PBRT_BVH_INTERSECTIONP_TRAVERSED_LEAF_NODE(const_cast<LinearBVHNode *>(node));
00443 for (uint32_t i = 0; i < node->nPrimitives; ++i) {
00444 PBRT_BVH_INTERSECTIONP_PRIMITIVE_TEST(const_cast<Primitive *>(primitives[node->primitivesOffset + i].GetPtr()));
00445 if (primitives[node->primitivesOffset+i]->IntersectP(ray)) {
00446 PBRT_BVH_INTERSECTIONP_PRIMITIVE_HIT(const_cast<Primitive *>(primitives[node->primitivesOffset+i].GetPtr()));
00447 return true;
00448 }
00449 else {
00450 PBRT_BVH_INTERSECTIONP_PRIMITIVE_MISSED(const_cast<Primitive *>(primitives[node->primitivesOffset + i].GetPtr()));
00451 }
00452 }
00453 if (todoOffset == 0) break;
00454 nodeNum = todo[--todoOffset];
00455 }
00456 else {
00457 PBRT_BVH_INTERSECTIONP_TRAVERSED_INTERIOR_NODE(const_cast<LinearBVHNode *>(node));
00458 if (dirIsNeg[node->axis]) {
00460 todo[todoOffset++] = nodeNum + 1;
00461 nodeNum = node->secondChildOffset;
00462 }
00463 else {
00464 todo[todoOffset++] = node->secondChildOffset;
00465 nodeNum = nodeNum + 1;
00466 }
00467 }
00468 }
00469 else {
00470 if (todoOffset == 0) break;
00471 nodeNum = todo[--todoOffset];
00472 }
00473 }
00474 PBRT_BVH_INTERSECTIONP_FINISHED();
00475 return false;
00476 }
00477
00478
00479 BVHAccel *CreateBVHAccelerator(const vector<Reference<Primitive> > &prims,
00480 const ParamSet &ps) {
00481 string splitMethod = ps.FindOneString("splitmethod", "sah");
00482 uint32_t maxPrimsInNode = ps.FindOneInt("maxnodeprims", 4);
00483 return new BVHAccel(prims, maxPrimsInNode, splitMethod);
00484 }
00485
00486