00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025 #include "pbrt.h"
00026 #include "transport.h"
00027 #include "scene.h"
00028 #include "mc.h"
00029 #include "octree.h"
00030
00031 struct IrradianceSample;
00032 struct IrradProcess;
00033
00034 class IrradianceCache : public SurfaceIntegrator {
00035 public:
00036
00037 IrradianceCache(int maxspec, int maxind, float maxerr, int nsamples);
00038 ~IrradianceCache();
00039 Spectrum Li(const Scene *scene, const RayDifferential &ray, const Sample *sample, float *alpha) const;
00040 void RequestSamples(Sample *sample, const Scene *scene);
00041 void Preprocess(const Scene *);
00042 private:
00043
00044 float maxError;
00045 int nSamples;
00046 int maxSpecularDepth, maxIndirectDepth;
00047 mutable int specularDepth;
00048
00049 int *lightSampleOffset, lightNumOffset;
00050 int *bsdfSampleOffset, *bsdfComponentOffset;
00051 mutable Octree<IrradianceSample, IrradProcess> *octree;
00052
00053 Spectrum IndirectLo(const Point &p, const Normal &n,
00054 const Vector &wo, BSDF *bsdf, BxDFType flags,
00055 const Sample *sample, const Scene *scene) const;
00056 bool InterpolateIrradiance(const Scene *scene,
00057 const Point &p, const Normal &n, Spectrum *E) const;
00058 };
00059 struct IrradianceSample {
00060
00061 IrradianceSample() { }
00062 IrradianceSample(const Spectrum &e, const Point &P, const Normal &N,
00063 float md) : E(e), n(N), p(P) {
00064 maxDist = md;
00065 }
00066 Spectrum E;
00067 Normal n;
00068 Point p;
00069 float maxDist;
00070 };
00071 struct IrradProcess {
00072
00073 IrradProcess(const Normal &N, float me) {
00074 n = N;
00075 maxError = me;
00076 nFound = samplesChecked = 0;
00077 sumWt = 0.;
00078 E = 0.;
00079 }
00080 void operator()(const Point &P, const IrradianceSample &sample) const;
00081 bool Successful() {
00082 return (sumWt > 0. && nFound > 0);
00083 }
00084 Spectrum GetIrradiance() const { return E / sumWt; }
00085 Normal n;
00086 float maxError;
00087 mutable int nFound, samplesChecked;
00088 mutable float sumWt;
00089 mutable Spectrum E;
00090 };
00091
00092 IrradianceCache::IrradianceCache(int maxspec, int maxind,
00093 float maxerr, int ns) {
00094 maxError = maxerr;
00095 nSamples = ns;
00096 maxSpecularDepth = maxspec;
00097 maxIndirectDepth = maxind;
00098 specularDepth = 0;
00099 }
00100 void IrradianceCache::RequestSamples(Sample *sample,
00101 const Scene *scene) {
00102
00103 u_int nLights = scene->lights.size();
00104 lightSampleOffset = new int[nLights];
00105 bsdfSampleOffset = new int[nLights];
00106 bsdfComponentOffset = new int[nLights];
00107 for (u_int i = 0; i < nLights; ++i) {
00108 const Light *light = scene->lights[i];
00109 int lightSamples =
00110 scene->sampler->RoundSize(light->nSamples);
00111 lightSampleOffset[i] = sample->Add2D(lightSamples);
00112 bsdfSampleOffset[i] = sample->Add2D(lightSamples);
00113 bsdfComponentOffset[i] = sample->Add1D(lightSamples);
00114 }
00115 lightNumOffset = -1;
00116 }
00117 Spectrum IrradianceCache::Li(const Scene *scene, const RayDifferential &ray,
00118 const Sample *sample, float *alpha) const {
00119 Intersection isect;
00120 Spectrum L(0.);
00121 if (scene->Intersect(ray, &isect)) {
00122 if (alpha) *alpha = 1.;
00123
00124 BSDF *bsdf = isect.GetBSDF(ray);
00125 Vector wo = -ray.d;
00126 const Point &p = bsdf->dgShading.p;
00127 const Normal &n = bsdf->dgShading.nn;
00128
00129 L += isect.Le(wo);
00130 L += UniformSampleAllLights(scene, p, n, wo, bsdf, sample,
00131 lightSampleOffset, bsdfSampleOffset,
00132 bsdfComponentOffset);
00133
00134 if (specularDepth++ < maxSpecularDepth) {
00135 Vector wi;
00136
00137 Spectrum f = bsdf->Sample_f(wo, &wi,
00138 BxDFType(BSDF_REFLECTION | BSDF_SPECULAR));
00139 if (!f.Black()) {
00140
00141 RayDifferential rd(p, wi);
00142 rd.hasDifferentials = true;
00143 rd.rx.o = p + isect.dg.dpdx;
00144 rd.ry.o = p + isect.dg.dpdy;
00145
00146 Normal dndx = bsdf->dgShading.dndu * bsdf->dgShading.dudx +
00147 bsdf->dgShading.dndv * bsdf->dgShading.dvdx;
00148 Normal dndy = bsdf->dgShading.dndu * bsdf->dgShading.dudy +
00149 bsdf->dgShading.dndv * bsdf->dgShading.dvdy;
00150 Vector dwodx = -ray.rx.d - wo, dwody = -ray.ry.d - wo;
00151 float dDNdx = Dot(dwodx, n) + Dot(wo, dndx);
00152 float dDNdy = Dot(dwody, n) + Dot(wo, dndy);
00153 rd.rx.d = wi -
00154 dwodx + 2 * Vector(Dot(wo, n) * dndx +
00155 dDNdx * n);
00156 rd.ry.d = wi -
00157 dwody + 2 * Vector(Dot(wo, n) * dndy +
00158 dDNdy * n);
00159 L += scene->Li(rd, sample) * f * AbsDot(wi, n);
00160 }
00161 f = bsdf->Sample_f(wo, &wi,
00162 BxDFType(BSDF_TRANSMISSION | BSDF_SPECULAR));
00163 if (!f.Black()) {
00164
00165 RayDifferential rd(p, wi);
00166 rd.hasDifferentials = true;
00167 rd.rx.o = p + isect.dg.dpdx;
00168 rd.ry.o = p + isect.dg.dpdy;
00169
00170 float eta = bsdf->eta;
00171 Vector w = -wo;
00172 if (Dot(wo, n) < 0) eta = 1.f / eta;
00173
00174 Normal dndx = bsdf->dgShading.dndu * bsdf->dgShading.dudx + bsdf->dgShading.dndv * bsdf->dgShading.dvdx;
00175 Normal dndy = bsdf->dgShading.dndu * bsdf->dgShading.dudy + bsdf->dgShading.dndv * bsdf->dgShading.dvdy;
00176
00177 Vector dwodx = -ray.rx.d - wo, dwody = -ray.ry.d - wo;
00178 float dDNdx = Dot(dwodx, n) + Dot(wo, dndx);
00179 float dDNdy = Dot(dwody, n) + Dot(wo, dndy);
00180
00181 float mu = eta * Dot(w, n) - Dot(wi, n);
00182 float dmudx = (eta - (eta*eta*Dot(w,n))/Dot(wi, n)) * dDNdx;
00183 float dmudy = (eta - (eta*eta*Dot(w,n))/Dot(wi, n)) * dDNdy;
00184
00185 rd.rx.d = wi + eta * dwodx - Vector(mu * dndx + dmudx * n);
00186 rd.ry.d = wi + eta * dwody - Vector(mu * dndy + dmudy * n);
00187 L += scene->Li(rd, sample) * f * AbsDot(wi, n);
00188 }
00189 }
00190 --specularDepth;
00191
00192 Normal ng = isect.dg.nn;
00193 if (Dot(wo, ng) < 0.f) ng = -ng;
00194 BxDFType flags = BxDFType(BSDF_REFLECTION |
00195 BSDF_DIFFUSE |
00196 BSDF_GLOSSY);
00197 L += IndirectLo(p, ng, wo, bsdf, flags, sample, scene);
00198 flags = BxDFType(BSDF_TRANSMISSION |
00199 BSDF_DIFFUSE |
00200 BSDF_GLOSSY);
00201 L += IndirectLo(p, -ng, wo, bsdf, flags, sample, scene);
00202 }
00203 else {
00204
00205 if (alpha) *alpha = 0.;
00206 for (u_int i = 0; i < scene->lights.size(); ++i)
00207 L += scene->lights[i]->Le(ray);
00208 if (alpha && !L.Black()) *alpha = 1.;
00209 return L;
00210 }
00211 return L;
00212 }
00213 Spectrum IrradianceCache::IndirectLo(const Point &p,
00214 const Normal &n, const Vector &wo, BSDF *bsdf,
00215 BxDFType flags, const Sample *sample,
00216 const Scene *scene) const {
00217 if (bsdf->NumComponents(flags) == 0)
00218 return Spectrum(0.);
00219 Spectrum E;
00220 if (!InterpolateIrradiance(scene, p, n, &E)) {
00221
00222 u_int scramble[2] = { RandomUInt(), RandomUInt() };
00223 float sumInvDists = 0.;
00224 for (int i = 0; i < nSamples; ++i) {
00225
00226
00227 static StatsCounter nIrradiancePaths("Irradiance Cache",
00228 "Paths followed for irradiance estimates");
00229 ++nIrradiancePaths;
00230 float u[2];
00231 Sample02(i, scramble, u);
00232 Vector w = CosineSampleHemisphere(u[0], u[1]);
00233 RayDifferential r(p, bsdf->LocalToWorld(w));
00234 if (Dot(r.d, n) < 0) r.d = -r.d;
00235 Spectrum L(0.);
00236
00237 {
00238
00239 Spectrum pathThroughput = 1.;
00240 RayDifferential ray(r);
00241 bool specularBounce = false;
00242 for (int pathLength = 0; ; ++pathLength) {
00243
00244 Intersection isect;
00245 if (!scene->Intersect(ray, &isect))
00246 break;
00247 if (pathLength == 0)
00248 r.maxt = ray.maxt;
00249 pathThroughput *= scene->Transmittance(ray);
00250
00251 if (specularBounce)
00252 L += pathThroughput * isect.Le(-ray.d);
00253
00254 BSDF *bsdf = isect.GetBSDF(ray);
00255
00256 const Point &p = bsdf->dgShading.p;
00257 const Normal &n = bsdf->dgShading.nn;
00258 Vector wo = -ray.d;
00259 L += pathThroughput *
00260 UniformSampleOneLight(scene, p, n, wo, bsdf, sample);
00261 if (pathLength+1 == maxIndirectDepth) break;
00262
00263
00264 float bs1 = RandomFloat(), bs2 = RandomFloat(), bcs = RandomFloat();
00265 Vector wi;
00266 float pdf;
00267 BxDFType flags;
00268 Spectrum f = bsdf->Sample_f(wo, &wi, bs1, bs2, bcs,
00269 &pdf, BSDF_ALL, &flags);
00270 if (f.Black() || pdf == 0.)
00271 break;
00272 specularBounce = (flags & BSDF_SPECULAR) != 0;
00273 pathThroughput *= f * AbsDot(wi, n) / pdf;
00274 ray = RayDifferential(p, wi);
00275
00276 if (pathLength > 3) {
00277 float continueProbability = .5f;
00278 if (RandomFloat() > continueProbability)
00279 break;
00280 pathThroughput /= continueProbability;
00281 }
00282 }
00283 }
00284 E += L;
00285 float dist = r.maxt * r.d.Length();
00286 sumInvDists += 1.f / dist;
00287 }
00288 E *= M_PI / float(nSamples);
00289
00290
00291 static StatsCounter nSamplesComputed("Irradiance Cache",
00292 "Irradiance estimates computed");
00293 ++nSamplesComputed;
00294
00295 static float minMaxDist =
00296 .001f * powf(scene->WorldBound().Volume(), 1.f/3.f);
00297 static float maxMaxDist =
00298 .125f * powf(scene->WorldBound().Volume(), 1.f/3.f);
00299 float maxDist = nSamples / sumInvDists;
00300 if (minMaxDist > 0.f)
00301 maxDist = Clamp(maxDist, minMaxDist, maxMaxDist);
00302 maxDist *= maxError;
00303 BBox sampleExtent(p);
00304 sampleExtent.Expand(maxDist);
00305 octree->Add(IrradianceSample(E, p, n, maxDist),
00306 sampleExtent);
00307 }
00308 return .5f * bsdf->rho(wo, flags) * E;
00309 }
00310 void IrradianceCache::Preprocess(const Scene *scene) {
00311 BBox wb = scene->WorldBound();
00312 Vector delta = .01f * (wb.pMax - wb.pMin);
00313 wb.pMin -= delta;
00314 wb.pMax += delta;
00315 octree = new Octree<IrradianceSample, IrradProcess>(wb);
00316 }
00317 IrradianceCache::~IrradianceCache() {
00318 delete octree;
00319 }
00320 bool
00321 IrradianceCache::InterpolateIrradiance(const Scene *scene,
00322 const Point &p, const Normal &n, Spectrum *E) const {
00323 if (!octree) return false;
00324 IrradProcess proc(n, maxError);
00325 octree->Lookup(p, proc);
00326
00327 static StatsPercentage nSuccessfulLookups("Irradiance Cache",
00328 "Successful irradiance cache lookups");
00329 static StatsRatio nSamplesFound("Irradiance Cache",
00330 "Irradiance samples found per successful lookup");
00331 static StatsRatio nSamplesChecked("Irradiance Cache",
00332 "Irradiance samples checked per lookup");
00333 nSuccessfulLookups.Add(proc.Successful() ? 1 : 0, 1);
00334 nSamplesFound.Add(proc.nFound, 1);
00335 nSamplesChecked.Add(proc.samplesChecked, 1);
00336 if (!proc.Successful()) return false;
00337 *E = proc.GetIrradiance();
00338 return true;
00339 }
00340 void IrradProcess::operator()(const Point &p,
00341 const IrradianceSample &sample) const {
00342 ++samplesChecked;
00343
00344 if (Dot(n, sample.n) < 0.01f)
00345 return;
00346
00347 float d2 = DistanceSquared(p, sample.p);
00348 if (d2 > sample.maxDist * sample.maxDist)
00349 return;
00350
00351 Normal navg = sample.n + n;
00352 if (Dot(p - sample.p, navg) < -.01f)
00353 return;
00354
00355 float err = sqrtf(d2) / (sample.maxDist * Dot(n, sample.n));
00356 if (err < 1.) {
00357 ++nFound;
00358 float wt = (1.f - err) * (1.f - err);
00359 E += wt * sample.E;
00360 sumWt += wt;
00361 }
00362 }
00363 extern "C" DLLEXPORT SurfaceIntegrator *CreateSurfaceIntegrator(const ParamSet ¶ms) {
00364 float maxError = params.FindOneFloat("maxerror", .2f);
00365 int maxSpecularDepth = params.FindOneInt("maxspeculardepth", 5);
00366 int maxIndirectDepth = params.FindOneInt("maxindirectdepth", 3);
00367 int nSamples = params.FindOneInt("nsamples", 4096);
00368 return new IrradianceCache(maxSpecularDepth, maxIndirectDepth,
00369 maxError, nSamples);
00370 }