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