C++ MemoryArena::FreeAll方法代码示例

void IrradiancePrimeTask::Run() {
    if (!sampler) { progress.Update(); return; }
    MemoryArena arena;
    int sampleCount;
    RNG rng(29 * taskNum);
    int maxSamples = sampler->MaximumSampleCount();
    Sample *samples = origSample->Duplicate(maxSamples);
    while ((sampleCount = sampler->GetMoreSamples(samples, rng)) > 0) {
        for (int i = 0; i < sampleCount; ++i) {
            RayDifferential ray;
            camera->GenerateRayDifferential(samples[i], &ray);
            Intersection isect;
            if (scene->Intersect(ray, &isect))
                (void)irradianceCache->Li(scene, renderer, ray, isect, &samples[i], rng, arena);
        }
        arena.FreeAll();
    }
    delete[] samples;
    delete sampler;
    progress.Update();
}

void IGIIntegrator::Preprocess(const Scene *scene, const Camera *camera,
                               const Renderer *renderer) {
    if (scene->lights.size() == 0) return;
    MemoryArena arena;
    RNG rng;
    // Compute samples for emitted rays from lights
    vector<float> lightNum(nLightPaths * nLightSets);
    vector<float> lightSampPos(2 * nLightPaths * nLightSets, 0.f);
    vector<float> lightSampComp(nLightPaths * nLightSets, 0.f);
    vector<float> lightSampDir(2 * nLightPaths * nLightSets, 0.f);
    LDShuffleScrambled1D(nLightPaths, nLightSets, &lightNum[0], rng);
    LDShuffleScrambled2D(nLightPaths, nLightSets, &lightSampPos[0], rng);
    LDShuffleScrambled1D(nLightPaths, nLightSets, &lightSampComp[0], rng);
    LDShuffleScrambled2D(nLightPaths, nLightSets, &lightSampDir[0], rng);

    // Precompute information for light sampling densities
    Distribution1D *lightDistribution = ComputeLightSamplingCDF(scene);
    for (uint32_t s = 0; s < nLightSets; ++s) {
        for (uint32_t i = 0; i < nLightPaths; ++i) {
            // Follow path _i_ from light to create virtual lights
            int sampOffset = s*nLightPaths + i;

            // Choose light source to trace virtual light path from
            float lightPdf;
            int ln = lightDistribution->SampleDiscrete(lightNum[sampOffset],
                                                       &lightPdf);
            Light *light = scene->lights[ln];

            // Sample ray leaving light source for virtual light path
            RayDifferential ray;
            float pdf;
            LightSample ls(lightSampPos[2*sampOffset], lightSampPos[2*sampOffset+1],
                           lightSampComp[sampOffset]);
            Normal Nl;
            Spectrum alpha = light->Sample_L(scene, ls, lightSampDir[2*sampOffset],
                                             lightSampDir[2*sampOffset+1],
                                             camera->shutterOpen, &ray, &Nl, &pdf);
            if (pdf == 0.f || alpha.IsBlack()) continue;
            alpha /= pdf * lightPdf;
            Intersection isect;
            while (scene->Intersect(ray, &isect) && !alpha.IsBlack()) {
                // Create virtual light and sample new ray for path
                alpha *= renderer->Transmittance(scene, RayDifferential(ray), NULL,
                                                 rng, arena);
                Vector wo = -ray.d;
                BSDF *bsdf = isect.GetBSDF(ray, arena);

                // Create virtual light at ray intersection point
                Spectrum contrib = alpha * bsdf->rho(wo, rng) / M_PI;
                virtualLights[s].push_back(VirtualLight(isect.dg.p, isect.dg.nn, contrib,
                                                        isect.rayEpsilon));

                // Sample new ray direction and update weight for virtual light path
                Vector wi;
                float pdf;
                BSDFSample bsdfSample(rng);
                Spectrum fr = bsdf->Sample_f(wo, &wi, bsdfSample, &pdf);
                if (fr.IsBlack() || pdf == 0.f)
                    break;
                Spectrum contribScale = fr * AbsDot(wi, bsdf->dgShading.nn) / pdf;

                // Possibly terminate virtual light path with Russian roulette
                float rrProb = min(1.f, contribScale.y());
                if (rng.RandomFloat() > rrProb)
                    break;
                alpha *= contribScale / rrProb;
                ray = RayDifferential(isect.dg.p, wi, ray, isect.rayEpsilon);
            }
            arena.FreeAll();
        }
    }
    delete lightDistribution;
}

