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
00032 #include "data/sun_spect.h"
00033
00034 using namespace lux;
00035
00036 class SunBxDF : public BxDF
00037 {
00038 public:
00039 SunBxDF(float cosMax, float radius) : BxDF(BxDFType(BSDF_REFLECTION | BSDF_DIFFUSE)), cosThetaMax(cosMax), worldRadius(radius) {}
00040 SWCSpectrum f(const Vector &wo, const Vector &wi) const {return min(wo.z, wi.z) < cosThetaMax ? 0.f : 1.f;}
00041 SWCSpectrum Sample_f(const Vector &wo, Vector *wi, float u1, float u2, float *pdf, float *pdfBack = NULL) const
00042 {
00043 *wi = UniformSampleCone(u1, u2, cosThetaMax);
00044 *pdf = UniformConePdf(cosThetaMax);
00045 if (pdfBack)
00046 *pdfBack = Pdf(*wi, wo);
00047 return 1.f;
00048 }
00049 float Pdf(const Vector &wi, const Vector &wo) const
00050 {
00051 if (min(wi.z, wo.z) < cosThetaMax)
00052 return 0.;
00053 else
00054 return UniformConePdf(cosThetaMax);
00055 }
00056 private:
00057 float cosThetaMax, worldRadius;
00058 };
00059
00060
00061 SunLight::SunLight(const Transform &light2world,
00062 const float sunscale, const Vector &dir, float turb , float relSize, int ns)
00063 : Light(light2world, ns) {
00064 sundir = Normalize(LightToWorld(dir));
00065 turbidity = turb;
00066
00067 CoordinateSystem(sundir, &x, &y);
00068
00069
00070
00071 const float sunRadius = 695500;
00072 const float sunMeanDistance = 149600000;
00073 if(relSize*sunRadius <= sunMeanDistance) {
00074 cosThetaMax = sqrt(1.0f - pow(relSize*sunRadius/sunMeanDistance, 2));
00075 } else {
00076 std::stringstream ss;
00077 ss <<"Reducing relative sun size to "<< sunMeanDistance/sunRadius;
00078 luxError(LUX_LIMIT, LUX_WARNING, ss.str().c_str());
00079 cosThetaMax = 0.0f;
00080 }
00081
00082 float solidAngle = 2*M_PI*(1-cosThetaMax);
00083
00084 Vector wh = Normalize(sundir);
00085 phiS = SphericalPhi(wh);
00086 thetaS = SphericalTheta(wh);
00087
00088
00089 SPD *k_oCurve = new IrregularSPD(sun_k_oWavelengths,sun_k_oAmplitudes, 64);
00090 SPD *k_gCurve = new IrregularSPD(sun_k_gWavelengths, sun_k_gAmplitudes, 2);
00091 SPD *k_waCurve = new IrregularSPD(sun_k_waWavelengths,sun_k_waAmplitudes, 13);
00092
00093 SPD *solCurve = new RegularSPD(sun_sun_irradiance, 380, 770, 79);
00094
00095 float beta = 0.04608365822050 * turbidity - 0.04586025928522;
00096 float tauR, tauA, tauO, tauG, tauWA;
00097
00098 float m = 1.0/(cos(thetaS) + 0.000940 * pow(1.6386 - thetaS,-1.253));
00099
00100 int i;
00101 float lambda;
00102
00103 float Ldata[91];
00104 for(i = 0, lambda = 350; i < 91; i++, lambda+=5) {
00105
00106 tauR = exp( -m * 0.008735 * pow(lambda/1000, float(-4.08)));
00107
00108
00109
00110 const float alpha = 1.3;
00111 tauA = exp(-m * beta * pow(lambda/1000, -alpha));
00112
00113
00114 const float lOzone = .35;
00115 tauO = exp(-m * k_oCurve->sample(lambda) * lOzone);
00116
00117 tauG = exp(-1.41 * k_gCurve->sample(lambda) * m / pow(1 + 118.93 * k_gCurve->sample(lambda) * m, 0.45));
00118
00119
00120 const float w = 2.0;
00121 tauWA = exp(-0.2385 * k_waCurve->sample(lambda) * w * m /
00122 pow(1 + 20.07 * k_waCurve->sample(lambda) * w * m, 0.45));
00123
00124
00125 const float unitConv = 1./(solidAngle*1000000000.);
00126 Ldata[i] = (solCurve->sample(lambda) * tauR * tauA * tauO * tauG * tauWA * unitConv);
00127 }
00128 LSPD = new RegularSPD(Ldata, 350,800,91);
00129 LSPD->Scale(sunscale);
00130
00131 delete k_oCurve;
00132 delete k_gCurve;
00133 delete k_waCurve;
00134 delete solCurve;
00135 }
00136 SWCSpectrum SunLight::Le(const RayDifferential &r) const {
00137 Vector w = r.d;
00138 if(cosThetaMax < 1.0f && Dot(w,sundir) > cosThetaMax)
00139 return LSPD;
00140 else
00141 return 0.;
00142 }
00143 SWCSpectrum SunLight::Le(const Scene *scene, const Ray &r,
00144 const Normal &n, BSDF **bsdf, float *pdf, float *pdfDirect) const
00145 {
00146 if (cosThetaMax == 1.f || Dot(r.d, sundir) < cosThetaMax) {
00147 *bsdf = NULL;
00148 *pdf = 0.f;
00149 *pdfDirect = 0.f;
00150 return 0.f;
00151 }
00152 Point worldCenter;
00153 float worldRadius;
00154 scene->WorldBound().BoundingSphere(&worldCenter, &worldRadius);
00155 Vector toCenter(worldCenter - r.o);
00156 float approach = Dot(toCenter, sundir);
00157 float distance = approach + worldRadius;
00158 Point ps(r.o + distance * r.d);
00159 Normal ns(-sundir);
00160 DifferentialGeometry dg(r.o, ns, -x, y, Vector(0, 0, 0), Vector (0, 0, 0), 0, 0, NULL);
00161 *bsdf = BSDF_ALLOC(BSDF)(dg, ns);
00162 (*bsdf)->Add(BSDF_ALLOC(SunBxDF)(cosThetaMax, worldRadius));
00163 *pdf = 1.f / (M_PI * worldRadius * worldRadius);
00164 *pdfDirect = UniformConePdf(cosThetaMax) * AbsDot(r.d, ns) / DistanceSquared(r.o, ps);
00165 return LSPD;
00166 }
00167
00168 SWCSpectrum SunLight::Sample_L(const Point &p, float u1, float u2, float u3,
00169 Vector *wi, float *pdf, VisibilityTester *visibility) const {
00170 if(cosThetaMax == 1) {
00171 *pdf = 1.f;
00172 return Sample_L(p, wi, visibility);
00173 } else {
00174 *wi = UniformSampleCone(u1, u2, cosThetaMax, x, y, sundir);
00175 *pdf = UniformConePdf(cosThetaMax);
00176 visibility->SetRay(p, *wi);
00177 return LSPD;
00178 }
00179 }
00180 SWCSpectrum SunLight::Sample_L(const Point &p,
00181 Vector *wi, VisibilityTester *visibility) const {
00182 if(cosThetaMax == 1) {
00183 *wi = sundir;
00184 visibility->SetRay(p, *wi);
00185 return LSPD;
00186 } else {
00187 float pdf;
00188 SWCSpectrum Le = Sample_L(p, lux::random::floatValue(), lux::random::floatValue(),
00189 lux::random::floatValue(), wi, &pdf, visibility);
00190 if (pdf == 0.f) return Spectrum(0.f);
00191 return Le / pdf;
00192 }
00193 }
00194 float SunLight::Pdf(const Point &, const Vector &) const {
00195 if(cosThetaMax == 1) {
00196 return 0.;
00197 } else {
00198 return UniformConePdf(cosThetaMax);
00199 }
00200 }
00201 SWCSpectrum SunLight::Sample_L(const Scene *scene,
00202 float u1, float u2, float u3, float u4,
00203 Ray *ray, float *pdf) const {
00204
00205 Point worldCenter;
00206 float worldRadius;
00207 scene->WorldBound().BoundingSphere(&worldCenter, &worldRadius);
00208 float d1, d2;
00209 ConcentricSampleDisk(u1, u2, &d1, &d2);
00210 Point Pdisk =
00211 worldCenter + worldRadius * (d1 * x + d2 * y);
00212
00213 ray->o = Pdisk + worldRadius * sundir;
00214 ray->d = -UniformSampleCone(u3, u4, cosThetaMax, x, y, sundir);
00215 *pdf = UniformConePdf(cosThetaMax) / (M_PI * worldRadius * worldRadius);
00216 return LSPD;
00217 }
00218 SWCSpectrum SunLight::Sample_L(const Scene *scene, float u1, float u2, BSDF **bsdf, float *pdf) const
00219 {
00220 Point worldCenter;
00221 float worldRadius;
00222 scene->WorldBound().BoundingSphere(&worldCenter, &worldRadius);
00223 float d1, d2;
00224 ConcentricSampleDisk(u1, u2, &d1, &d2);
00225 Point ps = worldCenter + worldRadius * (sundir + d1 * x + d2 * y);
00226 Normal ns(-sundir);
00227 DifferentialGeometry dg(ps, ns, -x, y, Vector(0, 0, 0), Vector(0, 0, 0), 0, 0, NULL);
00228 *bsdf = BSDF_ALLOC(BSDF)(dg, ns);
00229 (*bsdf)->Add(BSDF_ALLOC(SunBxDF)(cosThetaMax, worldRadius));
00230 *pdf = 1.f / (M_PI * worldRadius * worldRadius);
00231 return LSPD;
00232 }
00233 SWCSpectrum SunLight::Sample_L(const Scene *scene, const Point &p, const Normal &n,
00234 float u1, float u2, float u3, BSDF **bsdf, float *pdf, float *pdfDirect,
00235 VisibilityTester *visibility) const
00236 {
00237 Vector wi;
00238 if(cosThetaMax == 1) {
00239 wi = sundir;
00240 *pdfDirect = 1.f;
00241 } else {
00242 wi = UniformSampleCone(u1, u2, cosThetaMax, x, y, sundir);
00243 *pdfDirect = UniformConePdf(cosThetaMax);
00244 }
00245 Point worldCenter;
00246 float worldRadius;
00247 scene->WorldBound().BoundingSphere(&worldCenter, &worldRadius);
00248 Vector toCenter(worldCenter - p);
00249 float approach = Dot(toCenter, sundir);
00250 float distance = approach + worldRadius;
00251 Point ps(p + distance * wi);
00252 Normal ns(-sundir);
00253 DifferentialGeometry dg(ps, ns, -x, y, Vector(0, 0, 0), Vector (0, 0, 0), 0, 0, NULL);
00254 *bsdf = BSDF_ALLOC(BSDF)(dg, ns);
00255 (*bsdf)->Add(BSDF_ALLOC(SunBxDF)(cosThetaMax, worldRadius));
00256 *pdf = 1.f / (M_PI * worldRadius * worldRadius);
00257 *pdfDirect *= AbsDot(wi, ns) / DistanceSquared(p, ps);
00258 visibility->SetSegment(p, ps);
00259 return LSPD;
00260 }
00261
00262 Light* SunLight::CreateLight(const Transform &light2world,
00263 const ParamSet ¶mSet) {
00264
00265 float scale = paramSet.FindOneFloat("gain", 1.f);
00266 int nSamples = paramSet.FindOneInt("nsamples", 1);
00267 Vector sundir = paramSet.FindOneVector("sundir", Vector(0,0,-1));
00268 float turb = paramSet.FindOneFloat("turbidity", 2.0f);
00269 float relSize = paramSet.FindOneFloat("relsize", 1.0f);
00270 return new SunLight(light2world, scale, sundir, turb, relSize, nSamples);
00271 }