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 "sun.h"
00025 #include "spd.h"
00026 #include "regular.h"
00027 #include "irregular.h"
00028 #include "mc.h"
00029 #include "paramset.h"
00030 #include "reflection/bxdf.h"
00031 #include "dynload.h"
00032
00033 #include "data/sun_spect.h"
00034
00035 using namespace lux;
00036
00037 class SunBxDF : public BxDF
00038 {
00039 public:
00040 SunBxDF(float sin2Max, float radius) : BxDF(BxDFType(BSDF_REFLECTION | BSDF_DIFFUSE)), sin2ThetaMax(sin2Max), cosThetaMax(sqrtf(1.f - sin2Max)), worldRadius(radius) {}
00041 virtual ~SunBxDF() { }
00042 virtual void f(const TsPack *tspack, const Vector &wo, const Vector &wi, SWCSpectrum *const f) const {
00043 if (wo.z < 0.f || wi.z < 0.f || (wo.x * wo.x + wo.y * wo.y) > sin2ThetaMax || (wi.x * wi.x + wi.y * wi.y) > sin2ThetaMax)
00044 return;
00045 *f += SWCSpectrum(1.f);
00046 }
00047 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
00048 {
00049 *wi = UniformSampleCone(u1, u2, cosThetaMax);
00050 *pdf = UniformConePdf(cosThetaMax);
00051 if (pdfBack)
00052 *pdfBack = Pdf(tspack, *wi, wo);
00053 *f = SWCSpectrum(1.f);
00054 return true;
00055 }
00056 virtual float Pdf(const TsPack *tspack, const Vector &wi, const Vector &wo) const
00057 {
00058 if (wo.z < 0.f || wi.z < 0.f || (wo.x * wo.x + wo.y * wo.y) > sin2ThetaMax || (wi.x * wi.x + wi.y * wi.y) > sin2ThetaMax)
00059 return 0.f;
00060 else
00061 return UniformConePdf(cosThetaMax);
00062 }
00063 private:
00064 float sin2ThetaMax, cosThetaMax, worldRadius;
00065 };
00066
00067
00068 SunLight::SunLight(const Transform &light2world,
00069 const float sunscale, const Vector &dir, float turb , float relSize, int ns)
00070 : Light(light2world, ns) {
00071 sundir = Normalize(LightToWorld(dir));
00072 turbidity = turb;
00073
00074 CoordinateSystem(sundir, &x, &y);
00075
00076
00077
00078 const float sunRadius = 695500.f;
00079 const float sunMeanDistance = 149600000.f;
00080 if (relSize * sunRadius <= sunMeanDistance) {
00081 sin2ThetaMax = relSize * sunRadius / sunMeanDistance;
00082 sin2ThetaMax *= sin2ThetaMax;
00083 cosThetaMax = sqrtf(1.f - sin2ThetaMax);
00084 } else {
00085 std::stringstream ss;
00086 ss << "Reducing relative sun size to " << sunMeanDistance / sunRadius;
00087 luxError(LUX_LIMIT, LUX_WARNING, ss.str().c_str());
00088 cosThetaMax = 0.f;
00089 sin2ThetaMax = 1.f;
00090 }
00091
00092 Vector wh = Normalize(sundir);
00093 phiS = SphericalPhi(wh);
00094 thetaS = SphericalTheta(wh);
00095
00096
00097 IrregularSPD k_oCurve(sun_k_oWavelengths, sun_k_oAmplitudes, 64);
00098 IrregularSPD k_gCurve(sun_k_gWavelengths, sun_k_gAmplitudes, 4);
00099 IrregularSPD k_waCurve(sun_k_waWavelengths, sun_k_waAmplitudes, 13);
00100
00101 RegularSPD solCurve(sun_solAmplitudes, 380, 750, 38);
00102
00103 float beta = 0.04608365822050f * turbidity - 0.04586025928522f;
00104 float tauR, tauA, tauO, tauG, tauWA;
00105
00106 float m = 1.f / (cos(thetaS) + 0.00094f * powf(1.6386f - thetaS, -1.253f));
00107
00108 int i;
00109 float lambda;
00110
00111 float Ldata[91];
00112 for(i = 0, lambda = 350.f; i < 91; ++i, lambda += 5.f) {
00113
00114 tauR = expf( -m * 0.008735f * powf(lambda / 1000.f, -4.08f));
00115
00116
00117
00118 const float alpha = 1.3f;
00119 tauA = expf(-m * beta * pow(lambda / 1000.f, -alpha));
00120
00121
00122 const float lOzone = .35f;
00123 tauO = expf(-m * k_oCurve.sample(lambda) * lOzone);
00124
00125 tauG = expf(-1.41f * k_gCurve.sample(lambda) * m / powf(1.f + 118.93f * k_gCurve.sample(lambda) * m, 0.45f));
00126
00127
00128 const float w = 2.0f;
00129 tauWA = expf(-0.2385f * k_waCurve.sample(lambda) * w * m /
00130 powf(1.f + 20.07f * k_waCurve.sample(lambda) * w * m, 0.45f));
00131
00132 Ldata[i] = (solCurve.sample(lambda) * tauR * tauA * tauO * tauG * tauWA);
00133 }
00134 LSPD = new RegularSPD(Ldata, 350,800,91);
00135 LSPD->Scale(sunscale);
00136 }
00137
00138 SWCSpectrum SunLight::Le(const TsPack *tspack, const RayDifferential &r) const {
00139 Vector w = r.d;
00140 if(cosThetaMax < 1.f && Dot(w,sundir) > cosThetaMax)
00141 return SWCSpectrum(tspack, LSPD);
00142 else
00143 return SWCSpectrum(0.f);
00144 }
00145
00146 SWCSpectrum SunLight::Le(const TsPack *tspack, const Scene *scene, const Ray &r,
00147 const Normal &n, BSDF **bsdf, float *pdf, float *pdfDirect) const
00148 {
00149 const float xD = Dot(r.d, x);
00150 const float yD = Dot(r.d, y);
00151 if (cosThetaMax == 1.f || Dot(r.d, sundir) < 0.f || (xD * xD + yD * yD) > sin2ThetaMax) {
00152 *bsdf = NULL;
00153 *pdf = 0.f;
00154 *pdfDirect = 0.f;
00155 return SWCSpectrum(0.f);
00156 }
00157 Point worldCenter;
00158 float worldRadius;
00159 scene->WorldBound().BoundingSphere(&worldCenter, &worldRadius);
00160 Vector toCenter(worldCenter - r.o);
00161 float approach = Dot(toCenter, sundir);
00162 float distance = (approach + worldRadius) / Dot(r.d, sundir);
00163 Point ps(r.o + distance * r.d);
00164 Normal ns(-sundir);
00165 DifferentialGeometry dg(ps, ns, -x, y, Normal(0, 0, 0), Normal (0, 0, 0), 0, 0, NULL);
00166 *bsdf = BSDF_ALLOC(tspack, SingleBSDF)(dg, ns,
00167 BSDF_ALLOC(tspack, SunBxDF)(sin2ThetaMax, worldRadius));
00168 if (!havePortalShape)
00169 *pdf = 1.f / (M_PI * worldRadius * worldRadius);
00170 else {
00171 *pdf = 0.f;
00172 for (int i = 0; i < nrPortalShapes; ++i) {
00173 Intersection isect;
00174 RayDifferential ray(ps, sundir);
00175 ray.mint = -INFINITY;
00176 if (PortalShapes[i]->Intersect(ray, &isect)) {
00177 float cosPortal = Dot(-sundir, isect.dg.nn);
00178 if (cosPortal > 0.f)
00179 *pdf += PortalShapes[i]->Pdf(isect.dg.p) / cosPortal;
00180 }
00181 }
00182 *pdf /= nrPortalShapes;
00183 }
00184 *pdfDirect = UniformConePdf(cosThetaMax) * AbsDot(r.d, ns) / DistanceSquared(r.o, ps);
00185 return SWCSpectrum(tspack, LSPD);
00186 }
00187
00188 bool SunLight::checkPortals(Ray portalRay) const {
00189 if (!havePortalShape)
00190 return true;
00191
00192 Ray isectRay(portalRay);
00193 Intersection isect;
00194 bool found = false;
00195 for (int i = 0; i < nrPortalShapes; ++i) {
00196
00197
00198 if (PortalShapes[i]->Intersect(isectRay, &isect)) {
00199
00200 if (Dot(portalRay.d, isect.dg.nn) < 0.f) {
00201 found = true;
00202 break;
00203 }
00204 }
00205 }
00206
00207 return found;
00208 }
00209
00210 SWCSpectrum SunLight::Sample_L(const TsPack *tspack, const Point &p, float u1, float u2, float u3,
00211 Vector *wi, float *pdf, VisibilityTester *visibility) const {
00212 if(cosThetaMax == 1.f) {
00213 *wi = sundir;
00214 *pdf = 1.f;
00215 } else {
00216 *wi = UniformSampleCone(u1, u2, cosThetaMax, x, y, sundir);
00217 *pdf = UniformConePdf(cosThetaMax);
00218 }
00219 visibility->SetRay(p, *wi, tspack->time);
00220
00221
00222
00223
00224
00225 return SWCSpectrum(tspack, LSPD);
00226 }
00227 float SunLight::Pdf(const Point &, const Vector &) const {
00228 if(cosThetaMax == 1)
00229 return 0.f;
00230 else
00231 return UniformConePdf(cosThetaMax);
00232 }
00233 float SunLight::Pdf(const Point &p, const Normal &n,
00234 const Point &po, const Normal &ns) const
00235 {
00236 const float cosTheta = AbsDot(Normalize(p - po), ns);
00237 if(cosTheta < cosThetaMax)
00238 return 0.f;
00239 else
00240 return UniformConePdf(cosThetaMax) * cosTheta / DistanceSquared(p, po);
00241 }
00242
00243 SWCSpectrum SunLight::Sample_L(const TsPack *tspack, const Scene *scene,
00244 float u1, float u2, float u3, float u4,
00245 Ray *ray, float *pdf) const {
00246 if (!havePortalShape) {
00247
00248 Point worldCenter;
00249 float worldRadius;
00250 scene->WorldBound().BoundingSphere(&worldCenter, &worldRadius);
00251 float d1, d2;
00252 ConcentricSampleDisk(u1, u2, &d1, &d2);
00253 Point Pdisk =
00254 worldCenter + worldRadius * (d1 * x + d2 * y);
00255
00256 ray->o = Pdisk + worldRadius * sundir;
00257 ray->d = -UniformSampleCone(u3, u4, cosThetaMax, x, y, sundir);
00258 *pdf = UniformConePdf(cosThetaMax) / (M_PI * worldRadius * worldRadius);
00259
00260 return SWCSpectrum(tspack, LSPD);
00261 } else {
00262
00263
00264 int shapeidx = 0;
00265 if(nrPortalShapes > 1)
00266 shapeidx = min<float>(nrPortalShapes - 1,
00267 Floor2Int(tspack->rng->floatValue() * nrPortalShapes));
00268
00269 DifferentialGeometry dg;
00270 dg.time = tspack->time;
00271 PortalShapes[shapeidx]->Sample(u1, u2, tspack->rng->floatValue(), &dg);
00272 ray->o = dg.p;
00273 float pdfPortal = PortalShapes[shapeidx]->Pdf(ray->o) / nrPortalShapes;
00274
00275 ray->d = -UniformSampleCone(u3, u4, cosThetaMax, x, y, sundir);
00276 float pdfSun = UniformConePdf(cosThetaMax);
00277
00278
00279 Point worldCenter;
00280 float worldRadius;
00281 scene->WorldBound().BoundingSphere(&worldCenter, &worldRadius);
00282 Vector centerSunDisk = Vector(worldCenter) + worldRadius;
00283
00284
00285
00286 float PnRd = Dot(sundir, ray->d);
00287 if (PnRd == 0.f)
00288 return SWCSpectrum(0.f);
00289
00290 float distanceSunDisk = centerSunDisk.Length();
00291 float t = - (Dot(sundir, Vector(ray->o)) + distanceSunDisk) / PnRd;
00292
00293
00294 ray->o = ray->o - ray->d * t;
00295
00296 *pdf = pdfSun * pdfPortal;
00297
00298 if (Dot(ray->d, dg.nn) <= 0.f)
00299 return SWCSpectrum(0.f);
00300 else
00301 return SWCSpectrum(tspack, LSPD);
00302 }
00303 }
00304
00305 bool SunLight::Sample_L(const TsPack *tspack, const Scene *scene, float u1, float u2, float u3, BSDF **bsdf, float *pdf, SWCSpectrum *Le) const
00306 {
00307 Point worldCenter;
00308 float worldRadius;
00309 scene->WorldBound().BoundingSphere(&worldCenter, &worldRadius);
00310
00311 Point ps;
00312 Normal ns(-sundir);
00313 if (!havePortalShape) {
00314 float d1, d2;
00315 ConcentricSampleDisk(u1, u2, &d1, &d2);
00316 ps = worldCenter + worldRadius * (sundir + d1 * x + d2 * y);
00317 *pdf = 1.f / (M_PI * worldRadius * worldRadius);
00318 } else {
00319
00320 int shapeIndex = 0;
00321 if (nrPortalShapes > 1) {
00322 u3 *= nrPortalShapes;
00323 shapeIndex = min(nrPortalShapes - 1, Floor2Int(u3));
00324 u3 -= shapeIndex;
00325 }
00326
00327 DifferentialGeometry dg;
00328 dg.time = tspack->time;
00329 PortalShapes[shapeIndex]->Sample(u1, u2, u3, &dg);
00330 ps = dg.p;
00331 const float cosPortal = Dot(ns, dg.nn);
00332 if (cosPortal <= 0.f) {
00333 *Le = SWCSpectrum(0.f);
00334 return false;
00335 }
00336
00337 *pdf = PortalShapes[shapeIndex]->Pdf(ps) / cosPortal;
00338 for (int i = 0; i < nrPortalShapes; ++i) {
00339 if (i == shapeIndex)
00340 continue;
00341 Intersection isect;
00342 RayDifferential ray(ps, sundir);
00343 ray.mint = -INFINITY;
00344 if (PortalShapes[i]->Intersect(ray, &isect)) {
00345 float cosP = Dot(ns, isect.dg.nn);
00346 if (cosP > 0.f)
00347 *pdf += PortalShapes[i]->Pdf(isect.dg.p) / cosP;
00348 }
00349 }
00350 *pdf /= nrPortalShapes;
00351 if (!(*pdf > 0.f)) {
00352 *Le = SWCSpectrum(0.f);
00353 return false;
00354 }
00355
00356 ps += (worldRadius + Dot(worldCenter - ps, sundir)) * sundir;
00357 }
00358
00359 DifferentialGeometry dg(ps, ns, -x, y, Normal(0, 0, 0), Normal(0, 0, 0), 0, 0, NULL);
00360 *bsdf = BSDF_ALLOC(tspack, SingleBSDF)(dg, ns,
00361 BSDF_ALLOC(tspack, SunBxDF)(sin2ThetaMax, worldRadius));
00362
00363 *Le = SWCSpectrum(tspack, LSPD);
00364 return true;
00365 }
00366
00367 bool SunLight::Sample_L(const TsPack *tspack, const Scene *scene, const Point &p, const Normal &n,
00368 float u1, float u2, float u3, BSDF **bsdf, float *pdf, float *pdfDirect,
00369 VisibilityTester *visibility, SWCSpectrum *Le) const
00370 {
00371 Vector wi;
00372 if(cosThetaMax == 1.f) {
00373 wi = sundir;
00374 *pdfDirect = 1.f;
00375
00376
00377
00378
00379
00380
00381
00382 } else {
00383 wi = UniformSampleCone(u1, u2, cosThetaMax, x, y, sundir);
00384 *pdfDirect = UniformConePdf(cosThetaMax);
00385
00386
00387
00388
00389
00390
00391
00392 }
00393
00394 Point worldCenter;
00395 float worldRadius;
00396 scene->WorldBound().BoundingSphere(&worldCenter, &worldRadius);
00397 Vector toCenter(worldCenter - p);
00398 float approach = Dot(toCenter, sundir);
00399 float distance = (approach + worldRadius) / Dot(wi, sundir);
00400 Point ps(p + distance * wi);
00401 Normal ns(-sundir);
00402
00403 DifferentialGeometry dg(ps, ns, -x, y, Normal(0, 0, 0), Normal (0, 0, 0), 0, 0, NULL);
00404 *bsdf = BSDF_ALLOC(tspack, SingleBSDF)(dg, ns,
00405 BSDF_ALLOC(tspack, SunBxDF)(sin2ThetaMax, worldRadius));
00406 if (!havePortalShape)
00407 *pdf = 1.f / (M_PI * worldRadius * worldRadius);
00408 else {
00409 *pdf = 0.f;
00410 for (int i = 0; i < nrPortalShapes; ++i) {
00411 Intersection isect;
00412 RayDifferential ray(ps, sundir);
00413 ray.mint = -INFINITY;
00414 if (PortalShapes[i]->Intersect(ray, &isect)) {
00415 float cosPortal = Dot(ns, isect.dg.nn);
00416 if (cosPortal > 0.f)
00417 *pdf += PortalShapes[i]->Pdf(isect.dg.p) / cosPortal;
00418 }
00419 }
00420 *pdf /= nrPortalShapes;
00421 }
00422 if (cosThetaMax < 1.f)
00423 *pdfDirect *= AbsDot(wi, ns) / DistanceSquared(p, ps);
00424 visibility->SetSegment(p, ps, tspack->time);
00425
00426 *Le = SWCSpectrum(tspack, LSPD);
00427 return true;
00428 }
00429
00430 Light* SunLight::CreateLight(const Transform &light2world,
00431 const ParamSet ¶mSet, const TextureParams &tp) {
00432
00433 float scale = paramSet.FindOneFloat("gain", 1.f);
00434 int nSamples = paramSet.FindOneInt("nsamples", 1);
00435 Vector sundir = paramSet.FindOneVector("sundir", Vector(0,0,-1));
00436 float turb = paramSet.FindOneFloat("turbidity", 2.0f);
00437 float relSize = paramSet.FindOneFloat("relsize", 1.0f);
00438 return new SunLight(light2world, scale, sundir, turb, relSize, nSamples);
00439 }
00440
00441 static DynamicLoader::RegisterLight<SunLight> r("sun");