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