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 "reflection/bxdf.h"
00029
00030 #include "data/skychroma_spect.h"
00031
00032 using namespace lux;
00033
00034 class SkyBxDF : public BxDF
00035 {
00036 public:
00037 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) {}
00038 SWCSpectrum f(const Vector &wo, const Vector &wi) const
00039 {
00040 Vector w(wi.x * X.x + wi.y * Y.x + wi.z * Z.x, wi.x * X.y + wi.y * Y.y + wi.z * Z.y, wi.x * X.z + wi.y * Y.z + wi.z * Z.z);
00041 w = -Normalize(WorldToLight(w));
00042 const float phi = SphericalPhi(w);
00043 const float theta = SphericalTheta(w);
00044 SWCSpectrum L;
00045 skyLight.GetSkySpectralRadiance(theta, phi, &L);
00046 return L;
00047 }
00048 private:
00049 const SkyLight &skyLight;
00050 const Transform &WorldToLight;
00051 Vector X, Y, Z;
00052 };
00053
00054
00055 SkyLight::~SkyLight() {
00056 }
00057
00058 SkyLight::SkyLight(const Transform &light2world,
00059 const float skyscale, int ns, Vector sd, float turb,
00060 float aconst, float bconst, float cconst, float dconst, float econst)
00061 : Light(light2world, ns) {
00062 skyScale = skyscale;
00063 sundir = sd;
00064 turbidity = turb;
00065
00066 InitSunThetaPhi();
00067
00068 float theta2 = thetaS*thetaS;
00069 float theta3 = theta2*thetaS;
00070 float T = turb;
00071 float T2 = turb*turb;
00072
00073 float chi = (4.0/9.0 - T / 120.0) * (M_PI - 2 * thetaS);
00074 zenith_Y = (4.0453 * T - 4.9710) * tan(chi) - .2155 * T + 2.4192;
00075 zenith_Y *= 1000;
00076
00077 zenith_x =
00078 (+0.00165*theta3 - 0.00374*theta2 + 0.00208*thetaS + 0) * T2 +
00079 (-0.02902*theta3 + 0.06377*theta2 - 0.03202*thetaS + 0.00394) * T +
00080 (+0.11693*theta3 - 0.21196*theta2 + 0.06052*thetaS + 0.25885);
00081
00082 zenith_y =
00083 (+0.00275*theta3 - 0.00610*theta2 + 0.00316*thetaS + 0) * T2 +
00084 (-0.04214*theta3 + 0.08970*theta2 - 0.04153*thetaS + 0.00515) * T +
00085 (+0.15346*theta3 - 0.26756*theta2 + 0.06669*thetaS + 0.26688);
00086
00087 perez_Y[1] = (0.17872 *T - 1.46303)*aconst;
00088 perez_Y[2] = (-0.35540 *T + 0.42749)*bconst;
00089 perez_Y[3] = (-0.02266 *T + 5.32505)*cconst;
00090 perez_Y[4] = (0.12064 *T - 2.57705)*dconst;
00091 perez_Y[5] = (-0.06696 *T + 0.37027)*econst;
00092
00093 perez_x[1] = (-0.01925 *T - 0.25922)*aconst;
00094 perez_x[2] = (-0.06651 *T + 0.00081)*bconst;
00095 perez_x[3] = (-0.00041 *T + 0.21247)*cconst;
00096 perez_x[4] = (-0.06409 *T - 0.89887)*dconst;
00097 perez_x[5] = (-0.00325 *T + 0.04517)*econst;
00098
00099 perez_y[1] = (-0.01669 *T - 0.26078)*aconst;
00100 perez_y[2] = (-0.09495 *T + 0.00921)*bconst;
00101 perez_y[3] = (-0.00792 *T + 0.21023)*cconst;
00102 perez_y[4] = (-0.04405 *T - 1.65369)*dconst;
00103 perez_y[5] = (-0.01092 *T + 0.05291)*econst;
00104 }
00105
00106 SWCSpectrum SkyLight::Le(const RayDifferential &r) const {
00107 Vector w = r.d;
00108
00109
00110 Vector wh = Normalize(WorldToLight(w));
00111 const float phi = SphericalPhi(wh);
00112 const float theta = SphericalTheta(wh);
00113
00114 SWCSpectrum L;
00115 GetSkySpectralRadiance(theta,phi,(SWCSpectrum * const)&L);
00116 L *= skyScale;
00117
00118 return L;
00119 }
00120 SWCSpectrum SkyLight::Le(const Scene *scene, const Ray &r,
00121 const Normal &n, BSDF **bsdf, float *pdf, float *pdfDirect) const
00122 {
00123 Point worldCenter;
00124 float worldRadius;
00125 scene->WorldBound().BoundingSphere(&worldCenter, &worldRadius);
00126 worldRadius *= 1.01f;
00127 Vector toCenter(worldCenter - r.o);
00128 float centerDistance = Dot(toCenter, toCenter);
00129 float approach = Dot(toCenter, r.d);
00130 float distance = approach + sqrtf(worldRadius * worldRadius - centerDistance + approach * approach);
00131 Point ps(r.o + distance * r.d);
00132 Normal ns(Normalize(worldCenter - ps));
00133 Vector dpdu, dpdv;
00134 CoordinateSystem(Vector(ns), &dpdu, &dpdv);
00135 DifferentialGeometry dg(ps, ns, dpdu, dpdv, Vector(0, 0, 0), Vector (0, 0, 0), 0, 0, NULL);
00136 *bsdf = BSDF_ALLOC(BSDF)(dg, ns);
00137 (*bsdf)->Add(BSDF_ALLOC(SkyBxDF)(*this, WorldToLight, dpdu, dpdv, Vector(ns)));
00138 *pdf = 1.f / (4.f * M_PI * worldRadius * worldRadius);
00139 *pdfDirect = AbsDot(r.d, n) * INV_TWOPI * AbsDot(r.d, ns) / DistanceSquared(r.o, ps);
00140 Vector wh = Normalize(WorldToLight(r.d));
00141 const float phi = SphericalPhi(wh);
00142 const float theta = SphericalTheta(wh);
00143 SWCSpectrum L;
00144 GetSkySpectralRadiance(theta, phi, &L);
00145 L *= skyScale;
00146 return L;
00147 }
00148
00149 SWCSpectrum SkyLight::Sample_L(const Point &p,
00150 const Normal &n, float u1, float u2, float u3,
00151 Vector *wi, float *pdf,
00152 VisibilityTester *visibility) const {
00153 if(!havePortalShape) {
00154
00155 float x, y, z;
00156 ConcentricSampleDisk(u1, u2, &x, &y);
00157 z = sqrtf(max(0.f, 1.f - x*x - y*y));
00158 if (u3 < .5) z *= -1;
00159 *wi = Vector(x, y, z);
00160
00161 *pdf = fabsf(wi->z) * INV_TWOPI;
00162
00163 Vector v1, v2;
00164 CoordinateSystem(Normalize(Vector(n)), &v1, &v2);
00165 *wi = Vector(v1.x * wi->x + v2.x * wi->y + n.x * wi->z,
00166 v1.y * wi->x + v2.y * wi->y + n.y * wi->z,
00167 v1.z * wi->x + v2.z * wi->y + n.z * wi->z);
00168 } else {
00169
00170 int shapeidx = 0;
00171 if(nrPortalShapes > 1)
00172 shapeidx = min<float>(nrPortalShapes - 1,
00173 Floor2Int(u3 * nrPortalShapes));
00174 Normal ns;
00175 Point ps;
00176 bool found = false;
00177 for (int i = 0; i < nrPortalShapes; ++i) {
00178 ps = PortalShapes[shapeidx]->Sample(p, u1, u2, &ns);
00179 *wi = Normalize(ps - p);
00180
00181 if (Dot(*wi, ns) < 0.f) {
00182 found = true;
00183 break;
00184 }
00185
00186 if (++shapeidx >= nrPortalShapes)
00187 shapeidx = 0;
00188 }
00189 if (found)
00190 *pdf = PortalShapes[shapeidx]->Pdf(p, *wi);
00191 else {
00192 *pdf = 0.f;
00193 return 0.f;
00194 }
00195 }
00196 visibility->SetRay(p, *wi);
00197 return Le(RayDifferential(p, *wi));
00198 }
00199 float SkyLight::Pdf(const Point &, const Normal &n,
00200 const Vector &wi) const {
00201 return AbsDot(n, wi) * INV_TWOPI;
00202 }
00203 SWCSpectrum SkyLight::Sample_L(const Point &p,
00204 float u1, float u2, float u3, Vector *wi, float *pdf,
00205 VisibilityTester *visibility) const {
00206 if(!havePortalShape) {
00207 *wi = UniformSampleSphere(u1, u2);
00208 *pdf = UniformSpherePdf();
00209 } else {
00210
00211 int shapeidx = 0;
00212 if(nrPortalShapes > 1)
00213 shapeidx = min<float>(nrPortalShapes - 1,
00214 Floor2Int(u3 * nrPortalShapes));
00215 Normal ns;
00216 Point ps;
00217 bool found = false;
00218 for (int i = 0; i < nrPortalShapes; ++i) {
00219 ps = PortalShapes[shapeidx]->Sample(p, u1, u2, &ns);
00220 *wi = Normalize(ps - p);
00221 if (Dot(*wi, ns) < 0.f) {
00222 found = true;
00223 break;
00224 }
00225
00226 if (++shapeidx >= nrPortalShapes)
00227 shapeidx = 0;
00228 }
00229 if (found)
00230 *pdf = PortalShapes[shapeidx]->Pdf(p, *wi);
00231 else {
00232 *pdf = 0.f;
00233 return 0.f;
00234 }
00235 }
00236 visibility->SetRay(p, *wi);
00237 return Le(RayDifferential(p, *wi));
00238 }
00239 float SkyLight::Pdf(const Point &, const Vector &) const {
00240 return 1.f / (4.f * M_PI);
00241 }
00242 SWCSpectrum SkyLight::Sample_L(const Scene *scene,
00243 float u1, float u2, float u3, float u4,
00244 Ray *ray, float *pdf) const {
00245 if(!havePortalShape) {
00246
00247 Point worldCenter;
00248 float worldRadius;
00249 scene->WorldBound().BoundingSphere(&worldCenter,
00250 &worldRadius);
00251 worldRadius *= 1.01f;
00252 Point p1 = worldCenter + worldRadius *
00253 UniformSampleSphere(u1, u2);
00254 Point p2 = worldCenter + worldRadius *
00255 UniformSampleSphere(u3, u4);
00256
00257 ray->o = p1;
00258 ray->d = Normalize(p2-p1);
00259
00260 Vector to_center = Normalize(worldCenter - p1);
00261 float costheta = AbsDot(to_center,ray->d);
00262 *pdf =
00263 costheta / ((4.f * M_PI * worldRadius * worldRadius));
00264 } else {
00265
00266
00267 int shapeidx = 0;
00268 if(nrPortalShapes > 1)
00269 shapeidx = min<float>(nrPortalShapes - 1,
00270 Floor2Int(lux::random::floatValue() * nrPortalShapes));
00271
00272 Normal ns;
00273 ray->o = PortalShapes[shapeidx]->Sample(u1, u2, &ns);
00274 ray->d = UniformSampleSphere(u3, u4);
00275 if (Dot(ray->d, ns) < 0.) ray->d *= -1;
00276
00277 *pdf = PortalShapes[shapeidx]->Pdf(ray->o) * INV_TWOPI / nrPortalShapes;
00278 }
00279
00280 return Le(RayDifferential(ray->o, -ray->d));
00281 }
00282 SWCSpectrum SkyLight::Sample_L(const Point &p,
00283 Vector *wi, VisibilityTester *visibility) const {
00284 float pdf;
00285 SWCSpectrum L = Sample_L(p, lux::random::floatValue(), lux::random::floatValue(),
00286 lux::random::floatValue(), wi, &pdf, visibility);
00287 if (pdf == 0.f) return Spectrum(0.f);
00288 return L / pdf;
00289 }
00290 SWCSpectrum SkyLight::Sample_L(const Scene *scene, float u1, float u2, BSDF **bsdf, float *pdf) const
00291 {
00292 Point worldCenter;
00293 float worldRadius;
00294 scene->WorldBound().BoundingSphere(&worldCenter, &worldRadius);
00295 worldRadius *= 1.01f;
00296 Point ps = worldCenter + worldRadius * UniformSampleSphere(u1, u2);
00297 Normal ns = Normal(Normalize(worldCenter - ps));
00298 Vector dpdu, dpdv;
00299 CoordinateSystem(Vector(ns), &dpdu, &dpdv);
00300 DifferentialGeometry dg(ps, ns, dpdu, dpdv, Vector(0, 0, 0), Vector (0, 0, 0), 0, 0, NULL);
00301 *bsdf = BSDF_ALLOC(BSDF)(dg, ns);
00302 (*bsdf)->Add(BSDF_ALLOC(SkyBxDF)(*this, WorldToLight, dpdu, dpdv, Vector(ns)));
00303 *pdf = 1.f / (4.f * M_PI * worldRadius * worldRadius);
00304 return SWCSpectrum(skyScale);
00305 }
00306 SWCSpectrum SkyLight::Sample_L(const Scene *scene, const Point &p, const Normal &n,
00307 float u1, float u2, float u3, BSDF **bsdf, float *pdf, float *pdfDirect,
00308 VisibilityTester *visibility) const
00309 {
00310 Vector wi;
00311 if(!havePortalShape) {
00312
00313 float x, y, z;
00314 ConcentricSampleDisk(u1, u2, &x, &y);
00315 z = sqrtf(max(0.f, 1.f - x*x - y*y));
00316 if (u3 < .5)
00317 z *= -1;
00318 wi = Vector(x, y, z);
00319
00320 *pdfDirect = fabsf(wi.z) * INV_TWOPI;
00321
00322 Vector v1, v2;
00323 CoordinateSystem(Normalize(Vector(n)), &v1, &v2);
00324 wi = Vector(v1.x * wi.x + v2.x * wi.y + n.x * wi.z,
00325 v1.y * wi.x + v2.y * wi.y + n.y * wi.z,
00326 v1.z * wi.x + v2.z * wi.y + n.z * wi.z);
00327 } else {
00328
00329 int shapeIndex = 0;
00330 if(nrPortalShapes > 1)
00331 shapeIndex = Floor2Int(u3 * nrPortalShapes);
00332 Normal ns;
00333 Point ps;
00334 bool found = false;
00335 for (int i = 0; i < nrPortalShapes; ++i) {
00336 ps = PortalShapes[shapeIndex]->Sample(p, u1, u2, &ns);
00337 wi = Normalize(ps - p);
00338 if (Dot(wi, ns) < 0.f) {
00339 found = true;
00340 break;
00341 }
00342 if (++shapeIndex >= nrPortalShapes)
00343 shapeIndex = 0;
00344 }
00345 if (found)
00346 *pdfDirect = PortalShapes[shapeIndex]->Pdf(p, wi);
00347 else {
00348 *pdf = 0.f;
00349 *pdfDirect = 0.f;
00350 return 0.f;
00351 }
00352 }
00353 Point worldCenter;
00354 float worldRadius;
00355 scene->WorldBound().BoundingSphere(&worldCenter, &worldRadius);
00356 worldRadius *= 1.01f;
00357 Vector toCenter(worldCenter - p);
00358 float centerDistance = Dot(toCenter, toCenter);
00359 float approach = Dot(toCenter, wi);
00360 float distance = approach + sqrtf(worldRadius * worldRadius - centerDistance + approach * approach);
00361 Point ps(p + distance * wi);
00362 Normal ns(Normalize(worldCenter - ps));
00363 Vector dpdu, dpdv;
00364 CoordinateSystem(Vector(ns), &dpdu, &dpdv);
00365 DifferentialGeometry dg(ps, ns, dpdu, dpdv, Vector(0, 0, 0), Vector (0, 0, 0), 0, 0, NULL);
00366 *bsdf = BSDF_ALLOC(BSDF)(dg, ns);
00367 (*bsdf)->Add(BSDF_ALLOC(SkyBxDF)(*this, WorldToLight, dpdu, dpdv, Vector(ns)));
00368 *pdf = 1.f / (4.f * M_PI * worldRadius * worldRadius);
00369 *pdfDirect *= AbsDot(wi, ns) / DistanceSquared(p, ps);
00370 visibility->SetSegment(p, ps);
00371 return SWCSpectrum(skyScale);
00372 }
00373
00374 Light* SkyLight::CreateLight(const Transform &light2world,
00375 const ParamSet ¶mSet) {
00376 float scale = paramSet.FindOneFloat("gain", 1.f);
00377 int nSamples = paramSet.FindOneInt("nsamples", 1);
00378 Vector sundir = paramSet.FindOneVector("sundir", Vector(0,0,1));
00379 Normalize(sundir);
00380 float turb = paramSet.FindOneFloat("turbidity", 2.0f);
00381 float aconst = paramSet.FindOneFloat("aconst", 1.0f);
00382 float bconst = paramSet.FindOneFloat("bconst", 1.0f);
00383 float cconst = paramSet.FindOneFloat("cconst", 1.0f);
00384 float dconst = paramSet.FindOneFloat("dconst", 1.0f);
00385 float econst = paramSet.FindOneFloat("econst", 1.0f);
00386
00387 return new SkyLight(light2world, scale, nSamples, sundir, turb, aconst, bconst, cconst, dconst, econst);
00388 }
00389
00390
00391 inline float RiAngleBetween(float thetav, float phiv, float theta, float phi)
00392 {
00393 float cospsi = sin(thetav) * sin(theta) * cos(phi-phiv) + cos(thetav) * cos(theta);
00394 if (cospsi > 1) return 0;
00395 if (cospsi < -1) return M_PI;
00396 return acos(cospsi);
00397 }
00398
00399
00400
00401
00402
00403
00404
00405
00406 void SkyLight::InitSunThetaPhi()
00407 {
00408 Vector wh = Normalize(sundir);
00409 phiS = SphericalPhi(wh);
00410 thetaS = SphericalTheta(wh);
00411 }
00412
00413
00414
00415
00416
00417
00418
00419 inline float SkyLight::PerezFunction(const float *lam, float theta, float gamma, float lvz) const
00420 {
00421 float den = ((1 + lam[1]*exp(lam[2])) *
00422 (1 + lam[3]*exp(lam[4]*thetaS) + lam[5]*cos(thetaS)*cos(thetaS)));
00423 float num = ((1 + lam[1]*exp(lam[2]/cos(theta))) *
00424 (1 + lam[3]*exp(lam[4]*gamma) + lam[5]*cos(gamma)*cos(gamma)));
00425 return lvz* num/den;
00426 }
00427
00428
00429 void SkyLight::GetSkySpectralRadiance(const float theta, const float phi, SWCSpectrum * const dst_spect) const
00430 {
00431
00432 const float theta_fin = min(theta,(M_PI * 0.5f) - 0.001f);
00433 const float gamma = RiAngleBetween(theta,phi,thetaS,phiS);
00434
00435
00436 const float x = PerezFunction(perez_x, theta_fin, gamma, zenith_x);
00437 const float y = PerezFunction(perez_y, theta_fin, gamma, zenith_y);
00438 const float Y = PerezFunction(perez_Y, theta_fin, gamma, zenith_Y);
00439
00440 ChromaticityToSpectrum(x,y,dst_spect);
00441 *dst_spect *= (Y / dst_spect->y() * 0.00000260f);
00442 }
00443
00444
00445 extern boost::thread_specific_ptr<SpectrumWavelengths> thread_wavelengths;
00446
00447
00448 void SkyLight::ChromaticityToSpectrum(const float x, const float y, SWCSpectrum * const dst_spect) const
00449 {
00450 const float den = 1.0f / (0.0241f + 0.2562f * x - 0.7341f * y);
00451 const float M1 = (-1.3515f - 1.7703f * x + 5.9114f * y) * den;
00452 const float M2 = ( 0.03f - 31.4424f * x + 30.0717f * y) * den;
00453
00454 for (unsigned int j = 0; j < WAVELENGTH_SAMPLES; ++j)
00455 {
00456 const float w = (thread_wavelengths->w[j] - 300.0f) * 0.1018867f;
00457 const int i = Floor2Int(w);
00458 const int i1 = i + 1;
00459
00460 const float b = w - float(i);
00461 const float a = 1.0f - b;
00462
00463 const float t0 = S0Amplitudes[i] * a + S0Amplitudes[i1] * b;
00464 const float t1 = S1Amplitudes[i] * a + S1Amplitudes[i1] * b;
00465 const float t2 = S2Amplitudes[i] * a + S2Amplitudes[i1] * b;
00466
00467 dst_spect->c[j] = t0 + M1 * t1 + M2 * t2;
00468 }
00469 }