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