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
00026 #include "stdafx.h"
00027 #include "renderers/metropolis.h"
00028 #include "renderers/samplerrenderer.h"
00029 #include "scene.h"
00030 #include "imageio.h"
00031 #include "spectrum.h"
00032 #include "camera.h"
00033 #include "film.h"
00034 #include "rng.h"
00035 #include "progressreporter.h"
00036 #include "paramset.h"
00037 #include "parallel.h"
00038 #include "probes.h"
00039 #include "intersection.h"
00040 #include "montecarlo.h"
00041 #include "samplers/lowdiscrepancy.h"
00042 #include "integrators/directlighting.h"
00043
00044
00045 struct PathSample {
00046 BSDFSample bsdfSample;
00047 float rrSample;
00048 };
00049
00050
00051 struct LightingSample {
00052 BSDFSample bsdfSample;
00053 float lightNum;
00054 LightSample lightSample;
00055 };
00056
00057
00058 struct MLTSample {
00059 MLTSample(int maxLength) {
00060 cameraPathSamples.resize(maxLength);
00061 lightPathSamples.resize(maxLength);
00062 lightingSamples.resize(maxLength);
00063 }
00064 CameraSample cameraSample;
00065 float lightNumSample, lightRaySamples[5];
00066 vector<PathSample> cameraPathSamples, lightPathSamples;
00067 vector<LightingSample> lightingSamples;
00068 };
00069
00070
00071 static void LargeStep(RNG &rng, MLTSample *sample, int maxDepth,
00072 float x, float y, float t0, float t1, bool bidirectional) {
00073
00074 sample->cameraSample.imageX = x;
00075 sample->cameraSample.imageY = y;
00076 sample->cameraSample.time = Lerp(rng.RandomFloat(), t0, t1);
00077 sample->cameraSample.lensU = rng.RandomFloat();
00078 sample->cameraSample.lensV = rng.RandomFloat();
00079 for (int i = 0; i < maxDepth; ++i) {
00080
00081 PathSample &cps = sample->cameraPathSamples[i];
00082 cps.bsdfSample.uComponent = rng.RandomFloat();
00083 cps.bsdfSample.uDir[0] = rng.RandomFloat();
00084 cps.bsdfSample.uDir[1] = rng.RandomFloat();
00085 cps.rrSample = rng.RandomFloat();
00086
00087
00088 LightingSample &ls = sample->lightingSamples[i];
00089 ls.bsdfSample.uComponent = rng.RandomFloat();
00090 ls.bsdfSample.uDir[0] = rng.RandomFloat();
00091 ls.bsdfSample.uDir[1] = rng.RandomFloat();
00092 ls.lightNum = rng.RandomFloat();
00093 ls.lightSample.uComponent = rng.RandomFloat();
00094 ls.lightSample.uPos[0] = rng.RandomFloat();
00095 ls.lightSample.uPos[1] = rng.RandomFloat();
00096 }
00097 if (bidirectional) {
00098
00099 sample->lightNumSample = rng.RandomFloat();
00100 for (int i = 0; i < 5; ++i)
00101 sample->lightRaySamples[i] = rng.RandomFloat();
00102 for (int i = 0; i < maxDepth; ++i) {
00103
00104 PathSample &lps = sample->lightPathSamples[i];
00105 lps.bsdfSample.uComponent = rng.RandomFloat();
00106 lps.bsdfSample.uDir[0] = rng.RandomFloat();
00107 lps.bsdfSample.uDir[1] = rng.RandomFloat();
00108 lps.rrSample = rng.RandomFloat();
00109 }
00110 }
00111 }
00112
00113
00114 static inline void mutate(RNG &rng, float *v, float min = 0.f,
00115 float max = 1.f) {
00116 if (min == max) { *v = min; return; }
00117 Assert(min < max);
00118 float a = 1.f / 1024.f, b = 1.f / 64.f;
00119 static const float logRatio = -logf(b/a);
00120 float delta = (max - min) * b * expf(logRatio * rng.RandomFloat());
00121 if (rng.RandomFloat() < 0.5f) {
00122 *v += delta;
00123 if (*v >= max) *v = min + (*v - max);
00124 }
00125 else {
00126 *v -= delta;
00127 if (*v < min) *v = max - (min - *v);
00128 }
00129 if (*v < min || *v >= max) *v = min;
00130 }
00131
00132
00133 static void SmallStep(RNG &rng, MLTSample *sample, int maxDepth,
00134 int x0, int x1, int y0, int y1, float t0, float t1,
00135 bool bidirectional) {
00136 mutate(rng, &sample->cameraSample.imageX, x0, x1);
00137 mutate(rng, &sample->cameraSample.imageY, y0, y1);
00138 mutate(rng, &sample->cameraSample.time, t0, t1);
00139 mutate(rng, &sample->cameraSample.lensU);
00140 mutate(rng, &sample->cameraSample.lensV);
00141
00142 for (int i = 0; i < maxDepth; ++i) {
00143
00144 PathSample &eps = sample->cameraPathSamples[i];
00145 mutate(rng, &eps.bsdfSample.uComponent);
00146 mutate(rng, &eps.bsdfSample.uDir[0]);
00147 mutate(rng, &eps.bsdfSample.uDir[1]);
00148 mutate(rng, &eps.rrSample);
00149
00150
00151 LightingSample &ls = sample->lightingSamples[i];
00152 mutate(rng, &ls.bsdfSample.uComponent);
00153 mutate(rng, &ls.bsdfSample.uDir[0]);
00154 mutate(rng, &ls.bsdfSample.uDir[1]);
00155 mutate(rng, &ls.lightNum);
00156 mutate(rng, &ls.lightSample.uComponent);
00157 mutate(rng, &ls.lightSample.uPos[0]);
00158 mutate(rng, &ls.lightSample.uPos[1]);
00159 }
00160
00161 if (bidirectional) {
00162 mutate(rng, &sample->lightNumSample);
00163 for (int i = 0; i < 5; ++i)
00164 mutate(rng, &sample->lightRaySamples[i]);
00165 for (int i = 0; i < maxDepth; ++i) {
00166
00167 PathSample &lps = sample->lightPathSamples[i];
00168 mutate(rng, &lps.bsdfSample.uComponent);
00169 mutate(rng, &lps.bsdfSample.uDir[0]);
00170 mutate(rng, &lps.bsdfSample.uDir[1]);
00171 mutate(rng, &lps.rrSample);
00172 }
00173 }
00174 }
00175
00176
00177 struct PathVertex {
00178 Intersection isect;
00179 Vector wPrev, wNext;
00180 BSDF *bsdf;
00181 bool specularBounce;
00182 int nSpecularComponents;
00183 Spectrum alpha;
00184 };
00185
00186
00187 static uint32_t GeneratePath(const RayDifferential &r, const Spectrum &alpha,
00188 const Scene *scene, MemoryArena &arena, const vector<PathSample> &samples,
00189 PathVertex *path, RayDifferential *escapedRay,
00190 Spectrum *escapedAlpha);
00191 inline float I(const Spectrum &L);
00192 class MLTTask : public Task {
00193 public:
00194 MLTTask(ProgressReporter &prog, uint32_t pfreq, uint32_t taskNum,
00195 float dx, float dy, int xx0, int xx1, int yy0, int yy1, float tt0, float tt1,
00196 float bb, const MLTSample &is, const Scene *sc, const Camera *c,
00197 MetropolisRenderer *renderer, Mutex *filmMutex,
00198 Distribution1D *lightDistribution);
00199 void Run();
00200
00201 private:
00202 ProgressReporter &progress;
00203 uint32_t progressUpdateFrequency, taskNum;
00204 float dx, dy;
00205 int currentPixelSample;
00206 int x0, x1, y0, y1;
00207 float t0, t1;
00208 float b;
00209 const MLTSample &initialSample;
00210 const Scene *scene;
00211 const Camera *camera;
00212 MetropolisRenderer *renderer;
00213 Mutex *filmMutex;
00214 Distribution1D *lightDistribution;
00215 };
00216
00217
00218
00219
00220 static uint32_t GeneratePath(const RayDifferential &r,
00221 const Spectrum &a, const Scene *scene, MemoryArena &arena,
00222 const vector<PathSample> &samples, PathVertex *path,
00223 RayDifferential *escapedRay, Spectrum *escapedAlpha) {
00224 PBRT_MLT_STARTED_GENERATE_PATH();
00225 RayDifferential ray = r;
00226 Spectrum alpha = a;
00227 if (escapedAlpha) *escapedAlpha = 0.f;
00228 uint32_t length = 0;
00229 for (; length < samples.size(); ++length) {
00230
00231 PathVertex &v = path[length];
00232 if (!scene->Intersect(ray, &v.isect)) {
00233
00234 if (escapedAlpha) *escapedAlpha = alpha;
00235 if (escapedRay) *escapedRay = ray;
00236 break;
00237 }
00238
00239
00240 v.alpha = alpha;
00241 BSDF *bsdf = v.isect.GetBSDF(ray, arena);
00242 v.bsdf = bsdf;
00243 v.wPrev = -ray.d;
00244
00245
00246 float pdf;
00247 BxDFType flags;
00248 Spectrum f = bsdf->Sample_f(-ray.d, &v.wNext, samples[length].bsdfSample,
00249 &pdf, BSDF_ALL, &flags);
00250 v.specularBounce = (flags & BSDF_SPECULAR) != 0;
00251 v.nSpecularComponents = bsdf->NumComponents(BxDFType(BSDF_SPECULAR |
00252 BSDF_REFLECTION | BSDF_TRANSMISSION));
00253 if (f.IsBlack() || pdf == 0.f)
00254 {
00255 PBRT_MLT_FINISHED_GENERATE_PATH();
00256 return length+1;
00257 }
00258
00259
00260 const Point &p = bsdf->dgShading.p;
00261 const Normal &n = bsdf->dgShading.nn;
00262 Spectrum pathScale = f * AbsDot(v.wNext, n) / pdf;
00263 float rrSurviveProb = min(1.f, pathScale.y());
00264 if (samples[length].rrSample > rrSurviveProb)
00265 {
00266 PBRT_MLT_FINISHED_GENERATE_PATH();
00267 return length+1;
00268 }
00269 alpha *= pathScale / rrSurviveProb;
00270
00271 ray = RayDifferential(p, v.wNext, ray, v.isect.rayEpsilon);
00272 }
00273 PBRT_MLT_FINISHED_GENERATE_PATH();
00274 return length;
00275 }
00276
00277
00278 Spectrum MetropolisRenderer::PathL(const MLTSample &sample,
00279 const Scene *scene, MemoryArena &arena, const Camera *camera,
00280 const Distribution1D *lightDistribution,
00281 PathVertex *cameraPath, PathVertex *lightPath,
00282 RNG &rng) const {
00283
00284 PBRT_STARTED_GENERATING_CAMERA_RAY((CameraSample *)(&sample.cameraSample));
00285 RayDifferential cameraRay;
00286 float cameraWt = camera->GenerateRayDifferential(sample.cameraSample,
00287 &cameraRay);
00288 cameraRay.ScaleDifferentials(1.f / sqrtf(nPixelSamples));
00289 PBRT_FINISHED_GENERATING_CAMERA_RAY((CameraSample *)(&sample.cameraSample), &cameraRay, cameraWt);
00290 RayDifferential escapedRay;
00291 Spectrum escapedAlpha;
00292 uint32_t cameraLength = GeneratePath(cameraRay, cameraWt, scene, arena,
00293 sample.cameraPathSamples, cameraPath, &escapedRay, &escapedAlpha);
00294 if (!bidirectional) {
00295
00296 return Lpath(scene, cameraPath, cameraLength, arena,
00297 sample.lightingSamples, rng, sample.cameraSample.time,
00298 lightDistribution, escapedRay, escapedAlpha);
00299 }
00300 else {
00301
00302
00303
00304 PBRT_MLT_STARTED_SAMPLE_LIGHT_FOR_BIDIR();
00305 float lightPdf, lightRayPdf;
00306 uint32_t lightNum = lightDistribution->SampleDiscrete(sample.lightNumSample,
00307 &lightPdf);
00308 const Light *light = scene->lights[lightNum];
00309 Ray lightRay;
00310 Normal Nl;
00311 LightSample lrs(sample.lightRaySamples[0], sample.lightRaySamples[1],
00312 sample.lightRaySamples[2]);
00313 Spectrum lightWt = light->Sample_L(scene, lrs, sample.lightRaySamples[3],
00314 sample.lightRaySamples[4], sample.cameraSample.time, &lightRay,
00315 &Nl, &lightRayPdf);
00316 PBRT_MLT_FINISHED_SAMPLE_LIGHT_FOR_BIDIR();
00317 if (lightWt.IsBlack() || lightRayPdf == 0.f) {
00318
00319 return Lpath(scene, cameraPath, cameraLength, arena,
00320 sample.lightingSamples, rng, sample.cameraSample.time,
00321 lightDistribution, escapedRay, escapedAlpha);
00322 }
00323 else {
00324
00325 lightWt *= AbsDot(Normalize(Nl), lightRay.d) / (lightPdf * lightRayPdf);
00326 uint32_t lightLength = GeneratePath(RayDifferential(lightRay), lightWt,
00327 scene, arena, sample.lightPathSamples, lightPath, NULL, NULL);
00328
00329 return Lbidir(scene, cameraPath, cameraLength, lightPath, lightLength,
00330 arena, sample.lightingSamples, rng, sample.cameraSample.time,
00331 lightDistribution, escapedRay, escapedAlpha);
00332 }
00333 }
00334 }
00335
00336
00337 Spectrum MetropolisRenderer::Lpath(const Scene *scene,
00338 const PathVertex *cameraPath, int cameraPathLength,
00339 MemoryArena &arena, const vector<LightingSample> &samples,
00340 RNG &rng, float time, const Distribution1D *lightDistribution,
00341 const RayDifferential &eRay, const Spectrum &eAlpha) const {
00342 PBRT_MLT_STARTED_LPATH();
00343 Spectrum L = 0.;
00344 bool previousSpecular = true, allSpecular = true;
00345 for (int i = 0; i < cameraPathLength; ++i) {
00346
00347 const PathVertex &vc = cameraPath[i];
00348 const Point &pc = vc.bsdf->dgShading.p;
00349 const Normal &nc = vc.bsdf->dgShading.nn;
00350
00351
00352 if (previousSpecular && (directLighting == NULL || !allSpecular))
00353 L += vc.alpha * vc.isect.Le(vc.wPrev);
00354
00355
00356 Spectrum Ld(0.f);
00357 if (directLighting == NULL || !allSpecular) {
00358
00359 const LightingSample &ls = samples[i];
00360 float lightPdf;
00361 uint32_t lightNum = lightDistribution->SampleDiscrete(ls.lightNum,
00362 &lightPdf);
00363 const Light *light = scene->lights[lightNum];
00364 PBRT_MLT_STARTED_ESTIMATE_DIRECT();
00365
00366 Ld = vc.alpha *
00367 EstimateDirect(scene, this, arena, light, pc, nc, vc.wPrev,
00368 vc.isect.rayEpsilon, time, vc.bsdf, rng,
00369 ls.lightSample, ls.bsdfSample,
00370 BxDFType(BSDF_ALL & ~BSDF_SPECULAR)) / lightPdf;
00371 PBRT_MLT_FINISHED_ESTIMATE_DIRECT();
00372 }
00373 previousSpecular = vc.specularBounce;
00374 allSpecular &= previousSpecular;
00375 L += Ld;
00376 }
00377
00378 if (!eAlpha.IsBlack() && previousSpecular &&
00379 (directLighting == NULL || !allSpecular))
00380 for (uint32_t i = 0; i < scene->lights.size(); ++i)
00381 L += eAlpha * scene->lights[i]->Le(eRay);
00382 PBRT_MLT_FINISHED_LPATH();
00383 return L;
00384 }
00385
00386
00387 Spectrum MetropolisRenderer::Lbidir(const Scene *scene,
00388 const PathVertex *cameraPath, int cameraPathLength,
00389 const PathVertex *lightPath, int lightPathLength,
00390 MemoryArena &arena, const vector<LightingSample> &samples,
00391 RNG &rng, float time, const Distribution1D *lightDistribution,
00392 const RayDifferential &eRay, const Spectrum &eAlpha) const {
00393 PBRT_MLT_STARTED_LBIDIR();
00394 Spectrum L = 0.;
00395 bool previousSpecular = true, allSpecular = true;
00396
00397 int nVerts = cameraPathLength + lightPathLength + 2;
00398 int *nSpecularVertices = ALLOCA(int, nVerts);
00399 memset(nSpecularVertices, 0, nVerts * sizeof(int));
00400 for (int i = 0; i < cameraPathLength; ++i)
00401 for (int j = 0; j < lightPathLength; ++j)
00402 if (cameraPath[i].specularBounce || lightPath[j].specularBounce)
00403 ++nSpecularVertices[i+j+2];
00404 for (int i = 0; i < cameraPathLength; ++i) {
00405
00406 const PathVertex &vc = cameraPath[i];
00407 const Point &pc = vc.bsdf->dgShading.p;
00408 const Normal &nc = vc.bsdf->dgShading.nn;
00409
00410
00411
00412
00413 if (previousSpecular && (directLighting == NULL || !allSpecular))
00414 L += vc.alpha * vc.isect.Le(vc.wPrev);
00415
00416
00417 Spectrum Ld(0.f);
00418 if (directLighting == NULL || !allSpecular) {
00419
00420 const LightingSample &ls = samples[i];
00421 float lightPdf;
00422 uint32_t lightNum = lightDistribution->SampleDiscrete(ls.lightNum,
00423 &lightPdf);
00424 const Light *light = scene->lights[lightNum];
00425 PBRT_MLT_STARTED_ESTIMATE_DIRECT();
00426
00427 Ld = vc.alpha *
00428 EstimateDirect(scene, this, arena, light, pc, nc, vc.wPrev,
00429 vc.isect.rayEpsilon, time, vc.bsdf, rng,
00430 ls.lightSample, ls.bsdfSample,
00431 BxDFType(BSDF_ALL & ~BSDF_SPECULAR)) / lightPdf;
00432 PBRT_MLT_FINISHED_ESTIMATE_DIRECT();
00433 }
00434 previousSpecular = vc.specularBounce;
00435 allSpecular &= previousSpecular;
00436 L += Ld / (i + 1 - nSpecularVertices[i+1]);
00437 if (!vc.specularBounce) {
00438
00439 for (int j = 0; j < lightPathLength; ++j) {
00440 const PathVertex &vl = lightPath[j];
00441 const Point &pl = vl.bsdf->dgShading.p;
00442 const Normal &nl = vl.bsdf->dgShading.nn;
00443 if (!vl.specularBounce) {
00444
00445 Vector w = Normalize(pl - pc);
00446 Spectrum fc = vc.bsdf->f(vc.wPrev, w) * (1 + vc.nSpecularComponents);
00447 Spectrum fl = vl.bsdf->f(-w, vl.wPrev) * (1 + vl.nSpecularComponents);
00448 if (fc.IsBlack() || fl.IsBlack()) continue;
00449 Ray r(pc, pl - pc, 1e-3f, .999f, time);
00450 if (!scene->IntersectP(r)) {
00451
00452 float pathWt = 1.f / (i + j + 2 - nSpecularVertices[i+j+2]);
00453 float G = AbsDot(nc, w) * AbsDot(nl, w) / DistanceSquared(pl, pc);
00454 L += (vc.alpha * fc * G * fl * vl.alpha) * pathWt;
00455 }
00456 }
00457 }
00458 }
00459 }
00460
00461 if (!eAlpha.IsBlack() && previousSpecular &&
00462 (directLighting == NULL || !allSpecular))
00463 for (uint32_t i = 0; i < scene->lights.size(); ++i)
00464 L += eAlpha * scene->lights[i]->Le(eRay);
00465 PBRT_MLT_FINISHED_LBIDIR();
00466 return L;
00467 }
00468
00469
00470 MetropolisRenderer::MetropolisRenderer(int perPixelSamples,
00471 int nboot, int dps, float lsp, bool dds, int mr, int md,
00472 Camera *c, bool db) {
00473 camera = c;
00474
00475 nPixelSamples = perPixelSamples;
00476 float largeStepProbability = lsp;
00477 largeStepsPerPixel = max(1u, RoundUpPow2(largeStepProbability * nPixelSamples));
00478 if (largeStepsPerPixel >= nPixelSamples) largeStepsPerPixel /= 2;
00479 Assert(largeStepsPerPixel >= 1 && largeStepsPerPixel < nPixelSamples);
00480 if ((nPixelSamples % largeStepsPerPixel) != 0) {
00481 int origPixelSamples = nPixelSamples;
00482 nPixelSamples += largeStepsPerPixel - (nPixelSamples % largeStepsPerPixel);
00483 Warning("Rounding up to %d Metropolis samples per pixel (from %d)",
00484 nPixelSamples, origPixelSamples);
00485 }
00486
00487 nBootstrap = nboot;
00488 nDirectPixelSamples = dps;
00489
00490 maxDepth = md;
00491 maxConsecutiveRejects = mr;
00492 nTasksFinished = 0;
00493 directLighting = dds ? new DirectLightingIntegrator(SAMPLE_ALL_UNIFORM, maxDepth) : NULL;
00494 bidirectional = db;
00495 }
00496
00497
00498 MetropolisRenderer::~MetropolisRenderer() {
00499 delete camera;
00500 delete directLighting;
00501 }
00502
00503
00504 MetropolisRenderer *CreateMetropolisRenderer(const ParamSet ¶ms,
00505 Camera *camera) {
00506 float largeStepProbability = params.FindOneFloat("largestepprobability", .25f);
00507 int perPixelSamples = params.FindOneInt("samplesperpixel", 100);
00508 int nBootstrap = params.FindOneInt("bootstrapsamples", 100000);
00509 int nDirectPixelSamples = params.FindOneInt("directsamples", 4);
00510 bool doDirectSeparately = params.FindOneBool("dodirectseparately", true);
00511 int mr = params.FindOneInt("maxconsecutiverejects", 512);
00512 int md = params.FindOneInt("maxdepth", 7);
00513 bool doBidirectional = params.FindOneBool("bidirectional", true);
00514
00515 if (PbrtOptions.quickRender) {
00516 perPixelSamples = max(1, perPixelSamples / 4);
00517 nBootstrap = max(1, nBootstrap / 4);
00518 nDirectPixelSamples = max(1, nDirectPixelSamples / 4);
00519 }
00520
00521 return new MetropolisRenderer(perPixelSamples, nBootstrap,
00522 nDirectPixelSamples, largeStepProbability, doDirectSeparately,
00523 mr, md, camera, doBidirectional);
00524 }
00525
00526
00527 void MetropolisRenderer::Render(const Scene *scene) {
00528 PBRT_MLT_STARTED_RENDERING();
00529 if (scene->lights.size() > 0) {
00530 int x0, x1, y0, y1;
00531 camera->film->GetPixelExtent(&x0, &x1, &y0, &y1);
00532 float t0 = camera->shutterOpen, t1 = camera->shutterClose;
00533 Distribution1D *lightDistribution = ComputeLightSamplingCDF(scene);
00534
00535 if (directLighting != NULL) {
00536 PBRT_MLT_STARTED_DIRECTLIGHTING();
00537
00538 if (nDirectPixelSamples > 0) {
00539 LDSampler sampler(x0, x1, y0, y1, nDirectPixelSamples, t0, t1);
00540 Sample *sample = new Sample(&sampler, directLighting, NULL, scene);
00541 vector<Task *> directTasks;
00542 int nDirectTasks = max(32 * NumSystemCores(),
00543 (camera->film->xResolution * camera->film->yResolution) / (16*16));
00544 nDirectTasks = RoundUpPow2(nDirectTasks);
00545 ProgressReporter directProgress(nDirectTasks, "Direct Lighting");
00546 for (int i = 0; i < nDirectTasks; ++i)
00547 directTasks.push_back(new SamplerRendererTask(scene, this, camera, directProgress,
00548 &sampler, sample, i, nDirectTasks));
00549 std::reverse(directTasks.begin(), directTasks.end());
00550 EnqueueTasks(directTasks);
00551 WaitForAllTasks();
00552 for (uint32_t i = 0; i < directTasks.size(); ++i)
00553 delete directTasks[i];
00554 delete sample;
00555 directProgress.Done();
00556 }
00557 camera->film->WriteImage();
00558 PBRT_MLT_FINISHED_DIRECTLIGHTING();
00559 }
00560
00561 PBRT_MLT_STARTED_BOOTSTRAPPING(nBootstrap);
00562 RNG rng(0);
00563 MemoryArena arena;
00564 vector<float> bootstrapI;
00565 vector<PathVertex> cameraPath(maxDepth, PathVertex());
00566 vector<PathVertex> lightPath(maxDepth, PathVertex());
00567 float sumI = 0.f;
00568 bootstrapI.reserve(nBootstrap);
00569 MLTSample sample(maxDepth);
00570 for (uint32_t i = 0; i < nBootstrap; ++i) {
00571
00572 float x = Lerp(rng.RandomFloat(), x0, x1);
00573 float y = Lerp(rng.RandomFloat(), y0, y1);
00574 LargeStep(rng, &sample, maxDepth, x, y, t0, t1, bidirectional);
00575 Spectrum L = PathL(sample, scene, arena, camera, lightDistribution,
00576 &cameraPath[0], &lightPath[0], rng);
00577
00578
00579 float I = ::I(L);
00580 sumI += I;
00581 bootstrapI.push_back(I);
00582 arena.FreeAll();
00583 }
00584 float b = sumI / nBootstrap;
00585 PBRT_MLT_FINISHED_BOOTSTRAPPING(b);
00586 Info("MLT computed b = %f", b);
00587
00588
00589 float contribOffset = rng.RandomFloat() * sumI;
00590 rng.Seed(0);
00591 sumI = 0.f;
00592 MLTSample initialSample(maxDepth);
00593 for (uint32_t i = 0; i < nBootstrap; ++i) {
00594 float x = Lerp(rng.RandomFloat(), x0, x1);
00595 float y = Lerp(rng.RandomFloat(), y0, y1);
00596 LargeStep(rng, &initialSample, maxDepth, x, y, t0, t1,
00597 bidirectional);
00598 sumI += bootstrapI[i];
00599 if (sumI > contribOffset)
00600 break;
00601 }
00602
00603
00604 uint32_t nTasks = largeStepsPerPixel;
00605 uint32_t largeStepRate = nPixelSamples / largeStepsPerPixel;
00606 Info("MLT running %d tasks, large step rate %d", nTasks, largeStepRate);
00607 ProgressReporter progress(nTasks * largeStepRate, "Metropolis");
00608 vector<Task *> tasks;
00609 Mutex *filmMutex = Mutex::Create();
00610 Assert(IsPowerOf2(nTasks));
00611 uint32_t scramble[2] = { rng.RandomUInt(), rng.RandomUInt() };
00612 uint32_t pfreq = (x1-x0) * (y1-y0);
00613 for (uint32_t i = 0; i < nTasks; ++i) {
00614 float d[2];
00615 Sample02(i, scramble, d);
00616 tasks.push_back(new MLTTask(progress, pfreq, i,
00617 d[0], d[1], x0, x1, y0, y1, t0, t1, b, initialSample,
00618 scene, camera, this, filmMutex, lightDistribution));
00619 }
00620 EnqueueTasks(tasks);
00621 WaitForAllTasks();
00622 for (uint32_t i = 0; i < tasks.size(); ++i)
00623 delete tasks[i];
00624 progress.Done();
00625 Mutex::Destroy(filmMutex);
00626 delete lightDistribution;
00627 }
00628 camera->film->WriteImage();
00629 PBRT_MLT_FINISHED_RENDERING();
00630 }
00631
00632
00633 inline float I(const Spectrum &L) {
00634 return L.y();
00635 }
00636
00637
00638 MLTTask::MLTTask(ProgressReporter &prog, uint32_t pfreq, uint32_t tn,
00639 float ddx, float ddy, int xx0, int xx1, int yy0, int yy1, float tt0, float tt1,
00640 float bb, const MLTSample &is, const Scene *sc, const Camera *c,
00641 MetropolisRenderer *ren, Mutex *fm, Distribution1D *ld)
00642 : progress(prog), initialSample(is) {
00643 progressUpdateFrequency = pfreq;
00644 taskNum = tn;
00645 dx = ddx;
00646 dy = ddy;
00647 x0 = xx0;
00648 x1 = xx1;
00649 y0 = yy0;
00650 y1 = yy1;
00651 t0 = tt0;
00652 t1 = tt1;
00653 currentPixelSample = 0;
00654 b = bb;
00655 scene = sc;
00656 camera = c;
00657 renderer = ren;
00658 filmMutex = fm;
00659 lightDistribution = ld;
00660 }
00661
00662
00663 void MLTTask::Run() {
00664 PBRT_MLT_STARTED_MLT_TASK(this);
00665
00666 PBRT_MLT_STARTED_TASK_INIT();
00667 uint32_t nPixels = (x1-x0) * (y1-y0);
00668 uint32_t nPixelSamples = renderer->nPixelSamples;
00669 uint32_t largeStepRate = nPixelSamples / renderer->largeStepsPerPixel;
00670 Assert(largeStepRate > 1);
00671 uint64_t nTaskSamples = uint64_t(nPixels) * uint64_t(largeStepRate);
00672 uint32_t consecutiveRejects = 0;
00673 uint32_t progressCounter = progressUpdateFrequency;
00674
00675
00676 MemoryArena arena;
00677 RNG rng(taskNum);
00678 vector<PathVertex> cameraPath(renderer->maxDepth, PathVertex());
00679 vector<PathVertex> lightPath(renderer->maxDepth, PathVertex());
00680 vector<MLTSample> samples(2, MLTSample(renderer->maxDepth));
00681 Spectrum L[2];
00682 float I[2];
00683 uint32_t current = 0, proposed = 1;
00684
00685
00686 samples[current] = initialSample;
00687 L[current] = renderer->PathL(initialSample, scene, arena, camera,
00688 lightDistribution, &cameraPath[0], &lightPath[0], rng);
00689 I[current] = ::I(L[current]);
00690 arena.FreeAll();
00691
00692
00693 uint32_t pixelNumOffset = 0;
00694 vector<int> largeStepPixelNum;
00695 largeStepPixelNum.reserve(nPixels);
00696 for (uint32_t i = 0; i < nPixels; ++i) largeStepPixelNum.push_back(i);
00697 Shuffle(&largeStepPixelNum[0], nPixels, 1, rng);
00698 PBRT_MLT_FINISHED_TASK_INIT();
00699 for (uint64_t s = 0; s < nTaskSamples; ++s) {
00700
00701 PBRT_MLT_STARTED_MUTATION();
00702 samples[proposed] = samples[current];
00703 bool largeStep = ((s % largeStepRate) == 0);
00704 if (largeStep) {
00705 int x = x0 + largeStepPixelNum[pixelNumOffset] % (x1 - x0);
00706 int y = y0 + largeStepPixelNum[pixelNumOffset] / (x1 - x0);
00707 LargeStep(rng, &samples[proposed], renderer->maxDepth,
00708 x + dx, y + dy, t0, t1, renderer->bidirectional);
00709 ++pixelNumOffset;
00710 }
00711 else
00712 SmallStep(rng, &samples[proposed], renderer->maxDepth,
00713 x0, x1, y0, y1, t0, t1, renderer->bidirectional);
00714 PBRT_MLT_FINISHED_MUTATION();
00715
00716
00717 L[proposed] = renderer->PathL(samples[proposed], scene, arena, camera,
00718 lightDistribution, &cameraPath[0], &lightPath[0], rng);
00719 I[proposed] = ::I(L[proposed]);
00720 arena.FreeAll();
00721
00722
00723 float a = min(1.f, I[proposed] / I[current]);
00724
00725
00726 PBRT_MLT_STARTED_SAMPLE_SPLAT();
00727 if (I[current] > 0.f) {
00728 if (!isinf(1.f / I[current])) {
00729 Spectrum contrib = (b / nPixelSamples) * L[current] / I[current];
00730 camera->film->Splat(samples[current].cameraSample,
00731 (1.f - a) * contrib);
00732 }
00733 }
00734 if (I[proposed] > 0.f) {
00735 if (!isinf(1.f / I[proposed])) {
00736 Spectrum contrib = (b / nPixelSamples) * L[proposed] / I[proposed];
00737 camera->film->Splat(samples[proposed].cameraSample,
00738 a * contrib);
00739 }
00740 }
00741 PBRT_MLT_FINISHED_SAMPLE_SPLAT();
00742
00743
00744 if (consecutiveRejects >= renderer->maxConsecutiveRejects ||
00745 rng.RandomFloat() < a) {
00746 PBRT_MLT_ACCEPTED_MUTATION(a, &samples[current], &samples[proposed]);
00747 current ^= 1;
00748 proposed ^= 1;
00749 consecutiveRejects = 0;
00750 }
00751 else
00752 {
00753 PBRT_MLT_REJECTED_MUTATION(a, &samples[current], &samples[proposed]);
00754 ++consecutiveRejects;
00755 }
00756 if (--progressCounter == 0) {
00757 progress.Update();
00758 progressCounter = progressUpdateFrequency;
00759 }
00760 }
00761 Assert(pixelNumOffset == nPixels);
00762
00763 PBRT_MLT_STARTED_DISPLAY_UPDATE();
00764 int ntf = AtomicAdd(&renderer->nTasksFinished, 1);
00765 int64_t totalSamples = int64_t(nPixels) * int64_t(nPixelSamples);
00766 float splatScale = float(double(totalSamples) / double(ntf * nTaskSamples));
00767 camera->film->UpdateDisplay(x0, y0, x1, y1, splatScale);
00768 if ((taskNum % 8) == 0) {
00769 MutexLock lock(*filmMutex);
00770 camera->film->WriteImage(splatScale);
00771 }
00772 PBRT_MLT_FINISHED_DISPLAY_UPDATE();
00773 PBRT_MLT_FINISHED_MLT_TASK(this);
00774 }
00775
00776
00777 Spectrum MetropolisRenderer::Li(const Scene *scene, const RayDifferential &ray,
00778 const Sample *sample, RNG &rng, MemoryArena &arena, Intersection *isect,
00779 Spectrum *T) const {
00780 Spectrum localT;
00781 if (!T) T = &localT;
00782 Intersection localIsect;
00783 if (!isect) isect = &localIsect;
00784 Spectrum Lo = 0.f;
00785 if (scene->Intersect(ray, isect))
00786 Lo = directLighting->Li(scene, this, ray, *isect, sample,
00787 rng, arena);
00788 else {
00789
00790 for (uint32_t i = 0; i < scene->lights.size(); ++i)
00791 Lo += scene->lights[i]->Le(ray);
00792 }
00793 return Lo;
00794 }
00795
00796
00797 Spectrum MetropolisRenderer::Transmittance(const Scene *scene, const RayDifferential &ray,
00798 const Sample *sample, RNG &rng, MemoryArena &arena) const {
00799 return 1.f;
00800 }
00801
00802