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