// SamplerRendererTask Definitions
void SamplerRendererTask::Run() {
    PBRT_STARTED_RENDERTASK(taskNum);
    // Get sub-_Sampler_ for _SamplerRendererTask_
    Sampler *sampler = mainSampler->GetSubSampler(taskNum, taskCount);
    if (!sampler)
    {
        reporter.Update();
        PBRT_FINISHED_RENDERTASK(taskNum);
        return;
    }

    // Declare local variables used for rendering loop
    MemoryArena arena;
    RNG rng(taskNum);

    // Allocate space for samples and intersections
    int maxSamples = sampler->MaximumSampleCount();
    Sample *samples = origSample->Duplicate(maxSamples);
    RayDifferential *rays = new RayDifferential[maxSamples];
    Spectrum *Ls = new Spectrum[maxSamples];
    Spectrum *Ts = new Spectrum[maxSamples];
    Intersection *isects = new Intersection[maxSamples];

    // Get samples from \use{Sampler} and update image
    int sampleCount;
    while ((sampleCount = sampler->GetMoreSamples(samples, rng)) > 0) {
        // Generate camera rays and compute radiance along rays
        for (int i = 0; i < sampleCount; ++i) {
            // Find camera ray for _sample[i]_
            PBRT_STARTED_GENERATING_CAMERA_RAY(&samples[i]);
            float rayWeight = camera->GenerateRayDifferential(samples[i], &rays[i]);
            PBRT_FINISHED_GENERATING_CAMERA_RAY(&samples[i], &rays[i], rayWeight);

            // Evaluate radiance along camera ray
            PBRT_STARTED_CAMERA_RAY_INTEGRATION(&rays[i], &samples[i]);
            if (rayWeight > 0.f)
                Ls[i] = rayWeight * renderer->Li(scene, rays[i], &samples[i], rng,
                                                 arena, &isects[i], &Ts[i]);
            else {
                Ls[i] = 0.f;
                Ts[i] = 1.f;
            }

            // Issue warning if unexpected radiance value returned
            if (Ls[i].HasNaNs()) {
                Error("Not-a-number radiance value returned "
                      "for image sample.  Setting to black.");
                Ls[i] = Spectrum(0.f);
            }
            else if (Ls[i].y() < -1e-5) {
                Error("Negative luminance value, %f, returned"
                      "for image sample.  Setting to black.", Ls[i].y());
                Ls[i] = Spectrum(0.f);
            }
            else if (isinf(Ls[i].y())) {
                Error("Infinite luminance value returned"
                      "for image sample.  Setting to black.");
                Ls[i] = Spectrum(0.f);
            }
            PBRT_FINISHED_CAMERA_RAY_INTEGRATION(&rays[i], &samples[i], &Ls[i]);
        }

        // Report sample results to _Sampler_, add contributions to image
        if (sampler->ReportResults(samples, rays, Ls, isects, sampleCount))
        {
            for (int i = 0; i < sampleCount; ++i)
            {
                PBRT_STARTED_ADDING_IMAGE_SAMPLE(&samples[i], &rays[i], &Ls[i], &Ts[i]);
                camera->film->AddSample(samples[i], Ls[i]);
                PBRT_FINISHED_ADDING_IMAGE_SAMPLE();
            }
        }

        // Free \use{MemoryArena} memory from computing image sample values
        arena.FreeAll();
    }

    // Clean up after \use{SamplerRendererTask} is done with its image region
    camera->film->UpdateDisplay(sampler->xPixelStart,
        sampler->yPixelStart, sampler->xPixelEnd+1, sampler->yPixelEnd+1);
    delete sampler;
    delete[] samples;
    delete[] rays;
    delete[] Ls;
    delete[] Ts;
    delete[] isects;
    reporter.Update();
    PBRT_FINISHED_RENDERTASK(taskNum);
}

