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 "mesh.h"
00025 #include "mc.h"
00026
00027 using namespace lux;
00028
00029 MeshWaldTriangle::MeshWaldTriangle(const Mesh *m, int n)
00030 : MeshBaryTriangle(m, n)
00031 {
00032
00033 const float l0 = DistanceSquared(mesh->p[v[0]], mesh->p[v[1]]);
00034 const float l1 = DistanceSquared(mesh->p[v[1]], mesh->p[v[2]]);
00035 const float l2 = DistanceSquared(mesh->p[v[2]], mesh->p[v[0]]);
00036 const float d0 = fabsf(l0 - l2);
00037 const float d1 = fabsf(l1 - l0);
00038 const float d2 = fabsf(l2 - l1);
00039 int *v_ = const_cast<int *>(v);
00040 if (d2 < d1 && d2 < d0) {
00041 swap(v_[0], v_[2]);
00042 swap(v_[1], v_[2]);
00043 } else if (d1 < d0) {
00044 swap(v_[0], v_[1]);
00045 swap(v_[2], v_[1]);
00046 }
00047
00048
00049
00050
00051 const Point &v0 = mesh->p[v[0]];
00052 const Point &v1 = mesh->p[v[1]];
00053 const Point &v2 = mesh->p[v[2]];
00054 Vector e1 = v1 - v0;
00055 Vector e2 = v2 - v0;
00056
00057 normalizedNormal = Normal(Normalize(Cross(e1, e2)));
00058
00059 if (isnan(normalizedNormal.x) || isnan(normalizedNormal.y) ||
00060 isnan(normalizedNormal.z)) {
00061 intersectionType = DEGENERATE;
00062 return;
00063 }
00064
00065
00066
00067
00068 if ((fabs(normalizedNormal.x) > fabs(normalizedNormal.y)) &&
00069 (fabs(normalizedNormal.x) > fabs(normalizedNormal.z)))
00070 intersectionType = DOMINANT_X;
00071 else if (fabs(normalizedNormal.y) > fabs(normalizedNormal.z))
00072 intersectionType = DOMINANT_Y;
00073 else
00074 intersectionType = DOMINANT_Z;
00075
00076 float ax, ay, bx, by, cx, cy;
00077 switch (intersectionType) {
00078 case DOMINANT_X: {
00079 const float invNormal = 1.f / normalizedNormal.x;
00080 nu = normalizedNormal.y * invNormal;
00081 nv = normalizedNormal.z * invNormal;
00082 nd = v0.x + nu * v0.y + nv * v0.z;
00083 ax = v0.y;
00084 ay = v0.z;
00085 bx = v2.y - ax;
00086 by = v2.z - ay;
00087 cx = v1.y - ax;
00088 cy = v1.z - ay;
00089 break;
00090 }
00091 case DOMINANT_Y: {
00092 const float invNormal = 1.f / normalizedNormal.y;
00093 nu = normalizedNormal.z * invNormal;
00094 nv = normalizedNormal.x * invNormal;
00095 nd = nv * v0.x + v0.y + nu * v0.z;
00096 ax = v0.z;
00097 ay = v0.x;
00098 bx = v2.z - ax;
00099 by = v2.x - ay;
00100 cx = v1.z - ax;
00101 cy = v1.x - ay;
00102 break;
00103 }
00104 case DOMINANT_Z: {
00105 const float invNormal = 1.f / normalizedNormal.z;
00106 nu = normalizedNormal.x * invNormal;
00107 nv = normalizedNormal.y * invNormal;
00108 nd = nu * v0.x + nv * v0.y + v0.z;
00109 ax = v0.x;
00110 ay = v0.y;
00111 bx = v2.x - ax;
00112 by = v2.y - ay;
00113 cx = v1.x - ax;
00114 cy = v1.y - ay;
00115 break;
00116 }
00117 default:
00118 BOOST_ASSERT(false);
00119
00120 return;
00121 }
00122
00123 float det = bx * cy - by * cx;
00124 float invDet = 1.f / det;
00125
00126 bnu = -by * invDet;
00127 bnv = bx * invDet;
00128 bnd = (by * ax - bx * ay) * invDet;
00129 cnu = cy * invDet;
00130 cnv = -cx * invDet;
00131 cnd = (cx * ay - cy * ax) * invDet;
00132
00133
00134
00135
00136
00137 float uvs[3][2];
00138 GetUVs(uvs);
00139
00140 const float du1 = uvs[0][0] - uvs[2][0];
00141 const float du2 = uvs[1][0] - uvs[2][0];
00142 const float dv1 = uvs[0][1] - uvs[2][1];
00143 const float dv2 = uvs[1][1] - uvs[2][1];
00144 const Vector dp1 = v0 - v2, dp2 = v1 - v2;
00145 const float determinant = du1 * dv2 - dv1 * du2;
00146 if (determinant == 0.f) {
00147
00148 CoordinateSystem(Vector(normalizedNormal), &dpdu, &dpdv);
00149 } else {
00150 const float invdet = 1.f / determinant;
00151 dpdu = ( dv2 * dp1 - dv1 * dp2) * invdet;
00152 dpdv = (-du2 * dp1 + du1 * dp2) * invdet;
00153 }
00154 }
00155
00156 bool MeshWaldTriangle::Intersect(const Ray &ray, Intersection *isect) const
00157 {
00158 float o0, o1, o2, d0, d1, d2;
00159 switch (intersectionType) {
00160 case DOMINANT_X: {
00161 o0 = ray.o.x;
00162 o1 = ray.o.y;
00163 o2 = ray.o.z;
00164 d0 = ray.d.x;
00165 d1 = ray.d.y;
00166 d2 = ray.d.z;
00167 break;
00168 }
00169 case DOMINANT_Y: {
00170 o0 = ray.o.y;
00171 o1 = ray.o.z;
00172 o2 = ray.o.x;
00173 d0 = ray.d.y;
00174 d1 = ray.d.z;
00175 d2 = ray.d.x;
00176 break;
00177 }
00178 case DOMINANT_Z: {
00179 o0 = ray.o.z;
00180 o1 = ray.o.x;
00181 o2 = ray.o.y;
00182 d0 = ray.d.z;
00183 d1 = ray.d.x;
00184 d2 = ray.d.y;
00185 break;
00186 }
00187 default:
00188 return false;
00189 }
00190 const float det = d0 + nu * d1 + nv * d2;
00191 if (det == 0.f)
00192 return false;
00193
00194 const float t = nd - o0 - nu * o1 - nv * o2;
00195 if (det > 0.f) {
00196 if (t <= det * ray.mint || t >= det * ray.maxt)
00197 return false;
00198 } else {
00199 if (t >= det * ray.mint || t <= det * ray.maxt)
00200 return false;
00201 }
00202
00203 const float hu = det * o1 + t * d1;
00204 const float hv = det * o2 + t * d2;
00205 const float uu = (hu * bnu + hv * bnv) / det + bnd;
00206 if (uu < 0.f)
00207 return false;
00208
00209 const float vv = (hu * cnu + hv * cnv) / det + cnd;
00210 if (vv < 0.f)
00211 return false;
00212
00213 const float b0 = 1.f - uu - vv;
00214 if (b0 < 0.f)
00215 return false;
00216
00217 const float tt = t / det;
00218
00219 float uvs[3][2];
00220 GetUVs(uvs);
00221
00222 const float tu = b0 * uvs[0][0] + uu * uvs[1][0] + vv * uvs[2][0];
00223 const float tv = b0 * uvs[0][1] + uu * uvs[1][1] + vv * uvs[2][1];
00224
00225 const Point pp(b0 * mesh->p[v[0]] + uu * mesh->p[v[1]] + vv * mesh->p[v[2]]);
00226
00227 isect->dg = DifferentialGeometry(pp, normalizedNormal, dpdu, dpdv,
00228 Normal(0, 0, 0), Normal(0, 0, 0), tu, tv, this);
00229 isect->Set(mesh->WorldToObject, this, mesh->GetMaterial().get());
00230 isect->dg.triangleBaryCoords[0] = b0;
00231 isect->dg.triangleBaryCoords[1] = uu;
00232 isect->dg.triangleBaryCoords[2] = vv;
00233 ray.maxt = tt;
00234
00235 return true;
00236 }
00237
00238 bool MeshWaldTriangle::IntersectP(const Ray &ray) const
00239 {
00240 float o0, o1, o2, d0, d1, d2;
00241 switch (intersectionType) {
00242 case DOMINANT_X: {
00243 o0 = ray.o.x;
00244 o1 = ray.o.y;
00245 o2 = ray.o.z;
00246 d0 = ray.d.x;
00247 d1 = ray.d.y;
00248 d2 = ray.d.z;
00249 break;
00250 }
00251 case DOMINANT_Y: {
00252 o0 = ray.o.y;
00253 o1 = ray.o.z;
00254 o2 = ray.o.x;
00255 d0 = ray.d.y;
00256 d1 = ray.d.z;
00257 d2 = ray.d.x;
00258 break;
00259 }
00260 case DOMINANT_Z: {
00261 o0 = ray.o.z;
00262 o1 = ray.o.x;
00263 o2 = ray.o.y;
00264 d0 = ray.d.z;
00265 d1 = ray.d.x;
00266 d2 = ray.d.y;
00267 break;
00268 }
00269 default:
00270 return false;
00271 }
00272 const float det = d0 + nu * d1 + nv * d2;
00273 if (det == 0.f)
00274 return false;
00275
00276 const float t = nd - o0 - nu * o1 - nv * o2;
00277 if (det > 0.f) {
00278 if (t <= det * ray.mint || t >= det * ray.maxt)
00279 return false;
00280 } else {
00281 if (t >= det * ray.mint || t <= det * ray.maxt)
00282 return false;
00283 }
00284
00285 const float hu = det * o1 + t * d1;
00286 const float hv = det * o2 + t * d2;
00287 const float uu = (hu * bnu + hv * bnv) / det + bnd;
00288 if (uu < 0.f)
00289 return false;
00290
00291 const float vv = (hu * cnu + hv * cnv) / det + cnd;
00292 if (vv < 0.f || uu + vv > 1.f)
00293 return false;
00294
00295 return true;
00296 }
00297
00298 void MeshWaldTriangle::Sample(float u1, float u2, float u3, DifferentialGeometry *dg) const {
00299 float b1, b2;
00300 UniformSampleTriangle(u1, u2, &b1, &b2);
00301
00302 const Point &p1 = mesh->p[v[0]];
00303 const Point &p2 = mesh->p[v[1]];
00304 const Point &p3 = mesh->p[v[2]];
00305 float b3 = 1.f - b1 - b2;
00306 dg->p = b1 * p1 + b2 * p2 + b3 * p3;
00307
00308 dg->nn = normalizedNormal;
00309 dg->dpdu = dpdu;
00310 dg->dpdv = dpdv;
00311
00312 float uv[3][2];
00313 GetUVs(uv);
00314 dg->u = b1 * uv[0][0] + b2 * uv[1][0] + b3 * uv[2][0];
00315 dg->v = b1 * uv[0][1] + b2 * uv[1][1] + b3 * uv[2][1];
00316 }
00317
00318 bool MeshWaldTriangle::isDegenerate() const {
00319 return intersectionType == DEGENERATE;
00320 }