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 "quadrilateral.h"
00025 #include "geometry/matrix3x3.h"
00026
00027 using namespace lux;
00028
00029 int MajorAxis(const Vector &v) {
00030 float absVx = fabsf(v.x);
00031 float absVy = fabsf(v.y);
00032 float absVz = fabsf(v.z);
00033
00034 if (absVx > absVy)
00035 return (absVx > absVz) ? 0 : 2;
00036 return (absVy > absVz) ? 1 : 2;
00037 }
00038
00039 BBox QuadMesh::ObjectBound() const {
00040 luxError(LUX_BUG,LUX_SEVERE, "Unimplemented QuadMesh::ObjectBound() method called");
00041 return BBox();
00042 }
00043
00044
00045
00046 Quadrilateral::Quadrilateral(const lux::Transform &o2w, bool ro,
00047 QuadMesh *m, int *indices) : Shape(o2w, ro) {
00048
00049 idx = indices;
00050 mesh = m;
00051 }
00052
00053 Quadrilateral::~Quadrilateral() {
00054 }
00055
00056 BBox Quadrilateral::ObjectBound() const {
00057
00058 const Point &p0 = mesh->p[idx[0]];
00059 const Point &p1 = mesh->p[idx[1]];
00060 const Point &p2 = mesh->p[idx[2]];
00061 const Point &p3 = mesh->p[idx[3]];
00062 return Union(BBox(WorldToObject(p0), WorldToObject(p1)),
00063 BBox(WorldToObject(p2), WorldToObject(p3)));
00064 }
00065
00066 BBox Quadrilateral::WorldBound() const {
00067
00068 const Point &p0 = mesh->p[idx[0]];
00069 const Point &p1 = mesh->p[idx[1]];
00070 const Point &p2 = mesh->p[idx[2]];
00071 const Point &p3 = mesh->p[idx[3]];
00072 return Union(BBox(p0, p1), BBox(p2, p3));
00073 }
00074
00075 bool Quadrilateral::Intersect(const Ray &ray, float *tHit,
00076 DifferentialGeometry *dg) const {
00077
00078
00079
00080
00081
00082
00083 const Point &p00 = mesh->p[idx[0]];
00084 const Point &p10 = mesh->p[idx[1]];
00085 const Point &p11 = mesh->p[idx[2]];
00086 const Point &p01 = mesh->p[idx[3]];
00087
00088
00089
00090 Vector e01 = p10 - p00;
00091 Vector e03 = p01 - p00;
00092 Vector P = Cross(ray.d, e03);
00093 float det = Dot(e01, P);
00094 if (fabsf(det) < 1e-7)
00095 return false;
00096
00097 float invdet = 1.f / det;
00098
00099 Vector T = ray.o - p00;
00100 float alpha = Dot(T, P) * invdet;
00101 if (alpha < 0 || alpha > 1)
00102 return false;
00103
00104 Vector Q = Cross(T, e01);
00105 float beta = Dot(ray.d, Q) * invdet;
00106 if (beta < 0 || beta > 1)
00107 return false;
00108
00109
00110
00111 if ((alpha + beta) > 1.f) {
00112 Vector e23 = p01 - p11;
00113 Vector e21 = p10 - p11;
00114 Vector P2 = Cross(ray.d, e21);
00115 float det2 = Dot(e23, P2);
00116 if (fabsf(det2) < 1e-7f)
00117 return false;
00118
00119
00120
00121 float invdet2 = (det2 < 0) ? -1 : 1;
00122 Vector T2 = ray.o - p11;
00123 float alpha2 = Dot(T2, P2) * invdet2;
00124 if (alpha2 < 0)
00125 return false;
00126 Vector Q2 = Cross(T2, e23);
00127 float beta2 = Dot(ray.d, Q2) * invdet2;
00128 if (beta2 < 0)
00129 return false;
00130 }
00131
00132
00133
00134 float t = Dot(e03, Q) * invdet;
00135 if (t < ray.mint || t > ray.maxt)
00136 return false;
00137
00138
00139 Vector e02 = p11 - p00;
00140 Vector N = Cross(e01, e03);
00141
00142 float a11 = 0.0f, b11 = 0.0f;
00143 int Nma = MajorAxis(N);
00144
00145 switch (Nma) {
00146 case 0: {
00147 float iNx = 1.f / N.x;
00148 a11 = (e02.y * e03.z - e02.z * e03.y) * iNx - 1.f;
00149 b11 = (e01.y * e02.z - e01.z * e02.y) * iNx - 1.f;
00150 break;
00151 }
00152 case 1: {
00153 float iNy = 1.f / N.y;
00154 a11 = (e02.z * e03.x - e02.x * e03.z) * iNy - 1.f;
00155 b11 = (e01.z * e02.x - e01.x * e02.z) * iNy - 1.f;
00156 break;
00157 }
00158 case 2: {
00159 float iNz = 1.f / N.z;
00160 a11 = (e02.x * e03.y - e02.y * e03.x) * iNz - 1.f;
00161 b11 = (e01.x * e02.y - e01.y * e02.x) * iNz - 1.f;
00162 break;
00163 }
00164 }
00165
00166
00167
00168 float u = 0.0f, v = 0.0f;
00169 if (fabsf(a11) < 1e-7f) {
00170 u = alpha;
00171 v = fabsf(b11) < 1e-7f ? beta : beta / (u * b11 + 1.f);
00172 } else if (fabsf(b11) < 1e-7f) {
00173 v = beta;
00174 u = alpha / (v * a11 + 1.f);
00175 } else {
00176 float A = -b11;
00177 float B = alpha*b11 - beta*a11 - 1.f;
00178 float C = alpha;
00179
00180 Quadratic(A, B, C, &u, &v);
00181 if ((u < 0) || (u > 1))
00182 u = v;
00183
00184 v = beta / (u * b11 + 1.f);
00185 }
00186
00187
00188
00189 Vector dpdu, dpdv;
00190 float uv[4][2];
00191 float A[3][3], InvA[3][3];
00192
00193 GetUVs(uv);
00194
00195 A[0][0] = uv[1][0] - uv[0][0];
00196 A[0][1] = uv[1][1] - uv[0][1];
00197 A[0][2] = uv[1][0]*uv[1][1] - uv[0][0]*uv[0][1];
00198 A[1][0] = uv[2][0] - uv[0][0];
00199 A[1][1] = uv[2][1] - uv[0][1];
00200 A[1][2] = uv[2][0]*uv[2][1] - uv[0][0]*uv[0][1];
00201 A[2][0] = uv[3][0] - uv[0][0];
00202 A[2][1] = uv[3][1] - uv[0][1];
00203 A[2][2] = uv[3][0]*uv[3][1] - uv[0][0]*uv[0][1];
00204
00205
00206 if (!Invert3x3(A, InvA)) {
00207
00208 Vector N = Cross(e01, e02);
00209 CoordinateSystem(Normalize(N), &dpdu, &dpdv);
00210 } else {
00211 dpdu = Vector(
00212 InvA[0][0] * e01.x + InvA[0][1] * e02.x + InvA[0][2] * e03.x,
00213 InvA[0][0] * e01.y + InvA[0][1] * e02.y + InvA[0][2] * e03.y,
00214 InvA[0][0] * e01.z + InvA[0][1] * e02.z + InvA[0][2] * e03.z);
00215 dpdv = Vector(
00216 InvA[1][0] * e01.x + InvA[1][1] * e02.x + InvA[1][2] * e03.x,
00217 InvA[1][0] * e01.y + InvA[1][1] * e02.y + InvA[1][2] * e03.y,
00218 InvA[1][0] * e01.z + InvA[1][1] * e02.z + InvA[1][2] * e03.z);
00219 }
00220
00221 *dg = DifferentialGeometry(ray(t), dpdu, dpdv,
00222 Normal(0, 0, 0), Normal(0, 0, 0),
00223 u, v, this);
00224
00225 *tHit = t;
00226
00227 return true;
00228 }
00229
00230 bool Quadrilateral::IntersectP(const Ray &ray) const {
00231
00232 float tHit;
00233 DifferentialGeometry dg;
00234
00235 return Intersect(ray, &tHit, &dg);
00236 }
00237
00238 float Quadrilateral::Area() const {
00239
00240
00241 const Point &p0 = mesh->p[idx[0]];
00242 const Point &p1 = mesh->p[idx[1]];
00243 const Point &p3 = mesh->p[idx[3]];
00244
00245 Vector P = p1 - p0;
00246 Vector Q = p3 - p1;
00247
00248 Vector PxQ = Cross(P, Q);
00249
00250 return 0.5f * PxQ.Length();
00251 }