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