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