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