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 "waldtrianglemesh.h"
00025 #include "mc.h"
00026 #include "paramset.h"
00027
00028 using namespace lux;
00029
00030
00031 WaldTriangleMesh::WaldTriangleMesh(const Transform &o2w, bool ro,
00032 int nt, int nv, const int *vi, const Point *P,
00033 const Normal *N, const Vector *S, const float *uv)
00034 : Shape(o2w, ro) {
00035 ntris = nt;
00036 nverts = nv;
00037 vertexIndex = new int[3 * ntris];
00038 memcpy(vertexIndex, vi, 3 * ntris * sizeof(int));
00039
00040
00041 if (uv) {
00042 uvs = new float[2*nverts];
00043 memcpy(uvs, uv, 2*nverts*sizeof(float));
00044 } else
00045 uvs = NULL;
00046
00047 p = new Point[nverts];
00048 if (N) {
00049 n = new Normal[nverts];
00050 memcpy(n, N, nverts*sizeof(Normal));
00051 } else
00052 n = NULL;
00053
00054 if (S) {
00055 s = new Vector[nverts];
00056 memcpy(s, S, nverts*sizeof(Vector));
00057 } else
00058 s = NULL;
00059
00060
00061 for (int i = 0; i < nverts; ++i)
00062 p[i] = ObjectToWorld(P[i]);
00063 }
00064
00065 WaldTriangleMesh::~WaldTriangleMesh() {
00066 delete[] vertexIndex;
00067 delete[] p;
00068 delete[] s;
00069 delete[] n;
00070 delete[] uvs;
00071 }
00072
00073 BBox WaldTriangleMesh::ObjectBound() const {
00074 BBox bobj;
00075 for (int i = 0; i < nverts; i++)
00076 bobj = Union(bobj, WorldToObject(p[i]));
00077 return bobj;
00078 }
00079
00080 BBox WaldTriangleMesh::WorldBound() const {
00081 BBox worldBounds;
00082 for (int i = 0; i < nverts; i++)
00083 worldBounds = Union(worldBounds, p[i]);
00084 return worldBounds;
00085 }
00086
00087 void
00088 WaldTriangleMesh::Refine(vector<boost::shared_ptr<Shape> > &refined)
00089 const {
00090 for (int i = 0; i < ntris; ++i) {
00091 boost::shared_ptr<Shape> o(new WaldTriangle(ObjectToWorld,
00092 reverseOrientation,
00093 (WaldTriangleMesh *)this,
00094 i));
00095 refined.push_back(o);
00096 }
00097 }
00098
00099 WaldTriangle::WaldTriangle(const Transform &o2w, bool ro,
00100 WaldTriangleMesh *m, int n)
00101 : Shape(o2w, ro) {
00102 mesh = m;
00103 v = &mesh->vertexIndex[3*n];
00104
00105
00106
00107
00108
00109
00110
00111
00112 const Point &v0 = mesh->p[v[0]];
00113 const Point &v1 = mesh->p[v[1]];
00114 const Point &v2 = mesh->p[v[2]];
00115 Vector e1 = v1 - v0;
00116 Vector e2 = v2 - v0;
00117
00118 Vector normal = Normalize(Cross(e1, e2));
00119
00120 if (isnan(normal.x) || isnan(normal.y) || isnan(normal.z)) {
00121 intersectionType = DEGENERATE;
00122 return;
00123 }
00124
00125
00126
00127
00128 if ((normal.y == 0.0f) && (normal.z == 0.0f))
00129 intersectionType = ORTHOGONAL_X;
00130 else if((normal.x == 0.0f) && (normal.z == 0.0f))
00131 intersectionType = ORTHOGONAL_Y;
00132 else if((normal.x == 0.0f) && (normal.y == 0.0f))
00133 intersectionType = ORTHOGONAL_Z;
00134 else if ((fabs(normal.x) > fabs(normal.y)) && (fabs(normal.x) > fabs(normal.z)))
00135 intersectionType = DOMINANT_X;
00136 else if (fabs(normal.y) > fabs(normal.z))
00137 intersectionType = DOMINANT_Y;
00138 else
00139 intersectionType = DOMINANT_Z;
00140
00141 float ax, ay, bx, by, cx, cy;
00142 switch (intersectionType) {
00143 case DOMINANT_X: {
00144 const float invNormal = 1.0f / normal.x;
00145 nu = normal.y * invNormal;
00146 nv = normal.z * invNormal;
00147 nd = v0.x + nu*v0.y + nv * v0.z;
00148 ax = v0.y;
00149 ay = v0.z;
00150 bx = v2.y - ax;
00151 by = v2.z - ay;
00152 cx = v1.y - ax;
00153 cy = v1.z - ay;
00154 break;
00155 }
00156 case DOMINANT_Y: {
00157 const float invNormal = 1.0f / normal.y;
00158 nu = normal.z * invNormal;
00159 nv = normal.x * invNormal;
00160 nd = nv * v0.x + v0.y + nu * v0.z;
00161 ax = v0.z;
00162 ay = v0.x;
00163 bx = v2.z - ax;
00164 by = v2.x - ay;
00165 cx = v1.z - ax;
00166 cy = v1.x - ay;
00167 break;
00168 }
00169 case DOMINANT_Z: {
00170 const float invNormal = 1.0f / normal.z;
00171 nu = normal.x * invNormal;
00172 nv = normal.y * invNormal;
00173 nd = nu * v0.x + nv*v0.y + v0.z;
00174 ax = v0.x;
00175 ay = v0.y;
00176 bx = v2.x - ax;
00177 by = v2.y - ay;
00178 cx = v1.x - ax;
00179 cy = v1.y - ay;
00180 break;
00181 }
00182 case ORTHOGONAL_X:
00183 nu = 0.0f;
00184 nv = 0.0f;
00185 nd = v0.x;
00186 ax = v0.y;
00187 ay = v0.z;
00188 bx = v2.y - ax;
00189 by = v2.z - ay;
00190 cx = v1.y - ax;
00191 cy = v1.z - ay;
00192 break;
00193 case ORTHOGONAL_Y:
00194 nu = 0.0f;
00195 nv = 0.0f;
00196 nd = v0.y;
00197 ax = v0.z;
00198 ay = v0.x;
00199 bx = v2.z - ax;
00200 by = v2.x - ay;
00201 cx = v1.z - ax;
00202 cy = v1.x - ay;
00203 break;
00204 case ORTHOGONAL_Z:
00205 nu = 0.0f;
00206 nv = 0.0f;
00207 nd = v0.z;
00208 ax = v0.x;
00209 ay = v0.y;
00210 bx = v2.x - ax;
00211 by = v2.y - ay;
00212 cx = v1.x - ax;
00213 cy = v1.y - ay;
00214 break;
00215 default:
00216 BOOST_ASSERT(false);
00217
00218 return;
00219 }
00220
00221 float det = bx * cy - by * cx;
00222 float invDet = 1.0f / det;
00223
00224 bnu = -by * invDet;
00225 bnv = bx * invDet;
00226 bnd = (by * ax - bx * ay) * invDet;
00227 cnu = cy * invDet;
00228 cnv = -cx * invDet;
00229 cnd = (cx * ay - cy * ax) * invDet;
00230
00231
00232
00233
00234
00235 float uvs[3][2];
00236 GetUVs(uvs);
00237
00238 const float du1 = uvs[0][0] - uvs[2][0];
00239 const float du2 = uvs[1][0] - uvs[2][0];
00240 const float dv1 = uvs[0][1] - uvs[2][1];
00241 const float dv2 = uvs[1][1] - uvs[2][1];
00242 const Vector dp1 = v0 - v2, dp2 = v1 - v2;
00243 const float determinant = du1 * dv2 - dv1 * du2;
00244 if (determinant == 0.f) {
00245
00246 CoordinateSystem(Normalize(Cross(e1, e2)), &dpdu, &dpdv);
00247 } else {
00248 const float invdet = 1.f / determinant;
00249 dpdu = ( dv2 * dp1 - dv1 * dp2) * invdet;
00250 dpdv = (-du2 * dp1 + du1 * dp2) * invdet;
00251 }
00252
00253
00254
00255
00256 normalizedNormal = Normal(Normalize(Cross(dpdu, dpdv)));
00257 if(mesh->n) {
00258 if(Dot(ObjectToWorld(mesh->n[v[0]]+mesh->n[v[1]]+mesh->n[v[2]]), normalizedNormal) < 0)
00259 normalizedNormal *= -1;
00260 } else {
00261 if(Dot(Cross(e1, e2), normalizedNormal) < 0)
00262 normalizedNormal *= -1;
00263 }
00264
00265
00266 if (this->reverseOrientation ^ this->transformSwapsHandedness)
00267 normalizedNormal *= -1.f;
00268 }
00269
00270 BBox WaldTriangle::ObjectBound() const {
00271
00272 const Point &p1 = mesh->p[v[0]];
00273 const Point &p2 = mesh->p[v[1]];
00274 const Point &p3 = mesh->p[v[2]];
00275 return Union(BBox(WorldToObject(p1), WorldToObject(p2)),
00276 WorldToObject(p3));
00277 }
00278
00279 BBox WaldTriangle::WorldBound() const {
00280
00281 const Point &p1 = mesh->p[v[0]];
00282 const Point &p2 = mesh->p[v[1]];
00283 const Point &p3 = mesh->p[v[2]];
00284 return Union(BBox(p1, p2), p3);
00285 }
00286
00287 bool WaldTriangle::Intersect(const Ray &ray, float *tHit,
00288 DifferentialGeometry *dg) const {
00289
00290
00291
00292
00293
00294 float uu, vv, t;
00295 switch (intersectionType) {
00296 case DOMINANT_X: {
00297 const float det = ray.d.x + nu * ray.d.y + nv * ray.d.z;
00298 if(det==0.0f)
00299 return false;
00300
00301 const float invDet = 1.0f / det;
00302 t = (nd - ray.o.x - nu * ray.o.y - nv * ray.o.z) * invDet;
00303
00304
00305
00306
00307
00308
00309 if (t < ray.mint || t > ray.maxt)
00310 return false;
00311
00312 const float hu = ray.o.y + t * ray.d.y;
00313 const float hv = ray.o.z + t * ray.d.z;
00314 uu = hu * bnu + hv * bnv + bnd;
00315
00316 if (uu < 0.0f)
00317 return false;
00318
00319 vv = hu * cnu + hv * cnv + cnd;
00320
00321 if (vv < 0.0f)
00322 return false;
00323 if (uu + vv > 1.0f)
00324 return false;
00325 break;
00326 }
00327 case DOMINANT_Y: {
00328 const float det = ray.d.y + nu * ray.d.z + nv * ray.d.x;
00329 if(det==0.0f)
00330 return false;
00331
00332 const float invDet = 1.0f / det;
00333 t = (nd - ray.o.y - nu * ray.o.z - nv * ray.o.x) * invDet;
00334
00335 if (t < ray.mint || t > ray.maxt)
00336 return false;
00337
00338 const float hu = ray.o.z + t * ray.d.z;
00339 const float hv = ray.o.x + t * ray.d.x;
00340 uu = hu * bnu + hv * bnv + bnd;
00341
00342 if (uu < 0.0f)
00343 return false;
00344
00345 vv = hu * cnu + hv * cnv + cnd;
00346
00347 if (vv < 0.0f)
00348 return false;
00349 if (uu + vv > 1.0f)
00350 return false;
00351 break;
00352 }
00353 case DOMINANT_Z: {
00354 const float det = ray.d.z + nu * ray.d.x + nv * ray.d.y;
00355 if(det==0.0f)
00356 return false;
00357
00358 const float invDet = 1.0f / det;
00359 t = (nd - ray.o.z - nu * ray.o.x - nv * ray.o.y) * invDet;
00360
00361 if (t < ray.mint || t > ray.maxt)
00362 return false;
00363
00364 const float hu = ray.o.x + t * ray.d.x;
00365 const float hv = ray.o.y + t * ray.d.y;
00366
00367 uu = hu * bnu + hv * bnv + bnd;
00368
00369 if (uu < 0.0f)
00370 return false;
00371
00372 vv = hu * cnu + hv * cnv + cnd;
00373
00374 if (vv < 0.0f)
00375 return false;
00376 if (uu + vv > 1.0f)
00377 return false;
00378 break;
00379 }
00380 case ORTHOGONAL_X: {
00381 if(ray.d.x == 0.0f)
00382 return false;
00383
00384 const float invDet = 1.0f / ray.d.x;
00385 t = (nd - ray.o.x) * invDet;
00386
00387 if (t < ray.mint || t > ray.maxt)
00388 return false;
00389
00390 const float hu = ray.o.y + t * ray.d.y;
00391 const float hv = ray.o.z + t * ray.d.z;
00392 uu = hu * bnu + hv * bnv + bnd;
00393
00394 if (uu < 0.0f)
00395 return false;
00396
00397 vv = hu * cnu + hv * cnv + cnd;
00398
00399 if (vv < 0.0f)
00400 return false;
00401 if (uu + vv > 1.0f)
00402 return false;
00403 break;
00404 }
00405 case ORTHOGONAL_Y: {
00406 if(ray.d.y == 0.0f)
00407 return false;
00408
00409 const float invDet = 1.0f / ray.d.y;
00410 t = (nd - ray.o.y) * invDet;
00411
00412 if (t < ray.mint || t > ray.maxt)
00413 return false;
00414
00415 const float hu = ray.o.z + t * ray.d.z;
00416 const float hv = ray.o.x + t * ray.d.x;
00417 uu = hu * bnu + hv * bnv + bnd;
00418
00419 if (uu < 0.0f)
00420 return false;
00421
00422 vv = hu * cnu + hv * cnv + cnd;
00423
00424 if (vv < 0.0f)
00425 return false;
00426 if (uu + vv > 1.0f)
00427 return false;
00428 break;
00429 }
00430 case ORTHOGONAL_Z: {
00431 if(ray.d.z == 0.0f)
00432 return false;
00433
00434 const float invDet = 1.0f / ray.d.z;
00435 t = (nd - ray.o.z) * invDet;
00436
00437 if (t < ray.mint || t > ray.maxt)
00438 return false;
00439
00440 const float hu = ray.o.x + t * ray.d.x;
00441 const float hv = ray.o.y + t * ray.d.y;
00442
00443 uu = hu * bnu + hv * bnv + bnd;
00444
00445 if (uu < 0.0f)
00446 return false;
00447
00448 vv = hu * cnu + hv * cnv + cnd;
00449
00450 if (vv < 0.0f)
00451 return false;
00452 if (uu + vv > 1.0f)
00453 return false;
00454 break;
00455 }
00456 case DEGENERATE:
00457 return false;
00458 default:
00459 BOOST_ASSERT(false);
00460
00461 return false;
00462 }
00463
00464
00465 float uvs[3][2];
00466 GetUVs(uvs);
00467
00468 const float b0 = 1.0f - uu - vv;
00469 const float tu = b0*uvs[0][0] + uu*uvs[1][0] + vv*uvs[2][0];
00470 const float tv = b0*uvs[0][1] + uu*uvs[1][1] + vv*uvs[2][1];
00471 *dg = DifferentialGeometry(ray(t),
00472 normalizedNormal,
00473 dpdu, dpdv,
00474 Vector(0, 0, 0), Vector(0, 0, 0),
00475 tu, tv, this);
00476
00477 *tHit = t;
00478 return true;
00479 }
00480
00481 bool WaldTriangle::IntersectP(const Ray &ray) const {
00482 float uu, vv, t;
00483 switch (intersectionType) {
00484 case DOMINANT_X: {
00485 const float det = ray.d.x + nu * ray.d.y + nv * ray.d.z;
00486 if(det==0.0f)
00487 return false;
00488
00489 const float invDet = 1.0f / det;
00490 t = (nd - ray.o.x - nu * ray.o.y - nv * ray.o.z) * invDet;
00491
00492 if (t < ray.mint || t > ray.maxt)
00493 return false;
00494
00495 const float hu = ray.o.y + t * ray.d.y;
00496 const float hv = ray.o.z + t * ray.d.z;
00497 uu = hu * bnu + hv * bnv + bnd;
00498
00499 if (uu < 0.0f)
00500 return false;
00501
00502 vv = hu * cnu + hv * cnv + cnd;
00503
00504 if (vv < 0.0f)
00505 return false;
00506 if (uu + vv > 1.0f)
00507 return false;
00508 break;
00509 }
00510 case DOMINANT_Y: {
00511 const float det = ray.d.y + nu * ray.d.z + nv * ray.d.x;
00512 if(det==0.0f)
00513 return false;
00514
00515 const float invDet = 1.0f / det;
00516 t = (nd - ray.o.y - nu * ray.o.z - nv * ray.o.x) * invDet;
00517
00518 if (t < ray.mint || t > ray.maxt)
00519 return false;
00520
00521 const float hu = ray.o.z + t * ray.d.z;
00522 const float hv = ray.o.x + t * ray.d.x;
00523 uu = hu * bnu + hv * bnv + bnd;
00524
00525 if (uu < 0.0f)
00526 return false;
00527
00528 vv = hu * cnu + hv * cnv + cnd;
00529
00530 if (vv < 0.0f)
00531 return false;
00532 if (uu + vv > 1.0f)
00533 return false;
00534 break;
00535 }
00536 case DOMINANT_Z: {
00537 const float det = ray.d.z + nu * ray.d.x + nv * ray.d.y;
00538 if(det==0.0f)
00539 return false;
00540
00541 const float invDet = 1.0f / det;
00542 t = (nd - ray.o.z - nu * ray.o.x - nv * ray.o.y) * invDet;
00543
00544 if (t < ray.mint || t > ray.maxt)
00545 return false;
00546
00547 const float hu = ray.o.x + t * ray.d.x;
00548 const float hv = ray.o.y + t * ray.d.y;
00549
00550 uu = hu * bnu + hv * bnv + bnd;
00551
00552 if (uu < 0.0f)
00553 return false;
00554
00555 vv = hu * cnu + hv * cnv + cnd;
00556
00557 if (vv < 0.0f)
00558 return false;
00559 if (uu + vv > 1.0f)
00560 return false;
00561 break;
00562 }
00563 case ORTHOGONAL_X: {
00564 if(ray.d.x == 0.0f)
00565 return false;
00566
00567 const float invDet = 1.0f / ray.d.x;
00568 t = (nd - ray.o.x) * invDet;
00569
00570 if (t < ray.mint || t > ray.maxt)
00571 return false;
00572
00573 const float hu = ray.o.y + t * ray.d.y;
00574 const float hv = ray.o.z + t * ray.d.z;
00575 uu = hu * bnu + hv * bnv + bnd;
00576
00577 if (uu < 0.0f)
00578 return false;
00579
00580 vv = hu * cnu + hv * cnv + cnd;
00581
00582 if (vv < 0.0f)
00583 return false;
00584 if (uu + vv > 1.0f)
00585 return false;
00586 break;
00587 }
00588 case ORTHOGONAL_Y: {
00589 if(ray.d.y == 0.0f)
00590 return false;
00591
00592 const float invDet = 1.0f / ray.d.y;
00593 t = (nd - ray.o.y) * invDet;
00594
00595 if (t < ray.mint || t > ray.maxt)
00596 return false;
00597
00598 const float hu = ray.o.z + t * ray.d.z;
00599 const float hv = ray.o.x + t * ray.d.x;
00600 uu = hu * bnu + hv * bnv + bnd;
00601
00602 if (uu < 0.0f)
00603 return false;
00604
00605 vv = hu * cnu + hv * cnv + cnd;
00606
00607 if (vv < 0.0f)
00608 return false;
00609 if (uu + vv > 1.0f)
00610 return false;
00611 break;
00612 }
00613 case ORTHOGONAL_Z: {
00614 if(ray.d.z == 0.0f)
00615 return false;
00616
00617 const float invDet = 1.0f / ray.d.z;
00618 t = (nd - ray.o.z) * invDet;
00619
00620 if (t < ray.mint || t > ray.maxt)
00621 return false;
00622
00623 const float hu = ray.o.x + t * ray.d.x;
00624 const float hv = ray.o.y + t * ray.d.y;
00625
00626 uu = hu * bnu + hv * bnv + bnd;
00627
00628 if (uu < 0.0f)
00629 return false;
00630
00631 vv = hu * cnu + hv * cnv + cnd;
00632
00633 if (vv < 0.0f)
00634 return false;
00635 if (uu + vv > 1.0f)
00636 return false;
00637 break;
00638 }
00639 case DEGENERATE:
00640 return false;
00641 default:
00642 BOOST_ASSERT(false);
00643
00644 return false;
00645 }
00646
00647 return true;
00648 }
00649
00650 float WaldTriangle::Area() const {
00651
00652 const Point &p1 = mesh->p[v[0]];
00653 const Point &p2 = mesh->p[v[1]];
00654 const Point &p3 = mesh->p[v[2]];
00655 return 0.5f * Cross(p2-p1, p3-p1).Length();
00656 }
00657
00658 Point WaldTriangle::Sample(float u1, float u2,
00659 Normal *Ns) const {
00660 float b1, b2;
00661 UniformSampleTriangle(u1, u2, &b1, &b2);
00662
00663 const Point &p1 = mesh->p[v[0]];
00664 const Point &p2 = mesh->p[v[1]];
00665 const Point &p3 = mesh->p[v[2]];
00666 Point p = b1 * p1 + b2 * p2 + (1.f - b1 - b2) * p3;
00667
00668 *Ns = normalizedNormal;
00669
00670 return p;
00671 }
00672
00673 Shape* WaldTriangleMesh::CreateShape(const Transform &o2w,
00674 bool reverseOrientation, const ParamSet ¶ms) {
00675 int nvi, npi, nuvi, nsi, nni;
00676 const int *vi = params.FindInt("indices", &nvi);
00677 const Point *P = params.FindPoint("P", &npi);
00678 const float *uvs = params.FindFloat("uv", &nuvi);
00679
00680 if (!uvs) uvs = params.FindFloat("st", &nuvi);
00681
00682 if (uvs && nuvi != npi * 2) {
00683 luxError(LUX_CONSISTENCY, LUX_ERROR, "Number of \"uv\"s for triangle mesh must match \"P\"s");
00684 uvs = NULL;
00685 }
00686 if (!vi || !P) return NULL;
00687
00688 const Vector *S = params.FindVector("S", &nsi);
00689 if (S && nsi != npi) {
00690 luxError(LUX_CONSISTENCY, LUX_ERROR, "Number of \"S\"s for triangle mesh must match \"P\"s");
00691 S = NULL;
00692 }
00693
00694 const Normal *N = params.FindNormal("N", &nni);
00695 if (N && nni != npi) {
00696 luxError(LUX_CONSISTENCY, LUX_ERROR, "Number of \"N\"s for triangle mesh must match \"P\"s");
00697 N = NULL;
00698 }
00699 if (uvs && N) {
00700
00701
00702 const int *vp = vi;
00703 for (int i = 0; i < nvi; i += 3, vp += 3) {
00704 float area = .5f * Cross(P[vp[0]]-P[vp[1]], P[vp[2]]-P[vp[1]]).Length();
00705 if (area < 1e-7) continue;
00706
00707
00708
00709
00710
00711
00712
00713
00714
00715
00716
00717
00718
00719 }
00720 }
00721 for (int i = 0; i < nvi; ++i)
00722 if (vi[i] >= npi) {
00723
00724
00725 std::stringstream ss;
00726 ss<<"waldtrianglemesh has out of-bounds vertex index "<<vi[i]<<" ("<<npi<<" \"P\" values were given";
00727 luxError(LUX_CONSISTENCY, LUX_ERROR, ss.str().c_str());
00728 return NULL;
00729 }
00730
00731 return new WaldTriangleMesh(o2w, reverseOrientation, nvi/3, npi, vi, P,
00732 N, S, uvs);
00733 }