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/paraboloid.h"
00028 #include "paramset.h"
00029
00030
00031 Paraboloid::Paraboloid(const Transform *o2w, const Transform *w2o, bool ro,
00032 float rad, float z0, float z1,
00033 float tm)
00034 : Shape(o2w, w2o, ro) {
00035 radius = rad;
00036 zmin = min(z0,z1);
00037 zmax = max(z0,z1);
00038 phiMax = Radians(Clamp(tm, 0.0f, 360.0f));
00039 }
00040
00041
00042 BBox Paraboloid::ObjectBound() const {
00043 Point p1 = Point( -radius, -radius, zmin );
00044 Point p2 = Point( radius, radius, zmax );
00045 return BBox( p1, p2 );
00046 }
00047
00048
00049 bool Paraboloid::Intersect(const Ray &r, float *tHit,
00050 float *rayEpsilon, DifferentialGeometry *dg) const {
00051 float phi;
00052 Point phit;
00053
00054 Ray ray;
00055 (*WorldToObject)(r, &ray);
00056
00057
00058 float k = zmax/(radius*radius);
00059 float A = k*(ray.d.x * ray.d.x + ray.d.y * ray.d.y);
00060 float B = 2*k*(ray.d.x * ray.o.x + ray.d.y * ray.o.y) -
00061 ray.d.z;
00062 float C = k*(ray.o.x * ray.o.x + ray.o.y * ray.o.y) -
00063 ray.o.z;
00064
00065
00066 float t0, t1;
00067 if (!Quadratic(A, B, C, &t0, &t1))
00068 return false;
00069
00070
00071 if (t0 > ray.maxt || t1 < ray.mint)
00072 return false;
00073 float thit = t0;
00074 if (t0 < ray.mint) {
00075 thit = t1;
00076 if (thit > ray.maxt) return false;
00077 }
00078
00079
00080 phit = ray(thit);
00081 phi = atan2f(phit.y, phit.x);
00082 if (phi < 0.) phi += 2.f*M_PI;
00083
00084
00085 if (phit.z < zmin || phit.z > zmax || phi > phiMax) {
00086 if (thit == t1) return false;
00087 thit = t1;
00088 if (t1 > ray.maxt) return false;
00089
00090 phit = ray(thit);
00091 phi = atan2f(phit.y, phit.x);
00092 if (phi < 0.) phi += 2.f*M_PI;
00093 if (phit.z < zmin || phit.z > zmax || phi > phiMax)
00094 return false;
00095 }
00096
00097
00098 float u = phi / phiMax;
00099 float v = (phit.z-zmin) / (zmax-zmin);
00100
00101
00102 Vector dpdu(-phiMax * phit.y, phiMax * phit.x, 0.);
00103 Vector dpdv = (zmax - zmin) *
00104 Vector(phit.x / (2.f * phit.z), phit.y / (2.f * phit.z), 1.);
00105
00106
00107 Vector d2Pduu = -phiMax * phiMax *
00108 Vector(phit.x, phit.y, 0);
00109 Vector d2Pduv = (zmax - zmin) * phiMax *
00110 Vector(-phit.y / (2.f * phit.z),
00111 phit.x / (2.f * phit.z),
00112 0);
00113 Vector d2Pdvv = -(zmax - zmin) * (zmax - zmin) *
00114 Vector(phit.x / (4.f * phit.z * phit.z),
00115 phit.y / (4.f * phit.z * phit.z),
00116 0.);
00117
00118
00119 float E = Dot(dpdu, dpdu);
00120 float F = Dot(dpdu, dpdv);
00121 float G = Dot(dpdv, dpdv);
00122 Vector N = Normalize(Cross(dpdu, dpdv));
00123 float e = Dot(N, d2Pduu);
00124 float f = Dot(N, d2Pduv);
00125 float g = Dot(N, d2Pdvv);
00126
00127
00128 float invEGF2 = 1.f / (E*G - F*F);
00129 Normal dndu = Normal((f*F - e*G) * invEGF2 * dpdu +
00130 (e*F - f*E) * invEGF2 * dpdv);
00131 Normal dndv = Normal((g*F - f*G) * invEGF2 * dpdu +
00132 (f*F - g*E) * invEGF2 * dpdv);
00133
00134
00135 const Transform &o2w = *ObjectToWorld;
00136 *dg = DifferentialGeometry(o2w(phit), o2w(dpdu), o2w(dpdv),
00137 o2w(dndu), o2w(dndv), u, v, this);
00138
00139
00140 *tHit = thit;
00141
00142
00143 *rayEpsilon = 5e-4f * *tHit;
00144 return true;
00145 }
00146
00147
00148
00149 bool Paraboloid::IntersectP(const Ray &r) const {
00150 float phi;
00151 Point phit;
00152
00153 Ray ray;
00154 (*WorldToObject)(r, &ray);
00155
00156
00157 float k = zmax/(radius*radius);
00158 float A = k*(ray.d.x * ray.d.x + ray.d.y * ray.d.y);
00159 float B = 2*k*(ray.d.x * ray.o.x + ray.d.y * ray.o.y) -
00160 ray.d.z;
00161 float C = k*(ray.o.x * ray.o.x + ray.o.y * ray.o.y) -
00162 ray.o.z;
00163
00164
00165 float t0, t1;
00166 if (!Quadratic(A, B, C, &t0, &t1))
00167 return false;
00168
00169
00170 if (t0 > ray.maxt || t1 < ray.mint)
00171 return false;
00172 float thit = t0;
00173 if (t0 < ray.mint) {
00174 thit = t1;
00175 if (thit > ray.maxt) return false;
00176 }
00177
00178
00179 phit = ray(thit);
00180 phi = atan2f(phit.y, phit.x);
00181 if (phi < 0.) phi += 2.f*M_PI;
00182
00183
00184 if (phit.z < zmin || phit.z > zmax || phi > phiMax) {
00185 if (thit == t1) return false;
00186 thit = t1;
00187 if (t1 > ray.maxt) return false;
00188
00189 phit = ray(thit);
00190 phi = atan2f(phit.y, phit.x);
00191 if (phi < 0.) phi += 2.f*M_PI;
00192 if (phit.z < zmin || phit.z > zmax || phi > phiMax)
00193 return false;
00194 }
00195 return true;
00196 }
00197
00198
00199 float Paraboloid::Area() const {
00200 return phiMax/12.0f *
00201 (powf(1+4*zmin, 1.5f) - powf(1+4*zmax, 1.5f));
00202 }
00203
00204
00205 Paraboloid *CreateParaboloidShape(const Transform *o2w, const Transform *w2o,
00206 bool reverseOrientation, const ParamSet ¶ms) {
00207 float radius = params.FindOneFloat("radius", 1);
00208 float zmin = params.FindOneFloat("zmin", 0);
00209 float zmax = params.FindOneFloat("zmax", 1);
00210 float phimax = params.FindOneFloat("phimax", 360);
00211 return new Paraboloid(o2w, w2o, reverseOrientation, radius, zmin, zmax, phimax);
00212 }
00213
00214