void DipoleSubsurfaceIntegrator::Preprocess(const Scene *scene,
        const Camera *camera, const Renderer *renderer) {
    if (scene->lights.size() == 0) return;
    vector<SurfacePoint> pts;
    // Get _SurfacePoint_s for translucent objects in scene
    if (filename != "") {
        // Initialize _SurfacePoint_s from file
        vector<float> fpts;
        if (ReadFloatFile(filename.c_str(), &fpts)) {
            if ((fpts.size() % 8) != 0)
                Error("Excess values (%d) in points file \"%s\"", int(fpts.size() % 8),
                      filename.c_str());
            for (u_int i = 0; i < fpts.size(); i += 8)
                pts.push_back(SurfacePoint(Point(fpts[i], fpts[i+1], fpts[i+2]),
                                           Normal(fpts[i+3], fpts[i+4], fpts[i+5]),
                                           fpts[i+6], fpts[i+7]));
        }
    }
    if (pts.size() == 0) {
        Point pCamera = camera->CameraToWorld(camera->shutterOpen,
                                              Point(0, 0, 0));
        FindPoissonPointDistribution(pCamera, camera->shutterOpen,
                                     minSampleDist, scene, &pts);
    }

    // Compute irradiance values at sample points
    RNG rng;
    MemoryArena arena;
    PBRT_SUBSURFACE_STARTED_COMPUTING_IRRADIANCE_VALUES();
    ProgressReporter progress(pts.size(), "Computing Irradiances");
    for (uint32_t i = 0; i < pts.size(); ++i) {
        SurfacePoint &sp = pts[i];
        Spectrum E(0.f);
        for (uint32_t j = 0; j < scene->lights.size(); ++j) {
            // Add irradiance from light at point
            const Light *light = scene->lights[j];
            Spectrum Elight = 0.f;
            int nSamples = RoundUpPow2(light->nSamples);
            uint32_t scramble[2] = { rng.RandomUInt(), rng.RandomUInt() };
            uint32_t compScramble = rng.RandomUInt();
            for (int s = 0; s < nSamples; ++s) {
                float lpos[2];
                Sample02(s, scramble, lpos);
                float lcomp = VanDerCorput(s, compScramble);
                LightSample ls(lpos[0], lpos[1], lcomp);
                Vector wi;
                float lightPdf;
                VisibilityTester visibility;
                Spectrum Li = light->Sample_L(sp.p, sp.rayEpsilon,
                    ls, camera->shutterOpen, &wi, &lightPdf, &visibility);
                if (Dot(wi, sp.n) <= 0.) continue;
                if (Li.IsBlack() || lightPdf == 0.f) continue;
                Li *= visibility.Transmittance(scene, renderer, NULL, rng, arena);
                if (visibility.Unoccluded(scene))
                    Elight += Li * AbsDot(wi, sp.n) / lightPdf;
            }
            E += Elight / nSamples;
        }
        if (E.y() > 0.f)
        {
            irradiancePoints.push_back(IrradiancePoint(sp, E));
            PBRT_SUBSURFACE_COMPUTED_IRRADIANCE_AT_POINT(&sp, &E);
        }
        arena.FreeAll();
        progress.Update();
    }
    progress.Done();
    PBRT_SUBSURFACE_FINISHED_COMPUTING_IRRADIANCE_VALUES();

    // Create octree of clustered irradiance samples
    octree = octreeArena.Alloc<SubsurfaceOctreeNode>();
    for (uint32_t i = 0; i < irradiancePoints.size(); ++i)
        octreeBounds = Union(octreeBounds, irradiancePoints[i].p);
    for (uint32_t i = 0; i < irradiancePoints.size(); ++i)
        octree->Insert(octreeBounds, &irradiancePoints[i], octreeArena);
    octree->InitHierarchy();
}

