00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 #include "stdafx.h"
00027 #include "shapes/hyperboloid.h"
00028 #include "paramset.h"
00029
00030
00031 Hyperboloid::Hyperboloid(const Transform *o2w, const Transform *w2o, bool ro,
00032 const Point &point1, const Point &point2, float tm)
00033 : Shape(o2w, w2o, ro) {
00034 p1 = point1;
00035 p2 = point2;
00036 phiMax = Radians(Clamp(tm, 0.0f, 360.0f));
00037 float radius1 = sqrtf(p1.x*p1.x + p1.y*p1.y);
00038 float radius2 = sqrtf(p2.x*p2.x + p2.y*p2.y);
00039 rmax = max(radius1, radius2);
00040 zmin = min(p1.z, p2.z);
00041 zmax = max(p1.z, p2.z);
00042
00043 if (p2.z == 0.) swap(p1, p2);
00044 Point pp = p1;
00045 float xy1, xy2;
00046 do {
00047 pp += 2.f * (p2-p1);
00048 xy1 = pp.x*pp.x + pp.y*pp.y;
00049 xy2 = p2.x*p2.x + p2.y*p2.y;
00050 a = (1.f/xy1 - (pp.z*pp.z)/(xy1*p2.z*p2.z)) /
00051 (1 - (xy2*pp.z*pp.z)/(xy1*p2.z*p2.z));
00052 c = (a * xy2 - 1) / (p2.z*p2.z);
00053 } while (isinf(a) || isnan(a));
00054 }
00055
00056
00057 BBox Hyperboloid::ObjectBound() const {
00058 Point p1 = Point( -rmax, -rmax, zmin );
00059 Point p2 = Point( rmax, rmax, zmax );
00060 return BBox( p1, p2 );
00061 }
00062
00063
00064 bool Hyperboloid::Intersect(const Ray &r, float *tHit,
00065 float *rayEpsilon, DifferentialGeometry *dg) const {
00066 float phi, v;
00067 Point phit;
00068
00069 Ray ray;
00070 (*WorldToObject)(r, &ray);
00071
00072
00073 float A = a*ray.d.x*ray.d.x +
00074 a*ray.d.y*ray.d.y -
00075 c*ray.d.z*ray.d.z;
00076 float B = 2.f * (a*ray.d.x*ray.o.x +
00077 a*ray.d.y*ray.o.y -
00078 c*ray.d.z*ray.o.z);
00079 float C = a*ray.o.x*ray.o.x +
00080 a*ray.o.y*ray.o.y -
00081 c*ray.o.z*ray.o.z - 1;
00082
00083
00084 float t0, t1;
00085 if (!Quadratic(A, B, C, &t0, &t1))
00086 return false;
00087
00088
00089 if (t0 > ray.maxt || t1 < ray.mint)
00090 return false;
00091 float thit = t0;
00092 if (t0 < ray.mint) {
00093 thit = t1;
00094 if (thit > ray.maxt) return false;
00095 }
00096
00097
00098 phit = ray(thit);
00099 v = (phit.z - p1.z)/(p2.z - p1.z);
00100 Point pr = (1.f-v) * p1 + v * p2;
00101 phi = atan2f(pr.x*phit.y - phit.x*pr.y,
00102 phit.x*pr.x + phit.y*pr.y);
00103 if (phi < 0)
00104 phi += 2*M_PI;
00105
00106
00107 if (phit.z < zmin || phit.z > zmax || phi > phiMax) {
00108 if (thit == t1) return false;
00109 thit = t1;
00110 if (t1 > ray.maxt) return false;
00111
00112 phit = ray(thit);
00113 v = (phit.z - p1.z)/(p2.z - p1.z);
00114 Point pr = (1.f-v) * p1 + v * p2;
00115 phi = atan2f(pr.x*phit.y - phit.x*pr.y,
00116 phit.x*pr.x + phit.y*pr.y);
00117 if (phi < 0)
00118 phi += 2*M_PI;
00119 if (phit.z < zmin || phit.z > zmax || phi > phiMax)
00120 return false;
00121 }
00122
00123
00124 float u = phi / phiMax;
00125
00126
00127 float cosphi = cosf(phi), sinphi = sinf(phi);
00128 Vector dpdu(-phiMax * phit.y, phiMax * phit.x, 0.);
00129 Vector dpdv((p2.x-p1.x) * cosphi - (p2.y-p1.y) * sinphi,
00130 (p2.x-p1.x) * sinphi + (p2.y-p1.y) * cosphi,
00131 p2.z-p1.z);
00132
00133
00134 Vector d2Pduu = -phiMax * phiMax *
00135 Vector(phit.x, phit.y, 0);
00136 Vector d2Pduv = phiMax *
00137 Vector(-dpdv.y, dpdv.x, 0.);
00138 Vector d2Pdvv(0, 0, 0);
00139
00140
00141 float E = Dot(dpdu, dpdu);
00142 float F = Dot(dpdu, dpdv);
00143 float G = Dot(dpdv, dpdv);
00144 Vector N = Normalize(Cross(dpdu, dpdv));
00145 float e = Dot(N, d2Pduu);
00146 float f = Dot(N, d2Pduv);
00147 float g = Dot(N, d2Pdvv);
00148
00149
00150 float invEGF2 = 1.f / (E*G - F*F);
00151 Normal dndu = Normal((f*F - e*G) * invEGF2 * dpdu +
00152 (e*F - f*E) * invEGF2 * dpdv);
00153 Normal dndv = Normal((g*F - f*G) * invEGF2 * dpdu +
00154 (f*F - g*E) * invEGF2 * dpdv);
00155
00156
00157 const Transform &o2w = *ObjectToWorld;
00158 *dg = DifferentialGeometry(o2w(phit), o2w(dpdu), o2w(dpdv),
00159 o2w(dndu), o2w(dndv), u, v, this);
00160
00161
00162 *tHit = thit;
00163
00164
00165 *rayEpsilon = 5e-4f * *tHit;
00166 return true;
00167 }
00168
00169
00170 bool Hyperboloid::IntersectP(const Ray &r) const {
00171 float phi, v;
00172 Point phit;
00173
00174 Ray ray;
00175 (*WorldToObject)(r, &ray);
00176
00177
00178 float A = a*ray.d.x*ray.d.x +
00179 a*ray.d.y*ray.d.y -
00180 c*ray.d.z*ray.d.z;
00181 float B = 2.f * (a*ray.d.x*ray.o.x +
00182 a*ray.d.y*ray.o.y -
00183 c*ray.d.z*ray.o.z);
00184 float C = a*ray.o.x*ray.o.x +
00185 a*ray.o.y*ray.o.y -
00186 c*ray.o.z*ray.o.z - 1;
00187
00188
00189 float t0, t1;
00190 if (!Quadratic(A, B, C, &t0, &t1))
00191 return false;
00192
00193
00194 if (t0 > ray.maxt || t1 < ray.mint)
00195 return false;
00196 float thit = t0;
00197 if (t0 < ray.mint) {
00198 thit = t1;
00199 if (thit > ray.maxt) return false;
00200 }
00201
00202
00203 phit = ray(thit);
00204 v = (phit.z - p1.z)/(p2.z - p1.z);
00205 Point pr = (1.f-v) * p1 + v * p2;
00206 phi = atan2f(pr.x*phit.y - phit.x*pr.y,
00207 phit.x*pr.x + phit.y*pr.y);
00208 if (phi < 0)
00209 phi += 2*M_PI;
00210
00211
00212 if (phit.z < zmin || phit.z > zmax || phi > phiMax) {
00213 if (thit == t1) return false;
00214 thit = t1;
00215 if (t1 > ray.maxt) return false;
00216
00217 phit = ray(thit);
00218 v = (phit.z - p1.z)/(p2.z - p1.z);
00219 Point pr = (1.f-v) * p1 + v * p2;
00220 phi = atan2f(pr.x*phit.y - phit.x*pr.y,
00221 phit.x*pr.x + phit.y*pr.y);
00222 if (phi < 0)
00223 phi += 2*M_PI;
00224 if (phit.z < zmin || phit.z > zmax || phi > phiMax)
00225 return false;
00226 }
00227 return true;
00228 }
00229
00230
00231 #define SQR(a) ((a)*(a))
00232 #define QUAD(a) ((SQR(a))*(SQR(a)))
00233 float Hyperboloid::Area() const {
00234 return phiMax/6.f *
00235 (2.f*QUAD(p1.x) - 2.f*p1.x*p1.x*p1.x*p2.x +
00236 2.f*QUAD(p2.x) +
00237 2.f*(p1.y*p1.y + p1.y*p2.y + p2.y*p2.y)*
00238 (SQR(p1.y - p2.y) + SQR(p1.z - p2.z)) +
00239 p2.x*p2.x*(5.f*p1.y*p1.y + 2.f*p1.y*p2.y -
00240 4.f*p2.y*p2.y + 2.f*SQR(p1.z - p2.z)) +
00241 p1.x*p1.x*(-4.f*p1.y*p1.y + 2.f*p1.y*p2.y +
00242 5.f*p2.y*p2.y + 2.f*SQR(p1.z - p2.z)) -
00243 2.f*p1.x*p2.x*(p2.x*p2.x - p1.y*p1.y +
00244 5.f*p1.y*p2.y - p2.y*p2.y - p1.z*p1.z +
00245 2.f*p1.z*p2.z - p2.z*p2.z));
00246 }
00247
00248
00249 #undef SQR
00250 #undef QUAD
00251 Shape *CreateHyperboloidShape(const Transform *o2w, const Transform *w2o,
00252 bool reverseOrientation, const ParamSet ¶ms) {
00253 Point p1 = params.FindOnePoint("p1", Point(0,0,0));
00254 Point p2 = params.FindOnePoint("p2", Point(1,1,1));
00255 float phimax = params.FindOneFloat("phimax", 360);
00256 return new Hyperboloid(o2w, w2o, reverseOrientation, p1, p2, phimax);
00257 }
00258
00259