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/kdtreeaccel.h"
00028 #include "paramset.h"
00029
00030
00031 struct KdAccelNode {
00032
00033 void initLeaf(uint32_t *primNums, int np, MemoryArena &arena);
00034 void initInterior(uint32_t axis, uint32_t ac, float s) {
00035 split = s;
00036 flags = axis;
00037 aboveChild |= (ac << 2);
00038 }
00039 float SplitPos() const { return split; }
00040 uint32_t nPrimitives() const { return nPrims >> 2; }
00041 uint32_t SplitAxis() const { return flags & 3; }
00042 bool IsLeaf() const { return (flags & 3) == 3; }
00043 uint32_t AboveChild() const { return aboveChild >> 2; }
00044 union {
00045 float split;
00046 uint32_t onePrimitive;
00047 uint32_t *primitives;
00048 };
00049
00050 private:
00051 union {
00052 uint32_t flags;
00053 uint32_t nPrims;
00054 uint32_t aboveChild;
00055 };
00056 };
00057
00058
00059 struct BoundEdge {
00060
00061 BoundEdge() { }
00062 BoundEdge(float tt, int pn, bool starting) {
00063 t = tt;
00064 primNum = pn;
00065 type = starting ? START : END;
00066 }
00067 bool operator<(const BoundEdge &e) const {
00068 if (t == e.t)
00069 return (int)type < (int)e.type;
00070 else return t < e.t;
00071 }
00072 float t;
00073 int primNum;
00074 enum { START, END } type;
00075 };
00076
00077
00078
00079
00080 KdTreeAccel::KdTreeAccel(const vector<Reference<Primitive> > &p,
00081 int icost, int tcost, float ebonus, int maxp,
00082 int md)
00083 : isectCost(icost), traversalCost(tcost), maxPrims(maxp), maxDepth(md),
00084 emptyBonus(ebonus) {
00085 PBRT_KDTREE_STARTED_CONSTRUCTION(this, p.size());
00086 for (uint32_t i = 0; i < p.size(); ++i)
00087 p[i]->FullyRefine(primitives);
00088
00089 nextFreeNode = nAllocedNodes = 0;
00090 if (maxDepth <= 0)
00091 maxDepth = Round2Int(8 + 1.3f * Log2Int(float(primitives.size())));
00092
00093
00094 vector<BBox> primBounds;
00095 primBounds.reserve(primitives.size());
00096 for (uint32_t i = 0; i < primitives.size(); ++i) {
00097 BBox b = primitives[i]->WorldBound();
00098 bounds = Union(bounds, b);
00099 primBounds.push_back(b);
00100 }
00101
00102
00103 BoundEdge *edges[3];
00104 for (int i = 0; i < 3; ++i)
00105 edges[i] = new BoundEdge[2*primitives.size()];
00106 uint32_t *prims0 = new uint32_t[primitives.size()];
00107 uint32_t *prims1 = new uint32_t[(maxDepth+1) * primitives.size()];
00108
00109
00110 uint32_t *primNums = new uint32_t[primitives.size()];
00111 for (uint32_t i = 0; i < primitives.size(); ++i)
00112 primNums[i] = i;
00113
00114
00115 buildTree(0, bounds, primBounds, primNums, primitives.size(),
00116 maxDepth, edges, prims0, prims1);
00117
00118
00119 delete[] primNums;
00120 for (int i = 0; i < 3; ++i)
00121 delete[] edges[i];
00122 delete[] prims0;
00123 delete[] prims1;
00124 PBRT_KDTREE_FINISHED_CONSTRUCTION(this);
00125 }
00126
00127
00128 void KdAccelNode::initLeaf(uint32_t *primNums, int np,
00129 MemoryArena &arena) {
00130 flags = 3;
00131 nPrims |= (np << 2);
00132
00133 if (np == 0)
00134 onePrimitive = 0;
00135 else if (np == 1)
00136 onePrimitive = primNums[0];
00137 else {
00138 primitives = arena.Alloc<uint32_t>(np);
00139 for (int i = 0; i < np; ++i)
00140 primitives[i] = primNums[i];
00141 }
00142 }
00143
00144
00145 KdTreeAccel::~KdTreeAccel() {
00146 FreeAligned(nodes);
00147 }
00148
00149
00150 void KdTreeAccel::buildTree(int nodeNum, const BBox &nodeBounds,
00151 const vector<BBox> &allPrimBounds, uint32_t *primNums,
00152 int nPrimitives, int depth, BoundEdge *edges[3],
00153 uint32_t *prims0, uint32_t *prims1, int badRefines) {
00154 Assert(nodeNum == nextFreeNode);
00155
00156 if (nextFreeNode == nAllocedNodes) {
00157 int nAlloc = max(2 * nAllocedNodes, 512);
00158 KdAccelNode *n = AllocAligned<KdAccelNode>(nAlloc);
00159 if (nAllocedNodes > 0) {
00160 memcpy(n, nodes, nAllocedNodes * sizeof(KdAccelNode));
00161 FreeAligned(nodes);
00162 }
00163 nodes = n;
00164 nAllocedNodes = nAlloc;
00165 }
00166 ++nextFreeNode;
00167
00168
00169 if (nPrimitives <= maxPrims || depth == 0) {
00170 PBRT_KDTREE_CREATED_LEAF(nPrimitives, maxDepth-depth);
00171 nodes[nodeNum].initLeaf(primNums, nPrimitives, arena);
00172 return;
00173 }
00174
00175
00176
00177
00178 int bestAxis = -1, bestOffset = -1;
00179 float bestCost = INFINITY;
00180 float oldCost = isectCost * float(nPrimitives);
00181 float totalSA = nodeBounds.SurfaceArea();
00182 float invTotalSA = 1.f / totalSA;
00183 Vector d = nodeBounds.pMax - nodeBounds.pMin;
00184
00185
00186 uint32_t axis = nodeBounds.MaximumExtent();
00187 int retries = 0;
00188 retrySplit:
00189
00190
00191 for (int i = 0; i < nPrimitives; ++i) {
00192 int pn = primNums[i];
00193 const BBox &bbox = allPrimBounds[pn];
00194 edges[axis][2*i] = BoundEdge(bbox.pMin[axis], pn, true);
00195 edges[axis][2*i+1] = BoundEdge(bbox.pMax[axis], pn, false);
00196 }
00197 sort(&edges[axis][0], &edges[axis][2*nPrimitives]);
00198
00199
00200 int nBelow = 0, nAbove = nPrimitives;
00201 for (int i = 0; i < 2*nPrimitives; ++i) {
00202 if (edges[axis][i].type == BoundEdge::END) --nAbove;
00203 float edget = edges[axis][i].t;
00204 if (edget > nodeBounds.pMin[axis] &&
00205 edget < nodeBounds.pMax[axis]) {
00206
00207 uint32_t otherAxis0 = (axis + 1) % 3, otherAxis1 = (axis + 2) % 3;
00208 float belowSA = 2 * (d[otherAxis0] * d[otherAxis1] +
00209 (edget - nodeBounds.pMin[axis]) *
00210 (d[otherAxis0] + d[otherAxis1]));
00211 float aboveSA = 2 * (d[otherAxis0] * d[otherAxis1] +
00212 (nodeBounds.pMax[axis] - edget) *
00213 (d[otherAxis0] + d[otherAxis1]));
00214 float pBelow = belowSA * invTotalSA;
00215 float pAbove = aboveSA * invTotalSA;
00216 float eb = (nAbove == 0 || nBelow == 0) ? emptyBonus : 0.f;
00217 float cost = traversalCost +
00218 isectCost * (1.f - eb) * (pBelow * nBelow + pAbove * nAbove);
00219
00220
00221 if (cost < bestCost) {
00222 bestCost = cost;
00223 bestAxis = axis;
00224 bestOffset = i;
00225 }
00226 }
00227 if (edges[axis][i].type == BoundEdge::START) ++nBelow;
00228 }
00229 Assert(nBelow == nPrimitives && nAbove == 0);
00230
00231
00232 if (bestAxis == -1 && retries < 2) {
00233 ++retries;
00234 axis = (axis+1) % 3;
00235 goto retrySplit;
00236 }
00237 if (bestCost > oldCost) ++badRefines;
00238 if ((bestCost > 4.f * oldCost && nPrimitives < 16) ||
00239 bestAxis == -1 || badRefines == 3) {
00240 PBRT_KDTREE_CREATED_LEAF(nPrimitives, maxDepth-depth);
00241 nodes[nodeNum].initLeaf(primNums, nPrimitives, arena);
00242 return;
00243 }
00244
00245
00246 int n0 = 0, n1 = 0;
00247 for (int i = 0; i < bestOffset; ++i)
00248 if (edges[bestAxis][i].type == BoundEdge::START)
00249 prims0[n0++] = edges[bestAxis][i].primNum;
00250 for (int i = bestOffset+1; i < 2*nPrimitives; ++i)
00251 if (edges[bestAxis][i].type == BoundEdge::END)
00252 prims1[n1++] = edges[bestAxis][i].primNum;
00253
00254
00255 float tsplit = edges[bestAxis][bestOffset].t;
00256 PBRT_KDTREE_CREATED_INTERIOR_NODE(bestAxis, tsplit);
00257 BBox bounds0 = nodeBounds, bounds1 = nodeBounds;
00258 bounds0.pMax[bestAxis] = bounds1.pMin[bestAxis] = tsplit;
00259 buildTree(nodeNum+1, bounds0,
00260 allPrimBounds, prims0, n0, depth-1, edges,
00261 prims0, prims1 + nPrimitives, badRefines);
00262 uint32_t aboveChild = nextFreeNode;
00263 nodes[nodeNum].initInterior(bestAxis, aboveChild, tsplit);
00264 buildTree(aboveChild, bounds1, allPrimBounds, prims1, n1,
00265 depth-1, edges, prims0, prims1 + nPrimitives, badRefines);
00266 }
00267
00268
00269 bool KdTreeAccel::Intersect(const Ray &ray,
00270 Intersection *isect) const {
00271 PBRT_KDTREE_INTERSECTION_TEST(const_cast<KdTreeAccel *>(this), const_cast<Ray *>(&ray));
00272
00273 float tmin, tmax;
00274 if (!bounds.IntersectP(ray, &tmin, &tmax))
00275 {
00276 PBRT_KDTREE_RAY_MISSED_BOUNDS();
00277 return false;
00278 }
00279
00280
00281 Vector invDir(1.f/ray.d.x, 1.f/ray.d.y, 1.f/ray.d.z);
00282 #define MAX_TODO 64
00283 KdToDo todo[MAX_TODO];
00284 int todoPos = 0;
00285
00286
00287 bool hit = false;
00288 const KdAccelNode *node = &nodes[0];
00289 while (node != NULL) {
00290
00291 if (ray.maxt < tmin) break;
00292 if (!node->IsLeaf()) {
00293 PBRT_KDTREE_INTERSECTION_TRAVERSED_INTERIOR_NODE(const_cast<KdAccelNode *>(node));
00294
00295
00296
00297 int axis = node->SplitAxis();
00298 float tplane = (node->SplitPos() - ray.o[axis]) * invDir[axis];
00299
00300
00301 const KdAccelNode *firstChild, *secondChild;
00302 int belowFirst = (ray.o[axis] < node->SplitPos()) ||
00303 (ray.o[axis] == node->SplitPos() && ray.d[axis] >= 0);
00304 if (belowFirst) {
00305 firstChild = node + 1;
00306 secondChild = &nodes[node->AboveChild()];
00307 }
00308 else {
00309 firstChild = &nodes[node->AboveChild()];
00310 secondChild = node + 1;
00311 }
00312
00313
00314 if (tplane > tmax || tplane <= 0)
00315 node = firstChild;
00316 else if (tplane < tmin)
00317 node = secondChild;
00318 else {
00319
00320 todo[todoPos].node = secondChild;
00321 todo[todoPos].tmin = tplane;
00322 todo[todoPos].tmax = tmax;
00323 ++todoPos;
00324 node = firstChild;
00325 tmax = tplane;
00326 }
00327 }
00328 else {
00329 PBRT_KDTREE_INTERSECTION_TRAVERSED_LEAF_NODE(const_cast<KdAccelNode *>(node), node->nPrimitives());
00330
00331 uint32_t nPrimitives = node->nPrimitives();
00332 if (nPrimitives == 1) {
00333 const Reference<Primitive> &prim = primitives[node->onePrimitive];
00334
00335 PBRT_KDTREE_INTERSECTION_PRIMITIVE_TEST(const_cast<Primitive *>(prim.GetPtr()));
00336 if (prim->Intersect(ray, isect))
00337 {
00338 PBRT_KDTREE_INTERSECTION_HIT(const_cast<Primitive *>(prim.GetPtr()));
00339 hit = true;
00340 }
00341 }
00342 else {
00343 uint32_t *prims = node->primitives;
00344 for (uint32_t i = 0; i < nPrimitives; ++i) {
00345 const Reference<Primitive> &prim = primitives[prims[i]];
00346
00347 PBRT_KDTREE_INTERSECTION_PRIMITIVE_TEST(const_cast<Primitive *>(prim.GetPtr()));
00348 if (prim->Intersect(ray, isect))
00349 {
00350 PBRT_KDTREE_INTERSECTION_HIT(const_cast<Primitive *>(prim.GetPtr()));
00351 hit = true;
00352 }
00353 }
00354 }
00355
00356
00357 if (todoPos > 0) {
00358 --todoPos;
00359 node = todo[todoPos].node;
00360 tmin = todo[todoPos].tmin;
00361 tmax = todo[todoPos].tmax;
00362 }
00363 else
00364 break;
00365 }
00366 }
00367 PBRT_KDTREE_INTERSECTION_FINISHED();
00368 return hit;
00369 }
00370
00371
00372 bool KdTreeAccel::IntersectP(const Ray &ray) const {
00373 PBRT_KDTREE_INTERSECTIONP_TEST(const_cast<KdTreeAccel *>(this), const_cast<Ray *>(&ray));
00374
00375 float tmin, tmax;
00376 if (!bounds.IntersectP(ray, &tmin, &tmax))
00377 {
00378 PBRT_KDTREE_RAY_MISSED_BOUNDS();
00379 return false;
00380 }
00381
00382
00383 Vector invDir(1.f/ray.d.x, 1.f/ray.d.y, 1.f/ray.d.z);
00384 #define MAX_TODO 64
00385 KdToDo todo[MAX_TODO];
00386 int todoPos = 0;
00387 const KdAccelNode *node = &nodes[0];
00388 while (node != NULL) {
00389 if (node->IsLeaf()) {
00390 PBRT_KDTREE_INTERSECTIONP_TRAVERSED_LEAF_NODE(const_cast<KdAccelNode *>(node), node->nPrimitives());
00391
00392 uint32_t nPrimitives = node->nPrimitives();
00393 if (nPrimitives == 1) {
00394 const Reference<Primitive> &prim = primitives[node->onePrimitive];
00395 PBRT_KDTREE_INTERSECTIONP_PRIMITIVE_TEST(const_cast<Primitive *>(prim.GetPtr()));
00396 if (prim->IntersectP(ray)) {
00397 PBRT_KDTREE_INTERSECTIONP_HIT(const_cast<Primitive *>(prim.GetPtr()));
00398 return true;
00399 }
00400 }
00401 else {
00402 uint32_t *prims = node->primitives;
00403 for (uint32_t i = 0; i < nPrimitives; ++i) {
00404 const Reference<Primitive> &prim = primitives[prims[i]];
00405 PBRT_KDTREE_INTERSECTIONP_PRIMITIVE_TEST(const_cast<Primitive *>(prim.GetPtr()));
00406 if (prim->IntersectP(ray)) {
00407 PBRT_KDTREE_INTERSECTIONP_HIT(const_cast<Primitive *>(prim.GetPtr()));
00408 return true;
00409 }
00410 }
00411 }
00412
00413
00414 if (todoPos > 0) {
00415 --todoPos;
00416 node = todo[todoPos].node;
00417 tmin = todo[todoPos].tmin;
00418 tmax = todo[todoPos].tmax;
00419 }
00420 else
00421 break;
00422 }
00423 else {
00424 PBRT_KDTREE_INTERSECTIONP_TRAVERSED_INTERIOR_NODE(const_cast<KdAccelNode *>(node));
00425
00426
00427
00428 int axis = node->SplitAxis();
00429 float tplane = (node->SplitPos() - ray.o[axis]) * invDir[axis];
00430
00431
00432 const KdAccelNode *firstChild, *secondChild;
00433 int belowFirst = (ray.o[axis] < node->SplitPos()) ||
00434 (ray.o[axis] == node->SplitPos() && ray.d[axis] >= 0);
00435 if (belowFirst) {
00436 firstChild = node + 1;
00437 secondChild = &nodes[node->AboveChild()];
00438 }
00439 else {
00440 firstChild = &nodes[node->AboveChild()];
00441 secondChild = node + 1;
00442 }
00443
00444
00445 if (tplane > tmax || tplane <= 0)
00446 node = firstChild;
00447 else if (tplane < tmin)
00448 node = secondChild;
00449 else {
00450
00451 todo[todoPos].node = secondChild;
00452 todo[todoPos].tmin = tplane;
00453 todo[todoPos].tmax = tmax;
00454 ++todoPos;
00455 node = firstChild;
00456 tmax = tplane;
00457 }
00458 }
00459 }
00460 PBRT_KDTREE_INTERSECTIONP_MISSED();
00461 return false;
00462 }
00463
00464
00465 KdTreeAccel *CreateKdTreeAccelerator(const vector<Reference<Primitive> > &prims,
00466 const ParamSet &ps) {
00467 int isectCost = ps.FindOneInt("intersectcost", 80);
00468 int travCost = ps.FindOneInt("traversalcost", 1);
00469 float emptyBonus = ps.FindOneFloat("emptybonus", 0.5f);
00470 int maxPrims = ps.FindOneInt("maxprims", 1);
00471 int maxDepth = ps.FindOneInt("maxdepth", -1);
00472 return new KdTreeAccel(prims, isectCost, travCost,
00473 emptyBonus, maxPrims, maxDepth);
00474 }
00475
00476