00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024 #include "sphere.h"
00025 #include "paramset.h"
00026 #include "dynload.h"
00027
00028 using namespace lux;
00029
00030
00031 Sphere::Sphere(const Transform &o2w, bool ro, float rad,
00032 float z0, float z1, float pm)
00033 : Shape(o2w, ro) {
00034 radius = rad;
00035 zmin = Clamp(min(z0, z1), -radius, radius);
00036 zmax = Clamp(max(z0, z1), -radius, radius);
00037 thetaMin = acosf(Clamp(zmin/radius, -1.f, 1.f));
00038 thetaMax = acosf(Clamp(zmax/radius, -1.f, 1.f));
00039 phiMax = Radians(Clamp(pm, 0.0f, 360.0f));
00040 }
00041 BBox Sphere::ObjectBound() const {
00042 return BBox(Point(-radius, -radius, zmin),
00043 Point(radius, radius, zmax));
00044 }
00045 bool Sphere::Intersect(const Ray &r, float *tHit,
00046 DifferentialGeometry *dg) const {
00047 float phi;
00048 Point phit;
00049
00050 Ray ray;
00051 WorldToObject(r, &ray);
00052
00053 float A = ray.d.x*ray.d.x + ray.d.y*ray.d.y +
00054 ray.d.z*ray.d.z;
00055 float B = 2 * (ray.d.x*ray.o.x + ray.d.y*ray.o.y +
00056 ray.d.z*ray.o.z);
00057 float C = ray.o.x*ray.o.x + ray.o.y*ray.o.y +
00058 ray.o.z*ray.o.z - radius*radius;
00059
00060 float t0, t1;
00061 if (!Quadratic(A, B, C, &t0, &t1))
00062 return false;
00063
00064 if (t0 > ray.maxt || t1 < ray.mint)
00065 return false;
00066 float thit = t0;
00067 if (t0 < ray.mint) {
00068 thit = t1;
00069 if (thit > ray.maxt) return false;
00070 }
00071
00072 phit = ray(thit);
00073 phi = atan2f(phit.y, phit.x);
00074 if (phi < 0.) phi += 2.f*M_PI;
00075
00076 if (phit.z < zmin || phit.z > zmax || phi > phiMax) {
00077 if (thit == t1) return false;
00078 if (t1 > ray.maxt) return false;
00079 thit = t1;
00080
00081 phit = ray(thit);
00082 phi = atan2f(phit.y, phit.x);
00083 if (phi < 0.) phi += 2.f*M_PI;
00084 if (phit.z < zmin || phit.z > zmax || phi > phiMax)
00085 return false;
00086 }
00087
00088 float u = phi / phiMax;
00089 float theta = acosf(phit.z / radius);
00090 float v = (theta - thetaMin) / (thetaMax - thetaMin);
00091
00092 float cosphi, sinphi;
00093 Vector dpdu, dpdv;
00094 float zradius = sqrtf(phit.x*phit.x + phit.y*phit.y);
00095 if (zradius == 0)
00096 {
00097
00098 cosphi = 0;
00099 sinphi = 1;
00100 dpdv = (thetaMax-thetaMin) *
00101 Vector(phit.z * cosphi, phit.z * sinphi,
00102 -radius * sinf(theta));
00103 Vector norm = Vector(phit);
00104 dpdu = Cross(dpdv, norm);
00105 }
00106 else
00107 {
00108 float invzradius = 1.f / zradius;
00109 cosphi = phit.x * invzradius;
00110 sinphi = phit.y * invzradius;
00111 dpdu = Vector(-phiMax * phit.y, phiMax * phit.x, 0);
00112 dpdv = (thetaMax-thetaMin) *
00113 Vector(phit.z * cosphi, phit.z * sinphi,
00114 -radius * sinf(theta));
00115 }
00116
00117 Vector d2Pduu = -phiMax * phiMax * Vector(phit.x,phit.y,0);
00118 Vector d2Pduv = (thetaMax - thetaMin) *
00119 phit.z * phiMax *
00120 Vector(-sinphi, cosphi, 0.);
00121 Vector d2Pdvv = -(thetaMax - thetaMin) *
00122 (thetaMax - thetaMin) *
00123 Vector(phit.x, phit.y, phit.z);
00124
00125 float E = Dot(dpdu, dpdu);
00126 float F = Dot(dpdu, dpdv);
00127 float G = Dot(dpdv, dpdv);
00128 Vector N = Normalize(Cross(dpdu, dpdv));
00129 float e = Dot(N, d2Pduu);
00130 float f = Dot(N, d2Pduv);
00131 float g = Dot(N, d2Pdvv);
00132
00133 float invEGF2 = 1.f / (E*G - F*F);
00134 Normal dndu((f*F - e*G) * invEGF2 * dpdu +
00135 (e*F - f*E) * invEGF2 * dpdv);
00136 Normal dndv((g*F - f*G) * invEGF2 * dpdu +
00137 (f*F - g*E) * invEGF2 * dpdv);
00138
00139 *dg = DifferentialGeometry(ObjectToWorld(phit),
00140 ObjectToWorld(dpdu),
00141 ObjectToWorld(dpdv),
00142 ObjectToWorld(dndu),
00143 ObjectToWorld(dndv),
00144 u, v, this);
00145
00146 *tHit = thit;
00147 return true;
00148 }
00149 bool Sphere::IntersectP(const Ray &r) const {
00150 float phi;
00151 Point phit;
00152
00153 Ray ray;
00154 WorldToObject(r, &ray);
00155
00156 float A = ray.d.x*ray.d.x + ray.d.y*ray.d.y +
00157 ray.d.z*ray.d.z;
00158 float B = 2 * (ray.d.x*ray.o.x + ray.d.y*ray.o.y +
00159 ray.d.z*ray.o.z);
00160 float C = ray.o.x*ray.o.x + ray.o.y*ray.o.y +
00161 ray.o.z*ray.o.z - radius*radius;
00162
00163 float t0, t1;
00164 if (!Quadratic(A, B, C, &t0, &t1))
00165 return false;
00166
00167 if (t0 > ray.maxt || t1 < ray.mint)
00168 return false;
00169 float thit = t0;
00170 if (t0 < ray.mint) {
00171 thit = t1;
00172 if (thit > ray.maxt) return false;
00173 }
00174
00175 phit = ray(thit);
00176 phi = atan2f(phit.y, phit.x);
00177 if (phi < 0.) phi += 2.f*M_PI;
00178
00179 if (phit.z < zmin || phit.z > zmax || phi > phiMax) {
00180 if (thit == t1) return false;
00181 if (t1 > ray.maxt) return false;
00182 thit = t1;
00183
00184 phit = ray(thit);
00185 phi = atan2f(phit.y, phit.x);
00186 if (phi < 0.) phi += 2.f*M_PI;
00187 if (phit.z < zmin || phit.z > zmax || phi > phiMax)
00188 return false;
00189 }
00190 return true;
00191 }
00192 float Sphere::Area() const {
00193 return phiMax * radius * (zmax-zmin);
00194 }
00195 Shape* Sphere::CreateShape(const Transform &o2w,
00196 bool reverseOrientation,
00197 const ParamSet ¶ms) {
00198 float radius = params.FindOneFloat("radius", 1.f);
00199 float zmin = params.FindOneFloat("zmin", -radius);
00200 float zmax = params.FindOneFloat("zmax", radius);
00201 float phimax = params.FindOneFloat("phimax", 360.f);
00202 return new Sphere(o2w, reverseOrientation, radius,
00203 zmin, zmax, phimax);
00204 }
00205
00206 static DynamicLoader::RegisterShape<Sphere> r("sphere");