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