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 #include "mc.h"
00027
00028 class Sphere: public Shape {
00029 public:
00030
00031 Sphere(const Transform &o2w, bool ro, float rad,
00032 float zmin, float zmax, float phiMax);
00033 BBox ObjectBound() const;
00034 bool Intersect(const Ray &ray, float *tHit,
00035 DifferentialGeometry *dg) const;
00036 bool IntersectP(const Ray &ray) const;
00037 float Area() const;
00038 Point Sample(float u1, float u2, Normal *ns) const {
00039 Point p = Point(0,0,0) + radius *
00040 UniformSampleSphere(u1, u2);
00041 *ns = Normalize(ObjectToWorld(Normal(p.x, p.y, p.z)));
00042 if (reverseOrientation) *ns *= -1.f;
00043 return ObjectToWorld(p);
00044 }
00045 Point Sample(const Point &p,
00046 float u1, float u2, Normal *ns) const {
00047
00048 Point Pcenter = ObjectToWorld(Point(0,0,0));
00049 Vector wc = Normalize(Pcenter - p);
00050 Vector wcX, wcY;
00051 CoordinateSystem(wc, &wcX, &wcY);
00052
00053 if (DistanceSquared(p, Pcenter) - radius*radius < 1e-4f)
00054 return Sample(u1, u2, ns);
00055
00056 float cosThetaMax = sqrtf(max(0.f, 1.f - radius*radius /
00057 DistanceSquared(p, Pcenter)));
00058 DifferentialGeometry dgSphere;
00059 float thit;
00060 Point ps;
00061 Ray r(p,
00062 UniformSampleCone(u1, u2, cosThetaMax, wcX, wcY, wc));
00063 if (!Intersect(r, &thit, &dgSphere)) {
00064 ps = Pcenter - radius * wc;
00065 } else {
00066 ps = r(thit);
00067 }
00068 *ns = Normal(Normalize(ps - Pcenter));
00069 if (reverseOrientation) *ns *= -1.f;
00070 return ps;
00071 }
00072 float Pdf(const Point &p, const Vector &wi) const {
00073 Point Pcenter = ObjectToWorld(Point(0,0,0));
00074
00075 if (DistanceSquared(p, Pcenter) - radius*radius < 1e-4f)
00076 return Shape::Pdf(p, wi);
00077
00078 float cosThetaMax = sqrtf(max(0.f, 1.f - radius*radius /
00079 DistanceSquared(p, Pcenter)));
00080 return UniformConePdf(cosThetaMax);
00081 }
00082 private:
00083
00084 float radius;
00085 float phiMax;
00086 float zmin, zmax;
00087 float thetaMin, thetaMax;
00088 };
00089
00090 Sphere::Sphere(const Transform &o2w, bool ro, float rad,
00091 float z0, float z1, float pm)
00092 : Shape(o2w, ro) {
00093 radius = rad;
00094 zmin = Clamp(min(z0, z1), -radius, radius);
00095 zmax = Clamp(max(z0, z1), -radius, radius);
00096 thetaMin = acosf(Clamp(zmin/radius, -1.f, 1.f));
00097 thetaMax = acosf(Clamp(zmax/radius, -1.f, 1.f));
00098 phiMax = Radians(Clamp(pm, 0.0f, 360.0f));
00099 }
00100 BBox Sphere::ObjectBound() const {
00101 return BBox(Point(-radius, -radius, zmin),
00102 Point(radius, radius, zmax));
00103 }
00104 bool Sphere::Intersect(const Ray &r, float *tHit,
00105 DifferentialGeometry *dg) const {
00106 float phi;
00107 Point phit;
00108
00109 Ray ray;
00110 WorldToObject(r, &ray);
00111
00112 float A = ray.d.x*ray.d.x + ray.d.y*ray.d.y +
00113 ray.d.z*ray.d.z;
00114 float B = 2 * (ray.d.x*ray.o.x + ray.d.y*ray.o.y +
00115 ray.d.z*ray.o.z);
00116 float C = ray.o.x*ray.o.x + ray.o.y*ray.o.y +
00117 ray.o.z*ray.o.z - radius*radius;
00118
00119 float t0, t1;
00120 if (!Quadratic(A, B, C, &t0, &t1))
00121 return false;
00122
00123 if (t0 > ray.maxt || t1 < ray.mint)
00124 return false;
00125 float thit = t0;
00126 if (t0 < ray.mint) {
00127 thit = t1;
00128 if (thit > ray.maxt) return false;
00129 }
00130
00131 phit = ray(thit);
00132 phi = atan2f(phit.y, phit.x);
00133 if (phi < 0.) phi += 2.f*M_PI;
00134
00135 if (phit.z < zmin || phit.z > zmax || phi > phiMax) {
00136 if (thit == t1) return false;
00137 if (t1 > ray.maxt) return false;
00138 thit = t1;
00139
00140 phit = ray(thit);
00141 phi = atan2f(phit.y, phit.x);
00142 if (phi < 0.) phi += 2.f*M_PI;
00143 if (phit.z < zmin || phit.z > zmax || phi > phiMax)
00144 return false;
00145 }
00146
00147 float u = phi / phiMax;
00148 float theta = acosf(phit.z / radius);
00149 float v = (theta - thetaMin) / (thetaMax - thetaMin);
00150
00151 float cosphi, sinphi;
00152 Vector dpdu, dpdv;
00153 float zradius = sqrtf(phit.x*phit.x + phit.y*phit.y);
00154 if (zradius == 0)
00155 {
00156
00157 cosphi = 0;
00158 sinphi = 1;
00159 dpdv = (thetaMax-thetaMin) *
00160 Vector(phit.z * cosphi, phit.z * sinphi,
00161 -radius * sinf(theta));
00162 Vector norm = Vector(phit);
00163 dpdu = Cross(dpdv, norm);
00164 }
00165 else
00166 {
00167 float invzradius = 1.f / zradius;
00168 cosphi = phit.x * invzradius;
00169 sinphi = phit.y * invzradius;
00170 dpdu = Vector(-phiMax * phit.y, phiMax * phit.x, 0);
00171 dpdv = (thetaMax-thetaMin) *
00172 Vector(phit.z * cosphi, phit.z * sinphi,
00173 -radius * sinf(theta));
00174 }
00175
00176 Vector d2Pduu = -phiMax * phiMax * Vector(phit.x,phit.y,0);
00177 Vector d2Pduv = (thetaMax - thetaMin) *
00178 phit.z * phiMax *
00179 Vector(-sinphi, cosphi, 0.);
00180 Vector d2Pdvv = -(thetaMax - thetaMin) *
00181 (thetaMax - thetaMin) *
00182 Vector(phit.x, phit.y, phit.z);
00183
00184 float E = Dot(dpdu, dpdu);
00185 float F = Dot(dpdu, dpdv);
00186 float G = Dot(dpdv, dpdv);
00187 Vector N = Normalize(Cross(dpdu, dpdv));
00188 float e = Dot(N, d2Pduu);
00189 float f = Dot(N, d2Pduv);
00190 float g = Dot(N, d2Pdvv);
00191
00192 float invEGF2 = 1.f / (E*G - F*F);
00193 Vector dndu = (f*F - e*G) * invEGF2 * dpdu +
00194 (e*F - f*E) * invEGF2 * dpdv;
00195 Vector dndv = (g*F - f*G) * invEGF2 * dpdu +
00196 (f*F - g*E) * invEGF2 * dpdv;
00197
00198 *dg = DifferentialGeometry(ObjectToWorld(phit),
00199 ObjectToWorld(dpdu),
00200 ObjectToWorld(dpdv),
00201 ObjectToWorld(dndu),
00202 ObjectToWorld(dndv),
00203 u, v, this);
00204
00205 *tHit = thit;
00206 return true;
00207 }
00208 bool Sphere::IntersectP(const Ray &r) const {
00209 float phi;
00210 Point phit;
00211
00212 Ray ray;
00213 WorldToObject(r, &ray);
00214
00215 float A = ray.d.x*ray.d.x + ray.d.y*ray.d.y +
00216 ray.d.z*ray.d.z;
00217 float B = 2 * (ray.d.x*ray.o.x + ray.d.y*ray.o.y +
00218 ray.d.z*ray.o.z);
00219 float C = ray.o.x*ray.o.x + ray.o.y*ray.o.y +
00220 ray.o.z*ray.o.z - radius*radius;
00221
00222 float t0, t1;
00223 if (!Quadratic(A, B, C, &t0, &t1))
00224 return false;
00225
00226 if (t0 > ray.maxt || t1 < ray.mint)
00227 return false;
00228 float thit = t0;
00229 if (t0 < ray.mint) {
00230 thit = t1;
00231 if (thit > ray.maxt) return false;
00232 }
00233
00234 phit = ray(thit);
00235 phi = atan2f(phit.y, phit.x);
00236 if (phi < 0.) phi += 2.f*M_PI;
00237
00238 if (phit.z < zmin || phit.z > zmax || phi > phiMax) {
00239 if (thit == t1) return false;
00240 if (t1 > ray.maxt) return false;
00241 thit = t1;
00242
00243 phit = ray(thit);
00244 phi = atan2f(phit.y, phit.x);
00245 if (phi < 0.) phi += 2.f*M_PI;
00246 if (phit.z < zmin || phit.z > zmax || phi > phiMax)
00247 return false;
00248 }
00249 return true;
00250 }
00251 float Sphere::Area() const {
00252 return phiMax * radius * (zmax-zmin);
00253 }
00254 extern "C" DLLEXPORT
00255 Shape *CreateShape(const Transform &o2w,
00256 bool reverseOrientation,
00257 const ParamSet ¶ms) {
00258 float radius = params.FindOneFloat("radius", 1.f);
00259 float zmin = params.FindOneFloat("zmin", -radius);
00260 float zmax = params.FindOneFloat("zmax", radius);
00261 float phimax = params.FindOneFloat("phimax", 360.f);
00262 return new Sphere(o2w, reverseOrientation, radius,
00263 zmin, zmax, phimax);
00264 }