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 "sky.h"
00025 #include "mc.h"
00026 #include "spectrumwavelengths.h"
00027 #include "paramset.h"
00028 #include "regular.h"
00029 #include "reflection/bxdf.h"
00030 #include "dynload.h"
00031
00032 #include "data/skychroma_spect.h"
00033
00034 using namespace lux;
00035
00036 class SkyBxDF : public BxDF
00037 {
00038 public:
00039 SkyBxDF(const SkyLight &sky, const Transform &WL, const Vector &x, const Vector &y, const Vector &z) : BxDF(BxDFType(BSDF_REFLECTION | BSDF_DIFFUSE)), skyLight(sky), WorldToLight(WL), X(x), Y(y), Z(z) {}
00040 virtual ~SkyBxDF() { }
00041 virtual void f(const TsPack *tspack, const Vector &wo, const Vector &wi, SWCSpectrum *const f) const
00042 {
00043 Vector w(wi.x * X + wi.y * Y + wi.z * Z);
00044 w = -Normalize(WorldToLight(w));
00045 const float phi = SphericalPhi(w);
00046 const float theta = SphericalTheta(w);
00047 SWCSpectrum L;
00048 skyLight.GetSkySpectralRadiance(tspack, theta, phi, &L);
00049 *f += L;
00050 }
00051 private:
00052 const SkyLight &skyLight;
00053 const Transform &WorldToLight;
00054 Vector X, Y, Z;
00055 };
00056
00057 class SkyPortalBxDF : public BxDF
00058 {
00059 public:
00060 SkyPortalBxDF(const SkyLight &sky, const Transform &WL, const Vector &x, const Vector &y, const Vector &z, const Point &p, const vector<boost::shared_ptr<Primitive> > &portalList, u_int portal, float u) : BxDF(BxDFType(BSDF_REFLECTION | BSDF_DIFFUSE)), skyLight(sky), WorldToLight(WL), X(x), Y(y), Z(z), ps(p), PortalShapes(portalList), shapeIndex(portal), u3(u) {}
00061 virtual ~SkyPortalBxDF() { }
00062 virtual bool Sample_f(const TsPack *tspack, const Vector &wo, Vector *wi, float u1, float u2, SWCSpectrum *const f,float *pdf, float *pdfBack = NULL, bool reverse = false) const
00063 {
00064 if (shapeIndex == ~0U || reverse)
00065 return false;
00066 DifferentialGeometry dg;
00067 dg.time = tspack->time;
00068 PortalShapes[shapeIndex]->Sample(ps, u1, u2, u3, &dg);
00069 Vector wiW = Normalize(dg.p - ps);
00070 Vector w = -Normalize(WorldToLight(wiW));
00071 const float phi = SphericalPhi(w);
00072 const float theta = SphericalTheta(w);
00073 skyLight.GetSkySpectralRadiance(tspack, theta, phi, f);
00074 wi->x = Dot(wiW, X);
00075 wi->y = Dot(wiW, Y);
00076 wi->z = Dot(wiW, Z);
00077 *wi = Normalize(*wi);
00078 *pdf = PortalShapes[shapeIndex]->Pdf(ps, dg.p) * DistanceSquared(ps, dg.p) / AbsDot(wiW, dg.nn);
00079 for (u_int i = 0; i < PortalShapes.size(); ++i) {
00080 if (i == shapeIndex)
00081 continue;
00082 Intersection isect;
00083 RayDifferential ray(ps, wiW);
00084 ray.mint = -INFINITY;
00085 if (PortalShapes[i]->Intersect(ray, &isect) && Dot(wiW, isect.dg.nn) > 0.f)
00086 *pdf += PortalShapes[i]->Pdf(ps, isect.dg.p) * DistanceSquared(ps, isect.dg.p) / AbsDot(wiW, isect.dg.nn);
00087 }
00088 *pdf /= PortalShapes.size();
00089 if (pdfBack)
00090 *pdfBack = 0.f;
00091 return true;
00092 }
00093 virtual void f(const TsPack *tspack, const Vector &wo, const Vector &wi, SWCSpectrum *const f) const
00094 {
00095 const Vector w(-Normalize(WorldToLight(wi.x * X + wi.y * Y + wi.z * Z)));
00096 const float phi = SphericalPhi(w);
00097 const float theta = SphericalTheta(w);
00098 SWCSpectrum L;
00099 skyLight.GetSkySpectralRadiance(tspack, theta, phi, &L);
00100 *f += L;
00101 }
00102 virtual float Pdf(const TsPack *tspack, const Vector &wi, const Vector &wo) const
00103 {
00104 const Vector w(wo.x * X + wo.y * Y + wo.z * Z);
00105 float pdf = 0.f;
00106 for (u_int i = 0; i < PortalShapes.size(); ++i) {
00107 Intersection isect;
00108 RayDifferential ray(ps, w);
00109 ray.mint = -INFINITY;
00110 if (PortalShapes[i]->Intersect(ray, &isect) && Dot(w, isect.dg.nn) > 0.f)
00111 pdf += PortalShapes[i]->Pdf(ps, isect.dg.p) * DistanceSquared(ps, isect.dg.p) / AbsDot(w, isect.dg.nn);
00112 }
00113 return pdf / PortalShapes.size();
00114 }
00115 private:
00116 const SkyLight &skyLight;
00117 const Transform &WorldToLight;
00118 Vector X, Y, Z;
00119 Point ps;
00120 const vector<boost::shared_ptr<Primitive> > &PortalShapes;
00121 u_int shapeIndex;
00122 float u3;
00123 };
00124
00125 static float PerezBase(const float lam[6], float theta, float gamma)
00126 {
00127 return (1.f + lam[1] * exp(lam[2] / cos(theta))) *
00128 (1.f + lam[3] * exp(lam[4] * gamma) + lam[5] * cos(gamma) * cos(gamma));
00129 }
00130
00131 static const RegularSPD S0(S0Amplitudes, 300.f, 830.f, 54);
00132 static const RegularSPD S1(S1Amplitudes, 300.f, 830.f, 54);
00133 static const RegularSPD S2(S2Amplitudes, 300.f, 830.f, 54);
00134 static const float S0Y = S0.Y();
00135 static const float S1Y = S1.Y();
00136 static const float S2Y = S2.Y();
00137
00138
00139 SkyLight::~SkyLight() {
00140 }
00141
00142 SkyLight::SkyLight(const Transform &light2world,
00143 const float skyscale, int ns, Vector sd, float turb,
00144 float aconst, float bconst, float cconst, float dconst, float econst)
00145 : Light(light2world, ns) {
00146 skyScale = skyscale;
00147 sundir = sd;
00148 turbidity = turb;
00149
00150 InitSunThetaPhi();
00151
00152 float theta2 = thetaS*thetaS;
00153 float theta3 = theta2*thetaS;
00154 float T = turb;
00155 float T2 = turb*turb;
00156
00157 float chi = (4.f / 9.f - T / 120.f) * (M_PI - 2 * thetaS);
00158 zenith_Y = (4.0453 * T - 4.9710) * tan(chi) - .2155 * T + 2.4192;
00159 zenith_Y *= 1000;
00160
00161 zenith_x =
00162 (0.00166f * theta3 - 0.00375f * theta2 + 0.00209f * thetaS) * T2 +
00163 (-0.02903f * theta3 + 0.06377f * theta2 - 0.03202f * thetaS + 0.00394f) * T +
00164 (0.11693f * theta3 - 0.21196f * theta2 + 0.06052f * thetaS + 0.25886f);
00165
00166 zenith_y =
00167 (0.00275f * theta3 - 0.00610f * theta2 + 0.00317f * thetaS) * T2 +
00168 (-0.04214f * theta3 + 0.08970f * theta2 - 0.04153f * thetaS + 0.00516f) * T +
00169 (0.15346f * theta3 - 0.26756f * theta2 + 0.06670f * thetaS + 0.26688f);
00170
00171 perez_Y[1] = (0.1787f * T - 1.4630f) * aconst;
00172 perez_Y[2] = (-0.3554f * T + 0.4275f) * bconst;
00173 perez_Y[3] = (-0.0227f * T + 5.3251f) * cconst;
00174 perez_Y[4] = (0.1206f * T - 2.5771f) * dconst;
00175 perez_Y[5] = (-0.0670f * T + 0.3703f) * econst;
00176
00177 perez_x[1] = (-0.0193f * T - 0.2592f) * aconst;
00178 perez_x[2] = (-0.0665f * T + 0.0008f) * bconst;
00179 perez_x[3] = (-0.0004f * T + 0.2125f) * cconst;
00180 perez_x[4] = (-0.0641f * T - 0.8989f) * dconst;
00181 perez_x[5] = (-0.0033f * T + 0.0452f) * econst;
00182
00183 perez_y[1] = (-0.0167f * T - 0.2608f) * aconst;
00184 perez_y[2] = (-0.0950f * T + 0.0092f) * bconst;
00185 perez_y[3] = (-0.0079f * T + 0.2102f) * cconst;
00186 perez_y[4] = (-0.0441f * T - 1.6537f) * dconst;
00187 perez_y[5] = (-0.0109f * T + 0.0529f) * econst;
00188
00189 zenith_Y /= PerezBase(perez_Y, 0, thetaS);
00190 zenith_x /= PerezBase(perez_x, 0, thetaS);
00191 zenith_y /= PerezBase(perez_y, 0, thetaS);
00192 }
00193
00194 SWCSpectrum SkyLight::Le(const TsPack *tspack, const RayDifferential &r) const {
00195 Vector w = r.d;
00196
00197
00198 Vector wh = Normalize(WorldToLight(w));
00199 const float phi = SphericalPhi(wh);
00200 const float theta = SphericalTheta(wh);
00201
00202 SWCSpectrum L;
00203 GetSkySpectralRadiance(tspack, theta,phi,(SWCSpectrum * const)&L);
00204 L *= skyScale;
00205
00206 return L;
00207 }
00208 SWCSpectrum SkyLight::Le(const TsPack *tspack, const Scene *scene, const Ray &r,
00209 const Normal &n, BSDF **bsdf, float *pdf, float *pdfDirect) const
00210 {
00211 Point worldCenter;
00212 float worldRadius;
00213 scene->WorldBound().BoundingSphere(&worldCenter, &worldRadius);
00214 Vector toCenter(worldCenter - r.o);
00215 float centerDistance = Dot(toCenter, toCenter);
00216 float approach = Dot(toCenter, r.d);
00217 float distance = approach + sqrtf(worldRadius * worldRadius - centerDistance + approach * approach);
00218 Point ps(r.o + distance * r.d);
00219 Normal ns(Normalize(worldCenter - ps));
00220 Vector dpdu, dpdv;
00221 CoordinateSystem(Vector(ns), &dpdu, &dpdv);
00222 DifferentialGeometry dg(ps, ns, dpdu, dpdv, Normal(0, 0, 0), Normal (0, 0, 0), 0, 0, NULL);
00223 dg.time = tspack->time;
00224 if (!havePortalShape) {
00225 *bsdf = BSDF_ALLOC(tspack, SingleBSDF)(dg, ns,
00226 BSDF_ALLOC(tspack, SkyBxDF)(*this, WorldToLight, dpdu, dpdv, Vector(ns)));
00227 *pdf = 1.f / (4.f * M_PI * worldRadius * worldRadius);
00228 *pdfDirect = AbsDot(r.d, n) * INV_TWOPI * AbsDot(r.d, ns) / DistanceSquared(r.o, ps);
00229 } else {
00230 *bsdf = BSDF_ALLOC(tspack, SingleBSDF)(dg, ns,
00231 BSDF_ALLOC(tspack, SkyPortalBxDF)(*this, WorldToLight, dpdu, dpdv, Vector(ns), ps, PortalShapes, ~0U, 0.f));
00232 *pdf = 0.f;
00233 *pdfDirect = 0.f;
00234 for (int i = 0; i < nrPortalShapes; ++i) {
00235 PortalShapes[i]->Sample(.5f, .5f, .5f, &dg);
00236 const Vector w(dg.p - ps);
00237 if (Dot(w, dg.nn) > 0.f) {
00238 const float distance = w.LengthSquared();
00239 *pdf += AbsDot(ns, w) / (sqrtf(distance) * distance);
00240 }
00241 Intersection isect;
00242 RayDifferential ray(r);
00243 ray.mint = -INFINITY;
00244 if (PortalShapes[i]->Intersect(ray, &isect) && Dot(r.d, isect.dg.nn) < 0.f)
00245 *pdfDirect += PortalShapes[i]->Pdf(r.o, isect.dg.p) * DistanceSquared(r.o, isect.dg.p) / DistanceSquared(r.o, ps) * AbsDot(r.d, ns) / AbsDot(r.d, isect.dg.nn);
00246 }
00247 *pdf *= INV_TWOPI / nrPortalShapes;
00248 *pdfDirect /= nrPortalShapes;
00249 }
00250 Vector wh = Normalize(WorldToLight(r.d));
00251 const float phi = SphericalPhi(wh);
00252 const float theta = SphericalTheta(wh);
00253 SWCSpectrum L;
00254 GetSkySpectralRadiance(tspack, theta, phi, &L);
00255 L *= skyScale;
00256 return L;
00257 }
00258
00259 SWCSpectrum SkyLight::Sample_L(const TsPack *tspack, const Point &p,
00260 const Normal &n, float u1, float u2, float u3,
00261 Vector *wi, float *pdf,
00262 VisibilityTester *visibility) const {
00263 if(!havePortalShape) {
00264
00265 float x, y, z;
00266 ConcentricSampleDisk(u1, u2, &x, &y);
00267 z = sqrtf(max(0.f, 1.f - x*x - y*y));
00268 if (u3 < .5) z *= -1;
00269 *wi = Vector(x, y, z);
00270
00271 *pdf = fabsf(wi->z) * INV_TWOPI;
00272
00273 Vector v1, v2;
00274 CoordinateSystem(Normalize(Vector(n)), &v1, &v2);
00275 *wi = Vector(v1.x * wi->x + v2.x * wi->y + n.x * wi->z,
00276 v1.y * wi->x + v2.y * wi->y + n.y * wi->z,
00277 v1.z * wi->x + v2.z * wi->y + n.z * wi->z);
00278 } else {
00279
00280 int shapeIndex = 0;
00281 if(nrPortalShapes > 1) {
00282 shapeIndex = min(nrPortalShapes - 1,
00283 Floor2Int(u3 * nrPortalShapes));
00284 u3 *= nrPortalShapes;
00285 u3 -= shapeIndex;
00286 }
00287 DifferentialGeometry dg;
00288 dg.time = tspack->time;
00289 PortalShapes[shapeIndex]->Sample(p, u1, u2, u3, &dg);
00290 Point ps = dg.p;
00291 *wi = Normalize(ps - p);
00292 if (Dot(*wi, dg.nn) < 0.f)
00293 *pdf = PortalShapes[shapeIndex]->Pdf(p, *wi) / nrPortalShapes;
00294 else {
00295 *pdf = 0.f;
00296 return 0.f;
00297 }
00298 }
00299 visibility->SetRay(p, *wi, tspack->time);
00300 return Le(tspack, RayDifferential(p, *wi));
00301 }
00302 float SkyLight::Pdf(const Point &p, const Normal &n,
00303 const Vector &wi) const {
00304 if (!havePortalShape)
00305 return AbsDot(n, wi) * INV_TWOPI;
00306 else {
00307 float pdf = 0.f;
00308 for (int i = 0; i < nrPortalShapes; ++i) {
00309 Intersection isect;
00310 RayDifferential ray(p, wi);
00311 if (PortalShapes[i]->Intersect(ray, &isect) && Dot(wi, isect.dg.nn) < .0f)
00312 pdf += PortalShapes[i]->Pdf(p, wi);
00313 }
00314 pdf /= nrPortalShapes;
00315 return pdf;
00316 }
00317 }
00318 float SkyLight::Pdf(const Point &p, const Normal &n,
00319 const Point &po, const Normal &ns) const
00320 {
00321 const Vector wi(po - p);
00322 if (!havePortalShape) {
00323 const float d2 = wi.LengthSquared();
00324 return AbsDot(n, wi) * INV_TWOPI * AbsDot(wi, ns) / (d2 * d2);
00325 } else {
00326 float pdf = 0.f;
00327 for (int i = 0; i < nrPortalShapes; ++i) {
00328 Intersection isect;
00329 RayDifferential ray(p, wi);
00330 ray.mint = -INFINITY;
00331 if (PortalShapes[i]->Intersect(ray, &isect) &&
00332 Dot(wi, isect.dg.nn) < 0.f)
00333 pdf += PortalShapes[i]->Pdf(p, isect.dg.p) * DistanceSquared(p, isect.dg.p) / DistanceSquared(p, po) * AbsDot(wi, ns) / AbsDot(wi, isect.dg.nn);
00334 }
00335 pdf /= nrPortalShapes;
00336 return pdf;
00337 }
00338 }
00339 SWCSpectrum SkyLight::Sample_L(const TsPack *tspack, const Point &p,
00340 float u1, float u2, float u3, Vector *wi, float *pdf,
00341 VisibilityTester *visibility) const {
00342 if(!havePortalShape) {
00343 *wi = UniformSampleSphere(u1, u2);
00344 *pdf = UniformSpherePdf();
00345 } else {
00346
00347 int shapeIndex = 0;
00348 if(nrPortalShapes > 1) {
00349 shapeIndex = min(nrPortalShapes - 1,
00350 Floor2Int(u3 * nrPortalShapes));
00351 u3 *= nrPortalShapes;
00352 u3 -= shapeIndex;
00353 }
00354 DifferentialGeometry dg;
00355 dg.time = tspack->time;
00356 PortalShapes[shapeIndex]->Sample(p, u1, u2, u3, &dg);
00357 Point ps = dg.p;
00358 *wi = Normalize(ps - p);
00359 if (Dot(*wi, dg.nn) < 0.f)
00360 *pdf = PortalShapes[shapeIndex]->Pdf(p, *wi) / nrPortalShapes;
00361 else {
00362 *pdf = 0.f;
00363 return 0.f;
00364 }
00365 }
00366 visibility->SetRay(p, *wi, tspack->time);
00367 return Le(tspack, RayDifferential(p, *wi));
00368 }
00369 float SkyLight::Pdf(const Point &, const Vector &) const {
00370 return 1.f / (4.f * M_PI);
00371 }
00372 SWCSpectrum SkyLight::Sample_L(const TsPack *tspack, const Scene *scene,
00373 float u1, float u2, float u3, float u4,
00374 Ray *ray, float *pdf) const {
00375 if (!havePortalShape) {
00376
00377 Point worldCenter;
00378 float worldRadius;
00379 scene->WorldBound().BoundingSphere(&worldCenter, &worldRadius);
00380 worldRadius *= 1.01f;
00381 Point p1 = worldCenter + worldRadius *
00382 UniformSampleSphere(u1, u2);
00383 Point p2 = worldCenter + worldRadius *
00384 UniformSampleSphere(u3, u4);
00385
00386 ray->o = p1;
00387 ray->d = Normalize(p2-p1);
00388
00389 Vector to_center = Normalize(worldCenter - p1);
00390 float costheta = AbsDot(to_center,ray->d);
00391 *pdf = costheta / ((4.f * M_PI * worldRadius * worldRadius));
00392 } else {
00393
00394
00395 int shapeidx = 0;
00396 if(nrPortalShapes > 1)
00397 shapeidx = min(nrPortalShapes - 1,
00398 Floor2Int(tspack->rng->floatValue() * nrPortalShapes));
00399
00400 DifferentialGeometry dg;
00401 dg.time = tspack->time;
00402 PortalShapes[shapeidx]->Sample(u1, u2, tspack->rng->floatValue(), &dg);
00403 ray->o = dg.p;
00404 ray->d = UniformSampleSphere(u3, u4);
00405 if (Dot(ray->d, dg.nn) < 0.) ray->d *= -1;
00406
00407 *pdf = PortalShapes[shapeidx]->Pdf(ray->o) * INV_TWOPI / nrPortalShapes;
00408 }
00409
00410 return Le(tspack, RayDifferential(ray->o, -ray->d));
00411 }
00412 bool SkyLight::Sample_L(const TsPack *tspack, const Scene *scene, float u1, float u2, float u3, BSDF **bsdf, float *pdf, SWCSpectrum *Le) const
00413 {
00414 Point worldCenter;
00415 float worldRadius;
00416 scene->WorldBound().BoundingSphere(&worldCenter, &worldRadius);
00417 if (!havePortalShape) {
00418 Point ps = worldCenter + worldRadius * UniformSampleSphere(u1, u2);
00419 Normal ns = Normal(Normalize(worldCenter - ps));
00420 Vector dpdu, dpdv;
00421 CoordinateSystem(Vector(ns), &dpdu, &dpdv);
00422 DifferentialGeometry dg(ps, ns, dpdu, dpdv, Normal(0, 0, 0), Normal (0, 0, 0), 0, 0, NULL);
00423 *bsdf = BSDF_ALLOC(tspack, SingleBSDF)(dg, ns,
00424 BSDF_ALLOC(tspack, SkyBxDF)(*this, WorldToLight, dpdu, dpdv, Vector(ns)));
00425 *pdf = 1.f / (4.f * M_PI * worldRadius * worldRadius);
00426 } else {
00427
00428 int shapeIndex = 0;
00429 if(nrPortalShapes > 1) {
00430 u3 *= nrPortalShapes;
00431 shapeIndex = min(nrPortalShapes - 1, Floor2Int(u3));
00432 u3 -= shapeIndex;
00433 }
00434 DifferentialGeometry dgs;
00435 dgs.time = tspack->time;
00436 PortalShapes[shapeIndex]->Sample(.5f, .5f, .5f, &dgs);
00437 Vector wi(UniformSampleHemisphere(u1, u2));
00438 wi = Normalize(wi.x * Normalize(dgs.dpdu) + wi.y * Normalize(dgs.dpdv) - wi.z * Vector(dgs.nn));
00439 Vector toCenter(worldCenter - dgs.p);
00440 float centerDistance = Dot(toCenter, toCenter);
00441 float approach = Dot(toCenter, wi);
00442 float distance = approach + sqrtf(worldRadius * worldRadius - centerDistance + approach * approach);
00443 Point ps(dgs.p + distance * wi);
00444 Normal ns(Normalize(worldCenter - ps));
00445 Vector dpdu, dpdv;
00446 CoordinateSystem(Vector(ns), &dpdu, &dpdv);
00447 DifferentialGeometry dg(ps, ns, dpdu, dpdv, Normal(0, 0, 0), Normal (0, 0, 0), 0, 0, NULL);
00448 *bsdf = BSDF_ALLOC(tspack, SingleBSDF)(dg, ns,
00449 BSDF_ALLOC(tspack, SkyPortalBxDF)(*this, WorldToLight, dpdu, dpdv, Vector(ns), ps, PortalShapes, shapeIndex, u3));
00450 *pdf = AbsDot(ns, wi) / (distance * distance);
00451 for (int i = 0; i < nrPortalShapes; ++i) {
00452 if (i == shapeIndex)
00453 continue;
00454 PortalShapes[i]->Sample(.5f, .5f, .5f, &dgs);
00455 wi = ps - dgs.p;
00456 if (Dot(wi, dg.nn) <= 0.f) {
00457 distance = wi.LengthSquared();
00458 *pdf += AbsDot(ns, wi) / (sqrtf(distance) * distance);
00459 }
00460 }
00461 *pdf *= INV_TWOPI / nrPortalShapes;
00462 }
00463 *Le = SWCSpectrum(skyScale);
00464 return true;
00465 }
00466 bool SkyLight::Sample_L(const TsPack *tspack, const Scene *scene, const Point &p, const Normal &n,
00467 float u1, float u2, float u3, BSDF **bsdf, float *pdf, float *pdfDirect,
00468 VisibilityTester *visibility, SWCSpectrum *Le) const
00469 {
00470 Vector wi;
00471 int shapeIndex = 0;
00472 Point worldCenter;
00473 float worldRadius;
00474 scene->WorldBound().BoundingSphere(&worldCenter, &worldRadius);
00475 if(!havePortalShape) {
00476
00477 float x, y, z;
00478 ConcentricSampleDisk(u1, u2, &x, &y);
00479 z = sqrtf(max(0.f, 1.f - x*x - y*y));
00480 if (u3 < .5f)
00481 z = -z;
00482 wi = Vector(x, y, z);
00483
00484 *pdfDirect = fabsf(wi.z) * INV_TWOPI;
00485
00486 Vector v1, v2;
00487 CoordinateSystem(Normalize(Vector(n)), &v1, &v2);
00488 wi = Vector(v1 * wi.x + v2 * wi.y + Vector(n) * wi.z);
00489 } else {
00490
00491 if (nrPortalShapes > 1) {
00492 u3 *= nrPortalShapes;
00493 shapeIndex = min(nrPortalShapes - 1, Floor2Int(u3));
00494 u3 -= shapeIndex;
00495 }
00496 DifferentialGeometry dg;
00497 dg.time = tspack->time;
00498 PortalShapes[shapeIndex]->Sample(p, u1, u2, u3, &dg);
00499 Point ps = dg.p;
00500 wi = Normalize(ps - p);
00501 if (Dot(wi, dg.nn) < 0.f) {
00502 *pdfDirect = PortalShapes[shapeIndex]->Pdf(p, ps);
00503 *pdfDirect *= DistanceSquared(p, dg.p) / AbsDot(wi, dg.nn);
00504 } else {
00505 *Le = 0.f;
00506 return false;
00507 }
00508 }
00509 Vector toCenter(worldCenter - p);
00510 float centerDistance = Dot(toCenter, toCenter);
00511 float approach = Dot(toCenter, wi);
00512 float distance = approach + sqrtf(worldRadius * worldRadius - centerDistance + approach * approach);
00513 Point ps(p + distance * wi);
00514 Normal ns(Normalize(worldCenter - ps));
00515 Vector dpdu, dpdv;
00516 CoordinateSystem(Vector(ns), &dpdu, &dpdv);
00517 DifferentialGeometry dg(ps, ns, dpdu, dpdv, Normal(0, 0, 0), Normal (0, 0, 0), 0, 0, NULL);
00518 if (!havePortalShape) {
00519 *bsdf = BSDF_ALLOC(tspack, SingleBSDF)(dg, ns,
00520 BSDF_ALLOC(tspack, SkyBxDF)(*this, WorldToLight, dpdu, dpdv, Vector(ns)));
00521 *pdf = 1.f / (4.f * M_PI * worldRadius * worldRadius);
00522 } else {
00523 *bsdf = BSDF_ALLOC(tspack, SingleBSDF)(dg, ns,
00524 BSDF_ALLOC(tspack, SkyPortalBxDF)(*this, WorldToLight, dpdu, dpdv, Vector(ns), ps, PortalShapes, shapeIndex, u3));
00525 *pdf = 0.f;
00526 DifferentialGeometry dgs;
00527 dgs.time = tspack->time;
00528 for (int i = 0; i < nrPortalShapes; ++i) {
00529 PortalShapes[i]->Sample(.5f, .5f, .5f, &dgs);
00530 Vector w(ps - dgs.p);
00531 if (Dot(wi, dg.nn) < 0.f) {
00532 float distance = w.LengthSquared();
00533 *pdf += AbsDot(ns, w) / (sqrtf(distance) * distance);
00534 }
00535 }
00536 *pdf *= INV_TWOPI / nrPortalShapes;
00537 }
00538 *pdfDirect *= AbsDot(wi, ns) / (distance * distance);
00539 if (havePortalShape) {
00540 for (int i = 0; i < nrPortalShapes; ++i) {
00541 if (i == shapeIndex)
00542 continue;
00543 Intersection isect;
00544 RayDifferential ray(p, wi);
00545 ray.mint = -INFINITY;
00546 if (PortalShapes[i]->Intersect(ray, &isect) &&
00547 Dot(wi, isect.dg.nn) < 0.f)
00548 *pdfDirect += PortalShapes[i]->Pdf(p, isect.dg.p) * DistanceSquared(p, isect.dg.p) / DistanceSquared(p, ps) * AbsDot(wi, ns) / AbsDot(wi, isect.dg.nn);
00549 }
00550 *pdfDirect /= nrPortalShapes;
00551 }
00552 visibility->SetSegment(p, ps, tspack->time);
00553 *Le = SWCSpectrum(skyScale);
00554 return true;
00555 }
00556
00557 Light* SkyLight::CreateLight(const Transform &light2world,
00558 const ParamSet ¶mSet, const TextureParams &tp) {
00559 float scale = paramSet.FindOneFloat("gain", 1.f);
00560 int nSamples = paramSet.FindOneInt("nsamples", 1);
00561 Vector sundir = paramSet.FindOneVector("sundir", Vector(0,0,1));
00562 Normalize(sundir);
00563 float turb = paramSet.FindOneFloat("turbidity", 2.0f);
00564 float aconst = paramSet.FindOneFloat("aconst", 1.0f);
00565 float bconst = paramSet.FindOneFloat("bconst", 1.0f);
00566 float cconst = paramSet.FindOneFloat("cconst", 1.0f);
00567 float dconst = paramSet.FindOneFloat("dconst", 1.0f);
00568 float econst = paramSet.FindOneFloat("econst", 1.0f);
00569
00570 return new SkyLight(light2world, scale, nSamples, sundir, turb, aconst, bconst, cconst, dconst, econst);
00571 }
00572
00573
00574 inline float RiAngleBetween(float thetav, float phiv, float theta, float phi)
00575 {
00576 float cospsi = sin(thetav) * sin(theta) * cos(phi-phiv) + cos(thetav) * cos(theta);
00577 if (cospsi > 1) return 0;
00578 if (cospsi < -1) return M_PI;
00579 return acos(cospsi);
00580 }
00581
00582
00583
00584
00585
00586
00587
00588
00589 void SkyLight::InitSunThetaPhi()
00590 {
00591 Vector wh = Normalize(sundir);
00592 phiS = SphericalPhi(wh);
00593 thetaS = SphericalTheta(wh);
00594 }
00595
00596
00597
00598
00599
00600
00601
00602
00603 void SkyLight::GetSkySpectralRadiance(const TsPack *tspack, const float theta, const float phi, SWCSpectrum * const dst_spect) const
00604 {
00605
00606 const float theta_fin = min(theta,(M_PI * 0.5f) - 0.001f);
00607 const float gamma = RiAngleBetween(theta,phi,thetaS,phiS);
00608
00609
00610 const float x = zenith_x * PerezBase(perez_x, theta_fin, gamma);
00611 const float y = zenith_y * PerezBase(perez_y, theta_fin, gamma);
00612 const float Y = zenith_Y * PerezBase(perez_Y, theta_fin, gamma);
00613
00614 ChromaticityToSpectrum(tspack, x, y, dst_spect);
00615 *dst_spect *= Y;
00616 }
00617
00618
00619 void SkyLight::ChromaticityToSpectrum(const TsPack *tspack, const float x, const float y, SWCSpectrum * const dst_spect) const
00620 {
00621 const float den = 1.0f / (0.0241f + 0.2562f * x - 0.7341f * y);
00622 const float M1 = (-1.3515f - 1.7703f * x + 5.9114f * y) * den;
00623 const float M2 = (0.03f - 31.4424f * x + 30.0717f * y) * den;
00624
00625 for (unsigned int j = 0; j < WAVELENGTH_SAMPLES; ++j) {
00626 const float w = tspack->swl->w[j];
00627 dst_spect->c[j] = S0.sample(w) + M1 * S1.sample(w) +
00628 M2 * S2.sample(w);
00629 }
00630 *dst_spect /= S0Y + M1 * S1Y + M2 * S2Y;
00631 }
00632
00633 static DynamicLoader::RegisterLight<SkyLight> r("sky");