00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023 #include "mesh.h"
00024 #include "geometry/matrix3x3.h"
00025
00026 using namespace lux;
00027
00028
00029 bool MeshQuadrilateral::IsDegenerate(const Point &p0, const Point &p1, const Point &p2, const Point &p3) {
00030
00031 Vector e0 = p1 - p0;
00032 Vector e1 = p2 - p1;
00033 Vector e2 = p3 - p2;
00034 Vector e3 = p0 - p3;
00035
00036 float el0 = e0.Length();
00037 float el1 = e1.Length();
00038 float el2 = e2.Length();
00039 float el3 = e3.Length();
00040
00041 return el0 < 1e-30 || el1 < 1e-30 || el2 < 1e-30 || el3 < 1e-30;
00042 }
00043
00044
00045
00046 bool MeshQuadrilateral::IsPlanar(const Point &p0, const Point &p1, const Point &p2, const Point &p3) {
00047
00048
00049 Vector e0 = p1 - p0;
00050 Vector e1 = p2 - p0;
00051
00052 Point p = p3;
00053
00054 if (1.f - fabsf(Dot(e0, e1)) < 1e-6) {
00055
00056 e1 = p3 - p0;
00057 p = p2;
00058 }
00059
00060 Vector n = Cross(e1, e0);
00061
00062 Vector x = p - p0;
00063
00064
00065 float D = fabsf(Dot(x, n));
00066
00067
00068
00069 return D < n.Length() * 1e-6f;
00070 }
00071
00072
00073 bool MeshQuadrilateral::IsConvex(const Point &p0, const Point &p1, const Point &p2, const Point &p3) {
00074
00075
00076 Vector b0 = Normalize(p1 - p0);
00077 Vector b1 = p3 - p0;
00078
00079 b1 = Normalize(b1 - b0 * Dot(b1, b0));
00080
00081 if (1.f - fabsf(Dot(b0, b1)) < 1e-6) {
00082
00083 b1 = p2 - p0;
00084
00085 b1 = Normalize(b1 - b0 * Dot(b1, b0));
00086 }
00087
00088
00089 Vector e[4];
00090
00091 e[0] = p1 - p0;
00092 e[1] = p2 - p1;
00093 e[2] = p3 - p2;
00094 e[3] = p0 - p3;
00095
00096
00097 for (int i = 0; i < 4; i++)
00098 e[i] = Vector(Dot(e[i], b0), Dot(e[i], b1), 0);
00099
00100
00101
00102 int altCount = 0;
00103 int curd, prevd;
00104
00105
00106
00107
00108 curd = 1;
00109 for (int i = 1; i <= 4; i++) {
00110 prevd = curd;
00111
00112
00113 curd = (e[i & 3].x < 1e-6) ? (e[i & 3].x > -1e-6f ? prevd : -1) : 1;
00114 altCount += prevd != curd ? 1 : 0;
00115 }
00116
00117 if (altCount != 2)
00118 return false;
00119
00120
00121
00122
00123 int curs, prevs;
00124 altCount = 0;
00125
00126 curs = (Cross(e[1], e[0]).z < 0) ? -1 : 1;
00127 for (int i = 1; i < 4; i++) {
00128 prevs = curs;
00129 curs = (Cross(e[(i + 1) & 3], e[i]).z < 0) ? -1 : 1;
00130 altCount += prevs != curs ? 1 : 0;
00131 }
00132
00133 return altCount == 0;
00134 }
00135
00136 int MeshQuadrilateral::MajorAxis(const Vector &v) {
00137 float absVx = fabsf(v.x);
00138 float absVy = fabsf(v.y);
00139 float absVz = fabsf(v.z);
00140
00141 if (absVx > absVy)
00142 return (absVx > absVz) ? 0 : 2;
00143 return (absVy > absVz) ? 1 : 2;
00144 }
00145
00146 void MeshQuadrilateral::ComputeV11BarycentricCoords(const Vector &e01,
00147 const Vector &e02, const Vector &e03, float *a11, float *b11) {
00148 const Vector N = Cross(e01, e03);
00149
00150 int Nma = MajorAxis(N);
00151
00152 switch (Nma) {
00153 case 0: {
00154 float iNx = 1.f / N.x;
00155 *a11 = (e02.y * e03.z - e02.z * e03.y) * iNx;
00156 *b11 = (e01.y * e02.z - e01.z * e02.y) * iNx;
00157 break;
00158 }
00159 case 1: {
00160 float iNy = 1.f / N.y;
00161 *a11 = (e02.z * e03.x - e02.x * e03.z) * iNy;
00162 *b11 = (e01.z * e02.x - e01.x * e02.z) * iNy;
00163 break;
00164 }
00165 case 2: {
00166 float iNz = 1.f / N.z;
00167 *a11 = (e02.x * e03.y - e02.y * e03.x) * iNz;
00168 *b11 = (e01.x * e02.y - e01.y * e02.x) * iNz;
00169 break;
00170 }
00171 default:
00172 BOOST_ASSERT(false);
00173
00174
00175 break;
00176 }
00177 }
00178
00179
00180 MeshQuadrilateral::MeshQuadrilateral(const Mesh *m, int n)
00181 : mesh(m), idx(&(mesh->quadVertexIndex[4 * n]))
00182 {
00183
00184 const Point &p0 = mesh->WorldToObject(mesh->p[idx[0]]);
00185 const Point &p1 = mesh->WorldToObject(mesh->p[idx[1]]);
00186 const Point &p2 = mesh->WorldToObject(mesh->p[idx[2]]);
00187 const Point &p3 = mesh->WorldToObject(mesh->p[idx[3]]);
00188
00189 if (IsDegenerate(p0, p1, p2, p3)) {
00190 luxError(LUX_CONSISTENCY, LUX_ERROR, string("Degenerate quadrilateral detected").c_str());
00191 idx = NULL;
00192 }
00193 else if (!IsPlanar(p0, p1, p2, p3)) {
00194 luxError(LUX_CONSISTENCY, LUX_ERROR, string("Non-planar quadrilateral detected").c_str());
00195 idx = NULL;
00196 }
00197 else if (!IsConvex(p0, p1, p2, p3)) {
00198 luxError(LUX_CONSISTENCY, LUX_ERROR, string("Non-convex quadrilateral detected").c_str());
00199 idx = NULL;
00200 }
00201
00202 if (!idx)
00203 return;
00204
00205
00206 for(int i = 0; i < 4; i++) {
00207
00208 const Point &p00 = mesh->p[idx[0]];
00209 const Point &p10 = mesh->p[idx[1]];
00210 const Point &p11 = mesh->p[idx[2]];
00211 const Point &p01 = mesh->p[idx[3]];
00212
00213
00214 const Vector e01 = p10 - p00;
00215 const Vector e03 = p01 - p00;
00216 const Vector e02 = p11 - p00;
00217 const Vector N = Cross(e01, e03);
00218
00219 float a11 = 0.0f;
00220 float b11 = 0.0f;
00221
00222 ComputeV11BarycentricCoords(e01, e02, e03, &a11, &b11);
00223
00224 if ((a11 > 1.0f) || (b11 > 1.0f)) {
00225
00226
00227
00228
00229
00230
00231 int* nonConstIdx = const_cast<int*>(idx);
00232 const int tmp = nonConstIdx[0];
00233 nonConstIdx[0] = nonConstIdx[1];
00234 nonConstIdx[1] = nonConstIdx[2];
00235 nonConstIdx[2] = nonConstIdx[3];
00236 nonConstIdx[3] = tmp;
00237 } else
00238 break;
00239 }
00240 }
00241
00242 BBox MeshQuadrilateral::ObjectBound() const {
00243
00244 if (!idx)
00245 return BBox();
00246
00247
00248 const Point &p0 = mesh->p[idx[0]];
00249 const Point &p1 = mesh->p[idx[1]];
00250 const Point &p2 = mesh->p[idx[2]];
00251 const Point &p3 = mesh->p[idx[3]];
00252
00253 return Union(BBox(mesh->WorldToObject(p0), mesh->WorldToObject(p1)),
00254 BBox(mesh->WorldToObject(p2), mesh->WorldToObject(p3)));
00255 }
00256
00257 BBox MeshQuadrilateral::WorldBound() const {
00258
00259 if (!idx)
00260 return BBox();
00261
00262
00263 const Point &p0 = mesh->p[idx[0]];
00264 const Point &p1 = mesh->p[idx[1]];
00265 const Point &p2 = mesh->p[idx[2]];
00266 const Point &p3 = mesh->p[idx[3]];
00267
00268 return Union(BBox(p0, p1), BBox(p2, p3));
00269 }
00270
00271 bool MeshQuadrilateral::Intersect(const Ray &ray, Intersection *isect) const {
00272
00273
00274
00275
00276
00277
00278 if (!idx)
00279 return false;
00280
00281
00282 const Point &p00 = mesh->p[idx[0]];
00283 const Point &p10 = mesh->p[idx[1]];
00284 const Point &p11 = mesh->p[idx[2]];
00285 const Point &p01 = mesh->p[idx[3]];
00286
00287
00288
00289 const Vector e01 = p10 - p00;
00290 const Vector e03 = p01 - p00;
00291 const Vector P = Cross(ray.d, e03);
00292 float det = Dot(e01, P);
00293 if (fabsf(det) < 1e-7)
00294 return false;
00295
00296 float invdet = 1.f / det;
00297
00298 const Vector T = ray.o - p00;
00299 float alpha = Dot(T, P) * invdet;
00300 if (alpha < 0)
00301 return false;
00302
00303 const Vector Q = Cross(T, e01);
00304 float beta = Dot(ray.d, Q) * invdet;
00305 if (beta < 0)
00306 return false;
00307
00308
00309
00310 if ((alpha + beta) > 1.f) {
00311 const Vector e23 = p01 - p11;
00312 const Vector e21 = p10 - p11;
00313 const Vector P2 = Cross(ray.d, e21);
00314 float det2 = Dot(e23, P2);
00315 if (fabsf(det2) < 1e-7f)
00316 return false;
00317
00318
00319
00320 float invdet2 = (det2 < 0) ? -1 : 1;
00321 const Vector T2 = ray.o - p11;
00322 float alpha2 = Dot(T2, P2) * invdet2;
00323 if (alpha2 < 0)
00324 return false;
00325 const Vector Q2 = Cross(T2, e23);
00326 float beta2 = Dot(ray.d, Q2) * invdet2;
00327 if (beta2 < 0)
00328 return false;
00329 }
00330
00331
00332
00333 float t = Dot(e03, Q) * invdet;
00334 if (t < ray.mint || t > ray.maxt)
00335 return false;
00336
00337
00338 const Vector e02 = p11 - p00;
00339
00340 float a11 = 0.f;
00341 float b11 = 0.f;
00342
00343 ComputeV11BarycentricCoords(e01, e02, e03, &a11, &b11);
00344
00345
00346 a11 = a11 - 1;
00347 b11 = b11 - 1;
00348
00349
00350
00351 float u = 0.0f, v = 0.0f;
00352 if (fabsf(a11) < 1e-7f) {
00353 u = alpha;
00354 v = fabsf(b11) < 1e-7f ? beta : beta / (u * b11 + 1.f);
00355 } else if (fabsf(b11) < 1e-7f) {
00356 v = beta;
00357 u = alpha / (v * a11 + 1.f);
00358 } else {
00359 float A = -b11;
00360 float B = alpha*b11 - beta*a11 - 1.f;
00361 float C = alpha;
00362
00363 Quadratic(A, B, C, &u, &v);
00364 if ((u < 0) || (u > 1))
00365 u = v;
00366
00367 v = beta / (u * b11 + 1.f);
00368 }
00369
00370
00371
00372 Vector dpdu, dpdv;
00373 float uv[4][2];
00374 float A[3][3], InvA[3][3];
00375
00376 GetUVs(uv);
00377
00378 A[0][0] = uv[1][0] - uv[0][0];
00379 A[0][1] = uv[1][1] - uv[0][1];
00380 A[0][2] = uv[1][0]*uv[1][1] - uv[0][0]*uv[0][1];
00381 A[1][0] = uv[2][0] - uv[0][0];
00382 A[1][1] = uv[2][1] - uv[0][1];
00383 A[1][2] = uv[2][0]*uv[2][1] - uv[0][0]*uv[0][1];
00384 A[2][0] = uv[3][0] - uv[0][0];
00385 A[2][1] = uv[3][1] - uv[0][1];
00386 A[2][2] = uv[3][0]*uv[3][1] - uv[0][0]*uv[0][1];
00387
00388
00389 if (!Invert3x3(A, InvA)) {
00390
00391 Vector N = Cross(e01, e02);
00392 CoordinateSystem(Normalize(N), &dpdu, &dpdv);
00393 } else {
00394 dpdu = Vector(
00395 InvA[0][0] * e01.x + InvA[0][1] * e02.x + InvA[0][2] * e03.x,
00396 InvA[0][0] * e01.y + InvA[0][1] * e02.y + InvA[0][2] * e03.y,
00397 InvA[0][0] * e01.z + InvA[0][1] * e02.z + InvA[0][2] * e03.z);
00398 dpdv = Vector(
00399 InvA[1][0] * e01.x + InvA[1][1] * e02.x + InvA[1][2] * e03.x,
00400 InvA[1][0] * e01.y + InvA[1][1] * e02.y + InvA[1][2] * e03.y,
00401 InvA[1][0] * e01.z + InvA[1][1] * e02.z + InvA[1][2] * e03.z);
00402 }
00403
00404
00405 Normal nn;
00406 if (mesh->n)
00407 nn = Normalize(mesh->ObjectToWorld(
00408 ((1.0f - u) * (1.0f - v)) * mesh->n[idx[0]] +
00409 (u * (1.0f - v)) * mesh->n[idx[1]] +
00410 (u * v) * mesh->n[idx[2]] +
00411 ((1.0f - u) * v) * mesh->n[idx[3]]));
00412 else {
00413 Vector N = Cross(e01, e02);
00414 nn = Normal(Normalize(N));
00415 }
00416
00417 if(isect) {
00418 isect->dg = DifferentialGeometry(ray(t),
00419 nn,
00420 dpdu, dpdv,
00421 Normal(0, 0, 0), Normal(0, 0, 0),
00422 u, v, this);
00423 isect->dg.AdjustNormal(mesh->reverseOrientation, mesh->transformSwapsHandedness);
00424 isect->Set(mesh->WorldToObject, this, mesh->GetMaterial().get());
00425 ray.maxt = t;
00426 }
00427
00428 return true;
00429 }
00430
00431 bool MeshQuadrilateral::IntersectP(const Ray &ray) const {
00432 return Intersect(ray, NULL);
00433 }
00434
00435 float MeshQuadrilateral::Area() const {
00436
00437 if (!idx)
00438 return 0.f;
00439
00440
00441 const Point &p0 = mesh->p[idx[0]];
00442 const Point &p1 = mesh->p[idx[1]];
00443 const Point &p3 = mesh->p[idx[3]];
00444
00445 Vector P = p1 - p0;
00446 Vector Q = p3 - p1;
00447
00448 Vector PxQ = Cross(P, Q);
00449
00450 return 0.5f * PxQ.Length();
00451 }
00452
00453 void MeshQuadrilateral::GetShadingGeometry(const Transform &obj2world,
00454 const DifferentialGeometry &dg,
00455 DifferentialGeometry *dgShading) const
00456 {
00457 if (!mesh->n) {
00458 *dgShading = dg;
00459 if(!mesh->uvs) {
00460
00461 const BBox bounds = MeshQuadrilateral::WorldBound();
00462 int maxExtent = bounds.MaximumExtent();
00463 float maxSize = bounds.pMax[maxExtent] - bounds.pMin[maxExtent];
00464 dgShading->dpdu *= (maxSize * .1f) / dgShading->dpdu.Length();
00465 dgShading->dpdv *= (maxSize * .1f) / dgShading->dpdv.Length();
00466 }
00467 return;
00468 }
00469
00470
00471 Normal ns = dg.nn;
00472 float lenDpDu = dg.dpdu.Length();
00473 float lenDpDv = dg.dpdv.Length();
00474 Vector ts = Normalize(Cross(dg.dpdu, ns));
00475 Vector ss = Cross(ts, ns);
00476
00477 if(mesh->uvs) {
00478 ss *= lenDpDu;
00479 ts *= lenDpDv;
00480 }
00481 else {
00482 const BBox bounds = MeshQuadrilateral::WorldBound();
00483 int maxExtent = bounds.MaximumExtent();
00484 float maxSize = bounds.pMax[maxExtent] - bounds.pMin[maxExtent];
00485 ss *= maxSize * .1f;
00486 ts *= maxSize * .1f;
00487 }
00488
00489
00490
00491 Normal dndu, dndv;
00492 float uv[4][2];
00493 float A[3][3], InvA[3][3];
00494
00495 GetUVs(uv);
00496
00497 A[0][0] = uv[1][0] - uv[0][0];
00498 A[0][1] = uv[1][1] - uv[0][1];
00499 A[0][2] = uv[1][0]*uv[1][1] - uv[0][0]*uv[0][1];
00500 A[1][0] = uv[2][0] - uv[0][0];
00501 A[1][1] = uv[2][1] - uv[0][1];
00502 A[1][2] = uv[2][0]*uv[2][1] - uv[0][0]*uv[0][1];
00503 A[2][0] = uv[3][0] - uv[0][0];
00504 A[2][1] = uv[3][1] - uv[0][1];
00505 A[2][2] = uv[3][0]*uv[3][1] - uv[0][0]*uv[0][1];
00506
00507
00508 if (!Invert3x3(A, InvA)) {
00509
00510 dndu = dndv = Normal(0, 0, 0);
00511 } else {
00512 const Normal &n00 = mesh->n[idx[0]];
00513 const Normal &n10 = mesh->n[idx[1]];
00514 const Normal &n11 = mesh->n[idx[2]];
00515 const Normal &n01 = mesh->n[idx[3]];
00516
00517 const Normal dn01 = n10 - n00;
00518 const Normal dn02 = n11 - n00;
00519 const Normal dn03 = n01 - n00;
00520
00521 dndu = Normal(
00522 InvA[0][0] * dn01.x + InvA[0][1] * dn02.x + InvA[0][2] * dn03.x,
00523 InvA[0][0] * dn01.y + InvA[0][1] * dn02.y + InvA[0][2] * dn03.y,
00524 InvA[0][0] * dn01.z + InvA[0][1] * dn02.z + InvA[0][2] * dn03.z);
00525 dndv = Normal(
00526 InvA[1][0] * dn01.x + InvA[1][1] * dn02.x + InvA[1][2] * dn03.x,
00527 InvA[1][0] * dn01.y + InvA[1][1] * dn02.y + InvA[1][2] * dn03.y,
00528 InvA[1][0] * dn01.z + InvA[1][1] * dn02.z + InvA[1][2] * dn03.z);
00529
00530 dndu = obj2world(dndu);
00531 dndv = obj2world(dndv);
00532 }
00533
00534 *dgShading = DifferentialGeometry(
00535 dg.p,
00536 ns,
00537 ss, ts,
00538 dndu, dndv,
00539 dg.u, dg.v, this);
00540
00541 dgShading->dudx = dg.dudx; dgShading->dvdx = dg.dvdx;
00542 dgShading->dudy = dg.dudy; dgShading->dvdy = dg.dvdy;
00543 dgShading->dpdx = dg.dpdx; dgShading->dpdy = dg.dpdy;
00544 }