void MLTTask::Run() {
    PBRT_MLT_STARTED_MLT_TASK(this);
    // Declare basic _MLTTask_ variables and prepare for sampling
    PBRT_MLT_STARTED_TASK_INIT();
    uint32_t nPixels = (x1-x0) * (y1-y0);
    uint32_t nPixelSamples = renderer->nPixelSamples;
    uint32_t largeStepRate = nPixelSamples / renderer->largeStepsPerPixel;
    Assert(largeStepRate > 1);
    uint64_t nTaskSamples = uint64_t(nPixels) * uint64_t(largeStepRate);
    uint32_t consecutiveRejects = 0;
    uint32_t progressCounter = progressUpdateFrequency;

    // Declare variables for storing and computing MLT samples
    MemoryArena arena;
    RNG rng(taskNum);
    vector<PathVertex> cameraPath(renderer->maxDepth, PathVertex());
    vector<PathVertex> lightPath(renderer->maxDepth, PathVertex());
    vector<MLTSample> samples(2, MLTSample(renderer->maxDepth));
    Spectrum L[2];
    float I[2];
    uint32_t current = 0, proposed = 1;

    // Compute _L[current]_ for initial sample
    samples[current] = initialSample;
    L[current] = renderer->PathL(initialSample, scene, arena, camera,
                     lightDistribution, &cameraPath[0], &lightPath[0], rng);
    I[current] = ::I(L[current]);
    arena.FreeAll();

    // Compute randomly permuted table of pixel indices for large steps
    uint32_t pixelNumOffset = 0;
    vector<int> largeStepPixelNum;
    largeStepPixelNum.reserve(nPixels);
    for (uint32_t i = 0; i < nPixels; ++i) largeStepPixelNum.push_back(i);
    Shuffle(&largeStepPixelNum[0], nPixels, 1, rng);
    PBRT_MLT_FINISHED_TASK_INIT();
    for (uint64_t s = 0; s < nTaskSamples; ++s) {
        // Compute proposed mutation to current sample
        PBRT_MLT_STARTED_MUTATION();
        samples[proposed] = samples[current];
        bool largeStep = ((s % largeStepRate) == 0);
        if (largeStep) {
            int x = x0 + largeStepPixelNum[pixelNumOffset] % (x1 - x0);
            int y = y0 + largeStepPixelNum[pixelNumOffset] / (x1 - x0);
            LargeStep(rng, &samples[proposed], renderer->maxDepth,
                      x + dx, y + dy, t0, t1, renderer->bidirectional);
            ++pixelNumOffset;
        }
        else
            SmallStep(rng, &samples[proposed], renderer->maxDepth,
                      x0, x1, y0, y1, t0, t1, renderer->bidirectional);
        PBRT_MLT_FINISHED_MUTATION();

        // Compute contribution of proposed sample
        L[proposed] = renderer->PathL(samples[proposed], scene, arena, camera,
                         lightDistribution, &cameraPath[0], &lightPath[0], rng);
        I[proposed] = ::I(L[proposed]);
        arena.FreeAll();

        // Compute acceptance probability for proposed sample
        float a = min(1.f, I[proposed] / I[current]);

        // Splat current and proposed samples to _Film_
        PBRT_MLT_STARTED_SAMPLE_SPLAT();
        if (I[current] > 0.f) {
            if (!isinf(1.f / I[current])) {
            Spectrum contrib =  (b / nPixelSamples) * L[current] / I[current];
            camera->film->Splat(samples[current].cameraSample,
                                (1.f - a) * contrib);
        }
        }
        if (I[proposed] > 0.f) {
            if (!isinf(1.f / I[proposed])) {
            Spectrum contrib =  (b / nPixelSamples) * L[proposed] / I[proposed];
            camera->film->Splat(samples[proposed].cameraSample,
                                a * contrib);
        }
        }
        PBRT_MLT_FINISHED_SAMPLE_SPLAT();

        // Randomly accept proposed path mutation (or not)
        if (consecutiveRejects >= renderer->maxConsecutiveRejects ||
            rng.RandomFloat() < a) {
            PBRT_MLT_ACCEPTED_MUTATION(a, &samples[current], &samples[proposed]);
            current ^= 1;
            proposed ^= 1;
            consecutiveRejects = 0;
        }
        else
        {
            PBRT_MLT_REJECTED_MUTATION(a, &samples[current], &samples[proposed]);
            ++consecutiveRejects;
        }
        if (--progressCounter == 0) {
            progress.Update();
            progressCounter = progressUpdateFrequency;
        }
    }
    Assert(pixelNumOffset == nPixels);
    // Update display for recently computed Metropolis samples
//.........这里部分代码省略.........

void CreateRadProbeTask::Run() {
    // Compute region in which to compute incident radiance probes
    int sx = pointNum % nProbes[0];
    int sy = (pointNum / nProbes[0]) % nProbes[1];
    int sz = pointNum / (nProbes[0] * nProbes[1]);
    Assert(sx >= 0 && sx < nProbes[0]);
    Assert(sy >= 0 && sy < nProbes[1]);
    Assert(sz >= 0 && sz < nProbes[2]);
    float tx0 = float(sx) / nProbes[0], tx1 = float(sx+1) / nProbes[0];
    float ty0 = float(sy) / nProbes[1], ty1 = float(sy+1) / nProbes[1];
    float tz0 = float(sz) / nProbes[2], tz1 = float(sz+1) / nProbes[2];
    BBox b(bbox.Lerp(tx0, ty0, tz0), bbox.Lerp(tx1, ty1, tz1));

    // Initialize common variables for _CreateRadProbeTask::Run()_
    RNG rng(pointNum);
    Spectrum *c_probe = new Spectrum[SHTerms(lmax)];
    MemoryArena arena;
    uint32_t nFound = 0, lastVisibleOffset = 0;
    for (int i = 0; i < 256; ++i) {
        if (nFound == 32) break;
        // Try to compute radiance probe contribution at _i_th sample point

        // Compute _i_th candidate point _p_ in cell's bounding box
        float dx = RadicalInverse(i+1, 2);
        float dy = RadicalInverse(i+1, 3);
        float dz = RadicalInverse(i+1, 5);
        Point p = b.Lerp(dx, dy, dz);

        // Skip point _p_ if not indirectly visible from camera
        if (scene->IntersectP(Ray(surfacePoints[lastVisibleOffset],
                                  p - surfacePoints[lastVisibleOffset],
                                  1e-4f, 1.f, time))) {
            uint32_t j;
            // See if point is visible to any element of _surfacePoints_
            for (j = 0; j < surfacePoints.size(); ++j)
                if (!scene->IntersectP(Ray(surfacePoints[j], p - surfacePoints[j],
                                           1e-4f, 1.f, time))) {
                    lastVisibleOffset = j;
                    break;
                }
            if (j == surfacePoints.size())
                continue;
        }
        ++nFound;

        // Compute SH coefficients of incident radiance at point _p_
        if (includeDirectInProbes) {
            for (int i = 0; i < SHTerms(lmax); ++i)
                c_probe[i] = 0.f;
            SHProjectIncidentDirectRadiance(p, 0.f, time, arena, scene,
                                            true, lmax, rng, c_probe);
            for (int i = 0; i < SHTerms(lmax); ++i)
                c_in[i] += c_probe[i];
        }

        if (includeIndirectInProbes) {
            for (int i = 0; i < SHTerms(lmax); ++i)
                c_probe[i] = 0.f;
            SHProjectIncidentIndirectRadiance(p, 0.f, time, renderer,
                origSample, scene, lmax, rng, nIndirSamples, c_probe);
            for (int i = 0; i < SHTerms(lmax); ++i)
                c_in[i] += c_probe[i];
        }
        arena.FreeAll();
    }
    // Compute final average value for probe and cleanup
    if (nFound > 0)
        for (int i = 0; i < SHTerms(lmax); ++i)
            c_in[i] /= nFound;
    delete[] c_probe;
        prog.Update();
}

void MetropolisRenderer::Render(const Scene *scene) {
    PBRT_MLT_STARTED_RENDERING();
    if (scene->lights.size() > 0) {
        int x0, x1, y0, y1;
        camera->film->GetPixelExtent(&x0, &x1, &y0, &y1);
        float t0 = camera->shutterOpen, t1 = camera->shutterClose;
        Distribution1D *lightDistribution = ComputeLightSamplingCDF(scene);

        if (directLighting != NULL) {
            PBRT_MLT_STARTED_DIRECTLIGHTING();
            // Compute direct lighting before Metropolis light transport
            if (nDirectPixelSamples > 0) {
                LDSampler sampler(x0, x1, y0, y1, nDirectPixelSamples, t0, t1);
                Sample *sample = new Sample(&sampler, directLighting, NULL, scene);
                vector<Task *> directTasks;
                int nDirectTasks = max(32 * NumSystemCores(),
                                 (camera->film->xResolution * camera->film->yResolution) / (16*16));
                nDirectTasks = RoundUpPow2(nDirectTasks);
                ProgressReporter directProgress(nDirectTasks, "Direct Lighting");
                for (int i = 0; i < nDirectTasks; ++i)
                    directTasks.push_back(new SamplerRendererTask(scene, this, camera, directProgress,
                                                                  &sampler, sample, false, i, nDirectTasks));
                std::reverse(directTasks.begin(), directTasks.end());
                EnqueueTasks(directTasks);
                WaitForAllTasks();
                for (uint32_t i = 0; i < directTasks.size(); ++i)
                    delete directTasks[i];
                delete sample;
                directProgress.Done();
            }
            camera->film->WriteImage();
            PBRT_MLT_FINISHED_DIRECTLIGHTING();
        }
        // Take initial set of samples to compute $b$
        PBRT_MLT_STARTED_BOOTSTRAPPING(nBootstrap);
        RNG rng(0);
        MemoryArena arena;
        vector<float> bootstrapI;
        vector<PathVertex> cameraPath(maxDepth, PathVertex());
        vector<PathVertex> lightPath(maxDepth, PathVertex());
        float sumI = 0.f;
        bootstrapI.reserve(nBootstrap);
        MLTSample sample(maxDepth);
        for (uint32_t i = 0; i < nBootstrap; ++i) {
            // Generate random sample and path radiance for MLT bootstrapping
            float x = Lerp(rng.RandomFloat(), x0, x1);
            float y = Lerp(rng.RandomFloat(), y0, y1);
            LargeStep(rng, &sample, maxDepth, x, y, t0, t1, bidirectional);
            Spectrum L = PathL(sample, scene, arena, camera, lightDistribution,
                               &cameraPath[0], &lightPath[0], rng);

            // Compute contribution for random sample for MLT bootstrapping
            float I = ::I(L);
            sumI += I;
            bootstrapI.push_back(I);
            arena.FreeAll();
        }
        float b = sumI / nBootstrap;
        PBRT_MLT_FINISHED_BOOTSTRAPPING(b);
        Info("MLT computed b = %f", b);

        // Select initial sample from bootstrap samples
        float contribOffset = rng.RandomFloat() * sumI;
        rng.Seed(0);
        sumI = 0.f;
        MLTSample initialSample(maxDepth);
        for (uint32_t i = 0; i < nBootstrap; ++i) {
            float x = Lerp(rng.RandomFloat(), x0, x1);
            float y = Lerp(rng.RandomFloat(), y0, y1);
            LargeStep(rng, &initialSample, maxDepth, x, y, t0, t1,
                      bidirectional);
            sumI += bootstrapI[i];
            if (sumI > contribOffset)
                break;
        }

        // Launch tasks to generate Metropolis samples
        uint32_t nTasks = largeStepsPerPixel;
        uint32_t largeStepRate = nPixelSamples / largeStepsPerPixel;
        Info("MLT running %d tasks, large step rate %d", nTasks, largeStepRate);
        ProgressReporter progress(nTasks * largeStepRate, "Metropolis");
        vector<Task *> tasks;
        Mutex *filmMutex = Mutex::Create();
        Assert(IsPowerOf2(nTasks));
        uint32_t scramble[2] = { rng.RandomUInt(), rng.RandomUInt() };
        uint32_t pfreq = (x1-x0) * (y1-y0);
        for (uint32_t i = 0; i < nTasks; ++i) {
            float d[2];
            Sample02(i, scramble, d);
            tasks.push_back(new MLTTask(progress, pfreq, i,
                d[0], d[1], x0, x1, y0, y1, t0, t1, b, initialSample,
                scene, camera, this, filmMutex, lightDistribution));
        }
        EnqueueTasks(tasks);
        WaitForAllTasks();
        for (uint32_t i = 0; i < tasks.size(); ++i)
            delete tasks[i];
        progress.Done();
        Mutex::Destroy(filmMutex);
        delete lightDistribution;
//.........这里部分代码省略.........

void SurfacePointTask::Run() {
    // Declare common variables for _SurfacePointTask::Run()_
    RNG rng(37 * taskNum);
    MemoryArena arena;
    vector<SurfacePoint> candidates;
    while (true) {
        int pathsTraced, raysTraced = 0;
        for (pathsTraced = 0; pathsTraced < 20000; ++pathsTraced) {
            // Follow ray path and attempt to deposit candidate sample points
            Vector dir = UniformSampleSphere(rng.RandomFloat(), rng.RandomFloat());
            Ray ray(origin, dir, 0.f, INFINITY, time);
            while (ray.depth < 30) {
                // Find ray intersection with scene geometry or bounding sphere
                ++raysTraced;
                bool hitOnSphere = false;
                auto optIsect = scene.Intersect(ray);
                if (!optIsect) {
                    optIsect = sphere.Intersect(ray);
                    if (!optIsect)
                        break;
                    hitOnSphere = true;
                }
                DifferentialGeometry &hitGeometry = optIsect->dg;
                hitGeometry.nn = Faceforward(hitGeometry.nn, -ray.d);

                // Store candidate sample point at ray intersection if appropriate
                if (!hitOnSphere && ray.depth >= 3 &&
                    optIsect->GetBSSRDF(RayDifferential(ray), arena) != NULL) {
                    float area = M_PI * (minSampleDist / 2.f) * (minSampleDist / 2.f);
                    candidates.push_back(SurfacePoint(hitGeometry.p, hitGeometry.nn,
                                                      area, optIsect->rayEpsilon));
                }

                // Generate random ray from intersection point
                Vector dir = UniformSampleSphere(rng.RandomFloat(), rng.RandomFloat());
                dir = Faceforward(dir, hitGeometry.nn);
                ray = Ray(hitGeometry.p, dir, ray, optIsect->rayEpsilon);
            }
            arena.FreeAll();
        }
        // Make first pass through candidate points with reader lock
        vector<bool> candidateRejected;
        candidateRejected.reserve(candidates.size());
        RWMutexLock lock(mutex, READ);
        for (uint32_t i = 0; i < candidates.size(); ++i) {
            PoissonCheck check(minSampleDist, candidates[i].p);
            octree.Lookup(candidates[i].p, check);
            candidateRejected.push_back(check.failed);
        }

        // Make second pass through points with writer lock and update octree
        lock.UpgradeToWrite();
        if (repeatedFails >= maxFails)
            return;
        totalPathsTraced += pathsTraced;
        totalRaysTraced += raysTraced;
        int oldMaxRepeatedFails = maxRepeatedFails;
        for (uint32_t i = 0; i < candidates.size(); ++i) {
            if (candidateRejected[i]) {
                // Update for rejected candidate point
                ++repeatedFails;
                maxRepeatedFails = max(maxRepeatedFails, repeatedFails);
                if (repeatedFails >= maxFails)
                    return;
            }
            else {
                // Recheck candidate point and possibly add to octree
                SurfacePoint &sp = candidates[i];
                PoissonCheck check(minSampleDist, sp.p);
                octree.Lookup(sp.p, check);
                if (check.failed) {
                    // Update for rejected candidate point
                    ++repeatedFails;
                    maxRepeatedFails = max(maxRepeatedFails, repeatedFails);
                    if (repeatedFails >= maxFails)
                        return;
                }
                else {
                    ++numPointsAdded;
                    repeatedFails = 0;
                    Vector delta(minSampleDist, minSampleDist, minSampleDist);
                    octree.Add(sp, BBox(sp.p-delta, sp.p+delta));
                    PBRT_SUBSURFACE_ADDED_POINT_TO_OCTREE(&sp, minSampleDist);
                    surfacePoints.push_back(sp);
                }
            }
        }

        // Stop following paths if not finding new points
        if (repeatedFails > oldMaxRepeatedFails) {
            int delta = repeatedFails - oldMaxRepeatedFails;
            prog.Update(delta);
        }
        if (totalPathsTraced > 50000 && numPointsAdded == 0) {
            Warning("There don't seem to be any objects with BSSRDFs "
                    "in this scene.  Giving up.");
            return;
        }
        candidates.erase(candidates.begin(), candidates.end());
    }
//.........这里部分代码省略.........

void MLTTask::Run() {
    PBRT_MLT_STARTED_MLT_TASK(this);
    // Declare basic _MLTTask_ variables and prepare for sampling
    RNG rng(taskNum);
    MemoryArena arena;
    vector<MLTSample> mltSamples(2, MLTSample(maxDepth));
    Spectrum sampleLs[2];
    uint32_t currentSample = 0, proposedSample = 1;
    mltSamples[currentSample] = initialSample;
    sampleLs[currentSample] = L(scene, renderer, camera, arena, rng, maxDepth,
                                ignoreDirect, mltSamples[currentSample]);
    int consecutiveRejects = 0;
    for (int sampleNum = 0; sampleNum < nSamples; ++sampleNum) {
        // Compute proposed mutation to current sample
        bool largeStep = rng.RandomFloat() < largeStepProbability;
        mltSamples[proposedSample] = mltSamples[currentSample];
        if (largeStep)
            LargeStep(rng, &mltSamples[proposedSample], maxDepth,
                      x0, x1, y0, y1, t0, t1);
        else
            SmallStep(rng, &mltSamples[proposedSample], maxDepth,
                      x0, x1, y0, y1, t0, t1);

        // Compute contribution of proposed sample and acceptance probability
        sampleLs[proposedSample] = L(scene, renderer, camera, arena, rng, maxDepth,
                                     ignoreDirect, mltSamples[proposedSample]);
        float currentI = I(sampleLs[currentSample], mltSamples[currentSample]);
        float proposedI = I(sampleLs[proposedSample], mltSamples[proposedSample]);
        float a = min(1.f, proposedI / currentI);
        float currentWeight = (1.f - a) /
                              (currentI / b + largeStepProbability) *
                              float(nPixels) / float(totalSamples);
        float proposedWeight = (a + (largeStep ? 1.f : 0.f)) /
                               (proposedI / b + largeStepProbability) *
                               float(nPixels) / float(totalSamples);

        // Splat current and proposed samples to _Film_
        if (currentWeight > 0.f && currentI > 0.f)
            camera->film->Splat(mltSamples[currentSample].cameraSample,
                sampleLs[currentSample] * currentWeight);
        if (proposedWeight > 0.f && proposedI > 0.f)
            camera->film->Splat(mltSamples[proposedSample].cameraSample,
               sampleLs[proposedSample] * proposedWeight);

        // Randomly accept proposed path mutation (or not)
        if (consecutiveRejects >= maxConsecutiveRejects ||
            rng.RandomFloat() < a) {
            PBRT_MLT_ACCEPTED_MUTATION(a, &mltSamples[currentSample], &mltSamples[proposedSample]);
            currentSample ^= 1;
            proposedSample ^= 1;
            consecutiveRejects = 0;
        }
        else {
            PBRT_MLT_REJECTED_MUTATION(a, &mltSamples[currentSample], &mltSamples[proposedSample]);
            ++consecutiveRejects;
        }
        arena.FreeAll();
    }
    // Update display for recently computed Metropolis samples
    float nf = AtomicAdd(nSamplesFinished, nSamples);
    float splatScale = float(totalSamples)/nf;
    camera->film->UpdateDisplay(x0, y0, x1, y1, splatScale);
    if ((taskNum % 32) == 0) {
        MutexLock lock(*filmMutex);
        camera->film->WriteImage(splatScale);
    }
    progress.Update();
    PBRT_MLT_FINISHED_MLT_TASK(this);
}

void MetropolisRenderer::Render(const Scene *scene) {
    int x0, x1, y0, y1;
    camera->film->GetPixelExtent(&x0, &x1, &y0, &y1);
    int nPixels = (x1-x0) * (y1-y0);
    float t0 = camera->shutterOpen;
    float t1 = camera->shutterClose;

    if (doDirectSeparately) {
        // Compute direct lighting before Metropolis light transport
        LDSampler sampler(x0, x1, y0, y1, directPixelSamples, t0, t1);
        Sample *sample = new Sample(&sampler, directLighting, NULL, scene);
        vector<Task *> directTasks;
        int nDirectTasks = max(32 * NumSystemCores(),
                         (camera->film->xResolution * camera->film->yResolution) / (16*16));
        nDirectTasks = RoundUpPow2(nDirectTasks);
        ProgressReporter directProgress(nDirectTasks, "Direct Lighting");
        for (int i = 0; i < nDirectTasks; ++i)
            directTasks.push_back(new SamplerRendererTask(scene, this, camera,
                &sampler, directProgress, sample, i, nDirectTasks));
        std::reverse(directTasks.begin(), directTasks.end());
        EnqueueTasks(directTasks);
        WaitForAllTasks();
        for (uint32_t i = 0; i < directTasks.size(); ++i)
            delete directTasks[i];
        delete sample;
        directProgress.Done();
    }
    // Take initial set of samples to compute $b$
    RNG rng(0);
    MemoryArena arena;
    vector<float> bootstrapSamples;
    float sumContrib = 0.f;
    bootstrapSamples.reserve(nBootstrap);
    MLTSample sample(maxDepth);
    for (int i = 0; i < nBootstrap; ++i) {
        // Compute contribution for random sample for MLT bootstrapping
        LargeStep(rng, &sample, maxDepth, x0, x1, y0, y1, t0, t1);
        float contrib = I(L(scene, this, camera, arena, rng, maxDepth,
                            doDirectSeparately, sample), sample);
        sumContrib += contrib;
        bootstrapSamples.push_back(contrib);
        arena.FreeAll();
    }
    float b = sumContrib / nBootstrap;

    // Select initial sample from bootstrap samples
    rng.Seed(0);
    float contribOffset = rng.RandomFloat() * sumContrib;
    sumContrib = 0.f;
    MLTSample initialSample(maxDepth);
    for (int i = 0; i < nBootstrap; ++i) {
        LargeStep(rng, &initialSample, maxDepth, x0, x1, y0, y1, t0, t1);
        sumContrib += bootstrapSamples[i];
        if (contribOffset < sumContrib)
            break;
    }

    // Launch tasks to generate Metropolis samples
    if (scene->lights.size() > 0) {
        int nTasks = int(nSamples / 50000);
        nTasks = max(nTasks, 32 * NumSystemCores());
        nTasks = min(nTasks, 32768);
        nSamples = (nSamples / nTasks) * nTasks;
        ProgressReporter progress(nTasks, "Metropolis");
        vector<Task *> tasks;
        Mutex *filmMutex = Mutex::Create();
        for (int i = 0; i < nTasks; ++i)
            tasks.push_back(new MLTTask(progress, i, int(nSamples/nTasks), nSamples,
                nPixels, x0, x1, y0, y1, t0, t1, b, largeStepProbability, initialSample,
                doDirectSeparately, maxConsecutiveRejects, maxDepth, scene, camera, this,
                &nSamplesFinished, filmMutex));
        EnqueueTasks(tasks);
        WaitForAllTasks();
        for (uint32_t i = 0; i < tasks.size(); ++i)
            delete tasks[i];
        progress.Done();
    }
    camera->film->WriteImage();
}

void LightShootingTask::Run() {
    // tady by mel byt kod z photon mappingu

    MemoryArena arena;
    uint32_t totalPaths = 0;
    RNG rng(seed);
    PermutedHalton halton(6, rng);
    while (true) {
        // Follow photon paths for a block of samples
        const uint32_t blockSize = 4096;
        for (uint32_t i = 0; i < blockSize; ++i) {
            float u[6];
            halton.Sample(++totalPaths, u);
            // Choose light to shoot photon from
            float lightPdf;
            int lightNum = lightDistribution->SampleDiscrete(u[0], &lightPdf);
            const Light *light = scene->lights[lightNum];

            // Generate _photonRay_ from light source and initialize _alpha_
            RayDifferential photonRay;
            float pdf;
            LightSample ls(u[1], u[2], u[3]);
            Normal Nl;
            Spectrum Le = light->Sample_L(scene, ls, u[4], u[5],time, &photonRay, &Nl, &pdf);
            if (pdf == 0.f || Le.IsBlack()) continue;
            Spectrum alpha = (AbsDot(Nl, photonRay.d) * Le) / (pdf * lightPdf);

            if (!alpha.IsBlack()) {
                // Follow photon path through scene and record intersections
                PBRT_PHOTON_MAP_STARTED_RAY_PATH(&photonRay, &alpha);
                bool specularPath = true;
                Intersection photonIsect;
                int nIntersections = 0;
                while (scene->Intersect(photonRay, &photonIsect)) {
                    ++nIntersections;
                    //MC tady by mel byt i kod pro volumetriku

                    // Handle photon/surface intersection
                    // alpha *= renderer->Transmittance(scene, photonRay, NULL, rng, arena);
                    BSDF *photonBSDF = photonIsect.GetBSDF(photonRay, arena);

                    Vector wo = -photonRay.d;
                    //MC tady se ukladaly photony takze tady bych mel ukladat samples do filmu kamery
                    //  // Deposit photon at surface
                    //Photon photon(photonIsect.dg.p, alpha, wo);
                    //tuhle metodu chci pouzit
                    //filmAddSample()

                    if (nIntersections >= maxDepth) break;

                    // Sample new photon ray direction
                    Vector wi;
                    float pdf;
                    BxDFType flags;
                    Spectrum fr = photonBSDF->Sample_f(wo, &wi, BSDFSample(rng),
                                                       &pdf, BSDF_ALL, &flags);

                    if (fr.IsBlack() || pdf == 0.f) break;
                    Spectrum anew = alpha * fr *
                                    AbsDot(wi, photonBSDF->dgShading.nn) / pdf;

                    // Possibly terminate photon path with Russian roulette
                    float continueProb = min(1.f, anew.y() / alpha.y());
                    if (rng.RandomFloat() > continueProb)
                        break;
                    alpha = anew / continueProb;
                    specularPath &= ((flags & BSDF_SPECULAR) != 0);

                    photonRay = RayDifferential(photonIsect.dg.p, wi, photonRay,
                                                photonIsect.rayEpsilon);

                }
                PBRT_PHOTON_MAP_FINISHED_RAY_PATH(&photonRay, &alpha);
            }

            arena.FreeAll();
        }

        //termination criteria ???
        if (totalPaths==maxPathCount) {
            break;
        }
    }
}