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