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 "torus.h"
00025 #include "paramset.h"
00026 #include "dynload.h"
00027
00028 using namespace lux;
00029
00030 int signR(double Z)
00031 {
00032 if (Z > 0.0) return 1;
00033 if (Z < 0.0) return -1;
00034
00035 return 0;
00036 }
00037
00038 double CBRT(double Z)
00039 {
00040 double ret;
00041 const double THIRD = 1./3.;
00042
00043
00044
00045
00046
00047
00048 ret = pow(fabs(Z),THIRD) * signR(Z);
00049 return ret;
00050 }
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072 int cubic(double A[4], double X[3])
00073 {
00074 const double PI = 3.1415926535897932;
00075 const double THIRD = 1./3.;
00076 double W, P, Q, DIS, PHI;
00077 int L;
00078
00079
00080
00081
00082
00083
00084 if (A[3] != 0.0) {
00085
00086 W = A[2]/A[3]*THIRD;
00087 double W2 = W*W;
00088 double P1 = -(A[1]/A[3]*THIRD - W2);
00089 P = P1*P1*P1;
00090 Q = -.5*(2.0*W2*W-(A[1]*W-A[0])/A[3]);
00091 DIS = Q*Q-P;
00092 if ( DIS < 0.0 ) {
00093
00094
00095 PHI = acos(min(1.0,max(-1.0,Q/sqrt(P))));
00096 P = 2.0*sqrt(P1);
00097 for (int i = 0; i < 3; i++)
00098 X[i] = P*cos((PHI+2*((double)i)*PI)*THIRD)-W;
00099 L = 3;
00100 }
00101 else {
00102
00103 DIS = sqrt(DIS);
00104 X[0] = CBRT(Q+DIS)+CBRT(Q-DIS)-W;
00105 L = 1;
00106 }
00107 }
00108 else if (A[2] != 0.0) {
00109
00110 P = 0.5*A[1]/A[2];
00111 DIS = P*P-A[0]/A[2];
00112 if (DIS > 0.0) {
00113
00114 DIS = sqrt(DIS);
00115 X[0] = -P - DIS;
00116 X[1] = -P + DIS;
00117 L = 2;
00118 }
00119 else {
00120
00121 return 0;
00122 }
00123 }
00124 else if (A[1] != 0.0) {
00125
00126 X[0] = A[0]/A[1];
00127 L = 1;
00128 }
00129 else {
00130
00131 return 0;
00132 }
00133
00134
00135
00136 for (int i = 0; i < L; i++) {
00137 X[i] -= (A[0]+X[i]*(A[1]+X[i]*(A[2]+X[i]*A[3]))) / (A[1]+X[i]*(2.0*A[2]+X[i]*3.0*A[3]));
00138 }
00139
00140 return L;
00141 }
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163 int quartic(double dd[5], double sol[4])
00164 {
00165 double AA[4], z[3];
00166 double a, b, c, d, f, p, q, r, zsol, xK2, xL, xK, sqp, sqm;
00167 int ncube;
00168 int Nsol = 0;
00169
00170 if (dd[4] == 0.0)
00171 {
00172 return 0;
00173 }
00174
00175 a = dd[4];
00176 b = dd[3];
00177 c = dd[2];
00178 d = dd[1];
00179 f = dd[0];
00180
00181 const double aa = a*a;
00182 const double bb = b*b;
00183
00184 p = (-3.0*b*b + 8.0 *a*c)/(8.0*aa);
00185 q = (bb*b - 4.0*a*b*c + 8.0*d*aa) / (8.0*aa*a);
00186 r = (-3.0*bb*bb + 16.0*a*bb*c - 64.0*aa*b*d + 256.0*aa*a*f)/(256.0*aa*aa);
00187
00188
00189 AA[3] = 8.0;
00190 AA[2] = -4.0*p;
00191 AA[1] = -8.0*r;
00192 AA[0] = 4.0*p*r - q*q;
00193
00194 ncube = cubic(AA, z);
00195
00196 zsol = -1.e99;
00197 for (int i = 0; i < ncube; i++)
00198 zsol = max(zsol, z[i]);
00199
00200 xK2 = 2.0*zsol - p;
00201 xK = sqrt(xK2);
00202 xL = q/(2.0*xK);
00203 sqp = xK2 - 4.0 * (zsol + xL);
00204 sqm = xK2 - 4.0 * (zsol - xL);
00205
00206 if (sqp >= 0.0) {
00207 sqp = sqrt(sqp);
00208
00209 if (sqm >= 0.0) {
00210 sqm = sqrt(sqm);
00211 sol[3] = 0.5 * (xK + sqp);
00212 sol[2] = 0.5 * (xK - sqp);
00213 sol[1] = 0.5 * (-xK + sqm);
00214 sol[0] = 0.5 * (-xK - sqm);
00215 Nsol = 4;
00216 }
00217 else {
00218 sol[1] = 0.5 * (xK + sqp);
00219 sol[0] = 0.5 * (xK - sqp);
00220 Nsol = 2;
00221 }
00222 }
00223 else {
00224 if (sqm < 0.0)
00225
00226 return 0;
00227
00228 sqm = sqrt(sqm);
00229
00230
00231 sol[1] = 0.5 * (-xK + sqm);
00232 sol[0] = 0.5 * (-xK - sqm);
00233 Nsol = 2;
00234 }
00235
00236 for (int i = 0; i < Nsol; i++)
00237 sol[i] -= b/(4.0*a);
00238
00239 return Nsol;
00240 }
00241
00242 bool Torus::FindIntersection(const Ray &ray, float *tHit, Point *pHit,
00243 float *phiHit, float *thetaHit) const {
00244
00245
00246
00247 const double r2 = minorRadius * minorRadius;
00248 const double R2 = majorRadius * majorRadius;
00249
00250 const double dd = Dot(ray.d, ray.d);
00251 const double pd = Dot(Vector(ray.o), ray.d);
00252 const double pp = Dot(Vector(ray.o), Vector(ray.o));
00253 const double prR = pp - r2 - R2;
00254
00255 double coef[5];
00256
00257 coef[4] = dd*dd;
00258 coef[3] = 4.0*dd*pd;
00259 coef[2] = 4.0*pd*pd + 2.0*dd*prR + 4.0*R2*ray.d.z*ray.d.z;
00260 coef[1] = 4.0*pd*prR + 8.0*R2*ray.o.z*ray.d.z;
00261 coef[0] = prR*prR + 4.0*R2*(ray.o.z*ray.o.z - r2);
00262
00263 double t[4];
00264 int Nsol;
00265
00266
00267
00268 Nsol = quartic(coef, t);
00269 if (Nsol < 1)
00270 return false;
00271
00272
00273 double tmax = t[Nsol-1];
00274 if (tmax < ray.mint)
00275 return false;
00276
00277
00278 double thit = t[0];
00279 int ti = 0;
00280 while (thit < ray.mint) {
00281 ti++;
00282 if (ti >= Nsol)
00283 return false;
00284 thit = t[ti];
00285 }
00286 if (thit > ray.maxt)
00287 return false;
00288
00289
00290 Point phit;
00291 float phi, theta;
00292 while (true) {
00293 phit = ray(thit);
00294 phi = atan2f(phit.y, phit.x);
00295 if (phi < 0.f) phi += 2.f*M_PI;
00296
00297
00298 float sintheta = Clamp(phit.z / minorRadius, -1.f, 1.f);
00299 theta = asinf(sintheta);
00300
00301
00302 if (phit.x*phit.x + phit.y*phit.y < R2)
00303 theta = M_PI - theta;
00304 if (theta < 0.f) theta += 2.f*M_PI;
00305
00306
00307 if (!(theta < thetaMin || theta > thetaMax || phi > phiMax))
00308 break;
00309
00310
00311 ti++;
00312 if (ti >= Nsol)
00313 return false;
00314 thit = t[ti];
00315 if (thit > ray.maxt)
00316 return false;
00317 }
00318
00319
00320
00321 const float costheta = cosf(theta);
00322 pHit->x = cosf(phi)*(majorRadius + minorRadius*costheta);
00323 pHit->y = sinf(phi)*(majorRadius + minorRadius*costheta);
00324 pHit->z = minorRadius*sinf(theta);
00325 *tHit = (*pHit - ray.o).Length() / ray.d.Length();
00326 *phiHit = phi;
00327 *thetaHit = theta;
00328
00329 return true;
00330 }
00331
00332
00333 Torus::Torus(const Transform &o2w, bool ro, float marad, float mirad,
00334 float tmi, float tma, float pm)
00335 : Shape(o2w, ro) {
00336 majorRadius = marad;
00337 minorRadius = mirad;
00338 thetaMin = Radians(Clamp(min(tmi, tma), 0.f, 360.f));
00339 thetaMax = Radians(Clamp(max(tmi, tma), 0.f, 360.f));
00340 phiMax = Radians(Clamp(pm, 0.0f, 360.0f));
00341
00342 if (thetaMin < M_PI && thetaMax > M_PI)
00343 zmin = -minorRadius;
00344 else
00345 zmin = minorRadius * min(cosf(thetaMin), cosf(thetaMax));
00346 zmax = minorRadius * max(cosf(thetaMin), cosf(thetaMax));
00347 }
00348
00349 BBox Torus::ObjectBound() const {
00350 return BBox(Point(-majorRadius-minorRadius, -majorRadius-minorRadius, -minorRadius),
00351 Point(majorRadius+minorRadius, majorRadius+minorRadius, minorRadius));
00352 }
00353
00354 bool Torus::Intersect(const Ray &r, float *tHit,
00355 DifferentialGeometry *dg) const {
00356 float phi, theta;
00357 float thit;
00358 Point phit;
00359
00360 Ray ray;
00361 WorldToObject(r, &ray);
00362
00363 if (!FindIntersection(ray, &thit, &phit, &phi, &theta))
00364 return false;
00365
00366
00367
00368
00369
00370
00371
00372 const float u = phi / phiMax;
00373 const float v = (theta - thetaMin) / (thetaMax - thetaMin);
00374
00375
00376 float cosphi, sinphi;
00377 Vector dpdu, dpdv;
00378
00379 const float costheta = cosf(theta);
00380
00381 float zradius = sqrtf(phit.x*phit.x + phit.y*phit.y);
00382 if (zradius == 0)
00383 {
00384
00385
00386 cosphi = 0;
00387 sinphi = 1;
00388 dpdv = (thetaMax-thetaMin) *
00389 Vector(-phit.z * cosphi, -phit.z * sinphi,
00390 minorRadius * costheta);
00391 Vector norm = Vector(phit);
00392 dpdu = Cross(dpdv, norm);
00393 }
00394 else
00395 {
00396 float invzradius = 1.f / zradius;
00397 cosphi = phit.x * invzradius;
00398 sinphi = phit.y * invzradius;
00399 dpdu = Vector(-phiMax * phit.y, phiMax * phit.x, 0);
00400 dpdv = (thetaMax-thetaMin) *
00401 Vector(-phit.z * cosphi, -phit.z * sinphi,
00402 minorRadius * costheta);
00403 }
00404
00405 Vector d2Pduu = -phiMax * phiMax * Vector(phit.x, phit.y, 0);
00406 Vector d2Pduv = (thetaMax - thetaMin) *
00407 phit.z * phiMax *
00408 Vector(sinphi, -cosphi, 0.);
00409 Vector d2Pdvv = -(thetaMax - thetaMin) *
00410 (thetaMax - thetaMin) *
00411 Vector(minorRadius*cosphi*costheta, minorRadius*sinphi*costheta, phit.z);
00412
00413 float E = Dot(dpdu, dpdu);
00414 float F = Dot(dpdu, dpdv);
00415 float G = Dot(dpdv, dpdv);
00416 Vector N = Normalize(Cross(dpdu, dpdv));
00417 float e = Dot(N, d2Pduu);
00418 float f = Dot(N, d2Pduv);
00419 float g = Dot(N, d2Pdvv);
00420
00421 float invEGF2 = 1.f / (E*G - F*F);
00422 Normal dndu((f*F - e*G) * invEGF2 * dpdu +
00423 (e*F - f*E) * invEGF2 * dpdv);
00424 Normal dndv((g*F - f*G) * invEGF2 * dpdu +
00425 (f*F - g*E) * invEGF2 * dpdv);
00426
00427 *dg = DifferentialGeometry(ObjectToWorld(phit),
00428 ObjectToWorld(dpdu),
00429 ObjectToWorld(dpdv),
00430 ObjectToWorld(dndu),
00431 ObjectToWorld(dndv),
00432 u, v, this);
00433
00434 *tHit = thit;
00435 return true;
00436 }
00437
00438 bool Torus::IntersectP(const Ray &r) const {
00439 float phi, theta;
00440 float thit;
00441 Point phit;
00442
00443
00444 Ray ray;
00445 WorldToObject(r, &ray);
00446
00447 return FindIntersection(ray, &thit, &phit, &phi, &theta);
00448 }
00449
00450 float Torus::Area() const {
00451 return phiMax * majorRadius * (thetaMax - thetaMin) * minorRadius;
00452 }
00453
00454 Shape* Torus::CreateShape(const Transform &o2w,
00455 bool reverseOrientation,
00456 const ParamSet ¶ms) {
00457 float majorRadius = params.FindOneFloat("majorradius", 1.f);
00458 float minorRadius = params.FindOneFloat("minorradius", .25f);
00459 float thetamin = params.FindOneFloat("thetamin", 0);
00460 float thetamax = params.FindOneFloat("thetamax", 360.f);
00461 float phimax = params.FindOneFloat("phimax", 360.f);
00462 return new Torus(o2w, reverseOrientation, majorRadius, minorRadius,
00463 thetamin, thetamax, phimax);
00464 }
00465
00466 static DynamicLoader::RegisterShape<Torus> r("torus");