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