Q135：PBRT-V3，随机渐进光子映射（Stochastic Progressive Photon Mapping）（16.2章节）

一、概述

一开始出现的是Photon Mapping，后来有了Progressive Photon Mapping，后来才有了Stochastic Progressive Photon Mapping。
关于这三个的前世今生等相互关联，参考：
photon mapping学习笔记
再谈光子映射

前面的都不是重点。
咱关注的是PBRT-V3上介绍的SPPM（Stochastic Progressive Photon Mapping）的原理及其C++代码实现。

二、SPPM的原理

SPPM是PPM的改版（或者说“升级版”）。SPPM的优势是没有内存的限制。

SPPM算法有如下几个步骤：

1，由相机产生visible points（“visible point”，即“（对相机来说）可见的点”。基本是一个像素点对应一个visible point，所以很省内存。这个过程中会计算visible point的“直接光照”）。
2，将所有visible points保存到一个grid中。
3，由光源发出光子光线给grid中的visible points“送光”。
4，步骤1、2、3结束，意味着一次迭代将要结束，将要开始下一次迭代。在开始下一次迭代之前，根据当前迭代的情况调整下一次迭代光子“送光”是的搜索半径。然后，回到步骤1开始下一次迭代。

5、所有迭代完成后，计算每一个visible point最终的“间接光照”，在“直接光照”的基础上添加“间接光照”，输出图形。

2.1 产生visible points

前面已经提到“visible points”是指“对相机来说，可见的点”。
那么，场景中那些点对相机来说是可见的呢？
从相机产生一条光线进入场景，“visible point”指的是：
1，光线直接diffuse surface相交的交点；
2，光线发生specular (reflection/transmission) 之后，和diffuse surface相交的交点。

一般情况下（即，不考虑光线在specular surfaces之间弹来弹去和打出场景），一个像素点，对应一条光线，对应一个visible point。

保存visible point的数据：

其中，特别说明一下Ld、radius。

Ld，visible point处的直接光照；

radius，visible point的“搜索半径”/“收光半径”。对于每一个visible point，当前只保存了“直接光照”。后续，来自光源的光子光线给每个visible point送来的光是属于“间接光照”。光子光线“送光”的方式：光线和场景中某surface相交产生一个交点，该交点附近的所有visible points都有可能收到当前光子光线送来的光，只要visible point到该交点的距离小于或者等于该visible point的“搜索半径”/“收光半径”。

2.2 构建grid

grid的每个cell中保存着一个visible point链表。

前面提到，每一个visible point都有自己的“搜索半径”/“收光半径”。
这个就相当于一个“以visible point为球心，以“收光半径”为半径“的球。
visible point对应的这个球可能和grid中的多个cell发生重叠，所以需要将这个visible point添加到有重叠的多个cell的“visible point链表”中。

另外，可能存在多个visible point对应的球面和同一个cell发生重叠，所以，一个cell的“visible point链表”中会保存所有“对应球面和该cell发生重叠”的visible point。

2.3 光子光线“送光”

来自光源的光子光线和场景中的某surface相交，得到一个交点。
将该交点坐标转换到“前面根据visible point构建”的grid中。
该交点经过转换后会对应grid中的某个cell。
该cell中保存了一个visible point的链表。
这条光子光线将给这个cell对应的visible point链表中的所有visible points“送光”。光子光线给visible point送来的光即为visible point处的“间接光照”。

光子光线“送光”是这样的：送完一家，送另一家。

这个有点类似于path tracer中path的构建。

关于“送光”的量，这里也是用beta表示（这个beta和visible point中保存的那个beta的物理意义是不一样的。两个beta在程序中的作用范围不一样，所以不会冲突）。
初始beta（即光子光线对应path的长度为1时）：

（这里面涉及到光子光线的采样，参考16.1.2章节，此处不表）

随着光子光线对应path的长度的增加，beta值发生改变。
这里涉及到一个问题：光子光线不能是无止尽地在场景中延伸，那么怎么终止光子光线呢？
一个光子撞击物体后，是继续传播呢，还是被物体吸收呢？这个用Russian roulette来决定。
β的计算如下：

2.4 调整visible point的“搜索半径”

对于每一迭代（对应若干条光子光线），相关数据的计算如下放截图：

当最后一次迭代结束时，计算光源给某个visible point送光所产生的效果（visible point的间接光照）：

正因为“搜索半径”的存在，所以SPPM是一个biased的算法。
当然，实际使用中，随着迭代次数的增加，“搜索半径”越来越小，从而biased的程度越来越小，得到的图形也越来越精确。

2.5 关于单次迭代的光子数和迭代次数

直接看张图吧：

如a图，单次迭代的光子数为10000个，迭代1000次。
（关于单次迭代的光子数，PBRT-V3官方代码中做法：若用户有设置单次迭代的光子数，就用这个数值；若用户没有设置单次迭代的光子数，就把像素点的个数作为单次迭代的光子数的默认值。
关于迭代次数，由于不涉及内存问题，可由用户随意设置）

三、SPPM的C++代码实现

PBRT-V3中对应的积分器是SPPMIntegrator。
SPPMIntegrator不是SamplerIntegrator，所以要实现自己的render()成员方法。

先截取需要特别说明的代码段，在贴出完整的render()成员方法。

产生visible point时：

给visible point送光时：

调整visible point的“搜索半径”时：

给visible point添加最终的间接光照时：

完整render()代码如下：

// SPPM Method Definitionsvoid SPPMIntegrator::Render(const Scene &scene) {    ProfilePhase p(Prof::IntegratorRender);    // Initialize _pixelBounds_ and _pixels_ array for SPPM    Bounds2i pixelBounds = camera->film->croppedPixelBounds;    int nPixels = pixelBounds.Area();    std::unique_ptr<SPPMPixel[]> pixels(new SPPMPixel[nPixels]);    for (int i = 0; i < nPixels; ++i) pixels[i].radius = initialSearchRadius;    const Float invSqrtSPP = 1.f / std::sqrt(nIterations);    pixelMemoryBytes = nPixels * sizeof(SPPMPixel);    // Compute _lightDistr_ for sampling lights proportional to power    std::unique_ptr<Distribution1D> lightDistr =        ComputeLightPowerDistribution(scene);    // Perform _nIterations_ of SPPM integration    HaltonSampler sampler(nIterations, pixelBounds);    // Compute number of tiles to use for SPPM camera pass    Vector2i pixelExtent = pixelBounds.Diagonal();    const int tileSize = 16;    Point2i nTiles((pixelExtent.x + tileSize - 1) / tileSize,                   (pixelExtent.y + tileSize - 1) / tileSize);    ProgressReporter progress(2 * nIterations, "Rendering");    for (int iter = 0; iter < nIterations; ++iter) {        // Generate SPPM visible points        std::vector<MemoryArena> perThreadArenas(MaxThreadIndex());        {            ProfilePhase _(Prof::SPPMCameraPass);            ParallelFor2D([&](Point2i tile) {                MemoryArena &arena = perThreadArenas[ThreadIndex];                // Follow camera paths for _tile_ in image for SPPM                int tileIndex = tile.y * nTiles.x + tile.x;                std::unique_ptr<Sampler> tileSampler = sampler.Clone(tileIndex);                // Compute _tileBounds_ for SPPM tile                int x0 = pixelBounds.pMin.x + tile.x * tileSize;                int x1 = std::min(x0 + tileSize, pixelBounds.pMax.x);                int y0 = pixelBounds.pMin.y + tile.y * tileSize;                int y1 = std::min(y0 + tileSize, pixelBounds.pMax.y);                Bounds2i tileBounds(Point2i(x0, y0), Point2i(x1, y1));                for (Point2i pPixel : tileBounds) {                    // Prepare _tileSampler_ for _pPixel_                    tileSampler->StartPixel(pPixel);                    tileSampler->SetSampleNumber(iter);                    // Generate camera ray for pixel for SPPM                    CameraSample cameraSample =                        tileSampler->GetCameraSample(pPixel);                    RayDifferential ray;                    Spectrum beta =                        camera->GenerateRayDifferential(cameraSample, &ray);                    ray.ScaleDifferentials(invSqrtSPP);                    // Follow camera ray path until a visible point is created                    // Get _SPPMPixel_ for _pPixel_                    Point2i pPixelO = Point2i(pPixel - pixelBounds.pMin);                    int pixelOffset =                        pPixelO.x +                        pPixelO.y * (pixelBounds.pMax.x - pixelBounds.pMin.x);                    SPPMPixel &pixel = pixels[pixelOffset];                    bool specularBounce = false;                    for (int depth = 0; depth < maxDepth; ++depth) {                        SurfaceInteraction isect;                        ++totalPhotonSurfaceInteractions;                        if (!scene.Intersect(ray, &isect)) {                            // Accumulate light contributions for ray with no                            // intersection                            for (const auto &light : scene.lights)                                pixel.Ld += beta * light->Le(ray);                            break;                        }                        // Process SPPM camera ray intersection                        // Compute BSDF at SPPM camera ray intersection                        isect.ComputeScatteringFunctions(ray, arena, true);                        if (!isect.bsdf) {                            ray = isect.SpawnRay(ray.d);                            --depth;                            continue;                        }                        const BSDF &bsdf = *isect.bsdf;                        // Accumulate direct illumination at SPPM camera ray                        // intersection                        Vector3f wo = -ray.d;                        if (depth == 0 || specularBounce)                            pixel.Ld += beta * isect.Le(wo);                        pixel.Ld +=                            beta * UniformSampleOneLight(isect, scene, arena,                                                         *tileSampler);                        // Possibly create visible point and end camera path                        bool isDiffuse = bsdf.NumComponents(BxDFType(                                             BSDF_DIFFUSE | BSDF_REFLECTION |                                             BSDF_TRANSMISSION)) > 0;                        bool isGlossy = bsdf.NumComponents(BxDFType(                                            BSDF_GLOSSY | BSDF_REFLECTION |                                            BSDF_TRANSMISSION)) > 0;                        if (isDiffuse || (isGlossy && depth == maxDepth - 1)) {                            pixel.vp = {isect.p, wo, &bsdf, beta};                            break;                        }                        // Spawn ray from SPPM camera path vertex                        if (depth < maxDepth - 1) {                            Float pdf;                            Vector3f wi;                            BxDFType type;                            Spectrum f =                                bsdf.Sample_f(wo, &wi, tileSampler->Get2D(),                                              &pdf, BSDF_ALL, &type);                            if (pdf == 0. || f.IsBlack()) break;                            specularBounce = (type & BSDF_SPECULAR) != 0;                            beta *= f * AbsDot(wi, isect.shading.n) / pdf;                            if (beta.y() < 0.25) {                                Float continueProb =                                    std::min((Float)1, beta.y());                                if (tileSampler->Get1D() > continueProb) break;                                beta /= continueProb;                            }                            ray = (RayDifferential)isect.SpawnRay(wi);                        }                    }                }            }, nTiles);        }        progress.Update();        // Create grid of all SPPM visible points        int gridRes[3];        Bounds3f gridBounds;        // Allocate grid for SPPM visible points        const int hashSize = nPixels;        std::vector<std::atomic<SPPMPixelListNode *>> grid(hashSize);        {            ProfilePhase _(Prof::SPPMGridConstruction);            // Compute grid bounds for SPPM visible points            Float maxRadius = 0.;            for (int i = 0; i < nPixels; ++i) {                const SPPMPixel &pixel = pixels[i];                if (pixel.vp.beta.IsBlack()) continue;                Bounds3f vpBound = Expand(Bounds3f(pixel.vp.p), pixel.radius);                gridBounds = Union(gridBounds, vpBound);                maxRadius = std::max(maxRadius, pixel.radius);            }            // Compute resolution of SPPM grid in each dimension            Vector3f diag = gridBounds.Diagonal();            Float maxDiag = MaxComponent(diag);            int baseGridRes = (int)(maxDiag / maxRadius);            CHECK_GT(baseGridRes, 0);            for (int i = 0; i < 3; ++i)                gridRes[i] = std::max((int)(baseGridRes * diag[i] / maxDiag), 1);            // Add visible points to SPPM grid            ParallelFor([&](int pixelIndex) {                MemoryArena &arena = perThreadArenas[ThreadIndex];                SPPMPixel &pixel = pixels[pixelIndex];                if (!pixel.vp.beta.IsBlack()) {                    // Add pixel's visible point to applicable grid cells                    Float radius = pixel.radius;                    Point3i pMin, pMax;                    ToGrid(pixel.vp.p - Vector3f(radius, radius, radius),                           gridBounds, gridRes, &pMin);                    ToGrid(pixel.vp.p + Vector3f(radius, radius, radius),                           gridBounds, gridRes, &pMax);                    for (int z = pMin.z; z <= pMax.z; ++z)                        for (int y = pMin.y; y <= pMax.y; ++y)                            for (int x = pMin.x; x <= pMax.x; ++x) {                                // Add visible point to grid cell $(x, y, z)$                                int h = hash(Point3i(x, y, z), hashSize);                                SPPMPixelListNode *node =                                    arena.Alloc<SPPMPixelListNode>();                                node->pixel = &pixel;                                // Atomically add _node_ to the start of                                // _grid[h]_'s linked list                                node->next = grid[h];                                while (grid[h].compare_exchange_weak(                                           node->next, node) == false)                                    ;                            }                    ReportValue(gridCellsPerVisiblePoint,                                (1 + pMax.x - pMin.x) * (1 + pMax.y - pMin.y) *                                    (1 + pMax.z - pMin.z));                }            }, nPixels, 4096);        }        // Trace photons and accumulate contributions        {            ProfilePhase _(Prof::SPPMPhotonPass);            std::vector<MemoryArena> photonShootArenas(MaxThreadIndex());            ParallelFor([&](int photonIndex) {                MemoryArena &arena = photonShootArenas[ThreadIndex];                // Follow photon path for _photonIndex_                uint64_t haltonIndex =                    (uint64_t)iter * (uint64_t)photonsPerIteration +                    photonIndex;                int haltonDim = 0;                // Choose light to shoot photon from                Float lightPdf;                Float lightSample = RadicalInverse(haltonDim++, haltonIndex);                int lightNum =                    lightDistr->SampleDiscrete(lightSample, &lightPdf);                const std::shared_ptr<Light> &light = scene.lights[lightNum];                // Compute sample values for photon ray leaving light source                Point2f uLight0(RadicalInverse(haltonDim, haltonIndex),                                RadicalInverse(haltonDim + 1, haltonIndex));                Point2f uLight1(RadicalInverse(haltonDim + 2, haltonIndex),                                RadicalInverse(haltonDim + 3, haltonIndex));                Float uLightTime =                    Lerp(RadicalInverse(haltonDim + 4, haltonIndex),                         camera->shutterOpen, camera->shutterClose);                haltonDim += 5;                // Generate _photonRay_ from light source and initialize _beta_                RayDifferential photonRay;                Normal3f nLight;                Float pdfPos, pdfDir;                Spectrum Le =                    light->Sample_Le(uLight0, uLight1, uLightTime, &photonRay,                                     &nLight, &pdfPos, &pdfDir);                if (pdfPos == 0 || pdfDir == 0 || Le.IsBlack()) return;                Spectrum beta = (AbsDot(nLight, photonRay.d) * Le) /                                (lightPdf * pdfPos * pdfDir);                if (beta.IsBlack()) return;                // Follow photon path through scene and record intersections                SurfaceInteraction isect;                for (int depth = 0; depth < maxDepth; ++depth) {                    if (!scene.Intersect(photonRay, &isect)) break;                    ++totalPhotonSurfaceInteractions;                    if (depth > 0) {                        // Add photon contribution to nearby visible points                        Point3i photonGridIndex;                        if (ToGrid(isect.p, gridBounds, gridRes,                                   &photonGridIndex)) {                            int h = hash(photonGridIndex, hashSize);                            // Add photon contribution to visible points in                            // _grid[h]_                            for (SPPMPixelListNode *node =                                     grid[h].load(std::memory_order_relaxed);                                 node != nullptr; node = node->next) {                                ++visiblePointsChecked;                                SPPMPixel &pixel = *node->pixel;                                Float radius = pixel.radius;                                if (DistanceSquared(pixel.vp.p, isect.p) >                                    radius * radius)                                    continue;                                // Update _pixel_ $\Phi$ and $M$ for nearby                                // photon                                Vector3f wi = -photonRay.d;                                Spectrum Phi =                                    beta * pixel.vp.bsdf->f(pixel.vp.wo, wi);                                for (int i = 0; i < Spectrum::nSamples; ++i)                                    pixel.Phi[i].Add(Phi[i]);                                ++pixel.M;                            }                        }                    }                    // Sample new photon ray direction                    // Compute BSDF at photon intersection point                    isect.ComputeScatteringFunctions(photonRay, arena, true,                                                     TransportMode::Importance);                    if (!isect.bsdf) {                        --depth;                        photonRay = isect.SpawnRay(photonRay.d);                        continue;                    }                    const BSDF &photonBSDF = *isect.bsdf;                    // Sample BSDF _fr_ and direction _wi_ for reflected photon                    Vector3f wi, wo = -photonRay.d;                    Float pdf;                    BxDFType flags;                    // Generate _bsdfSample_ for outgoing photon sample                    Point2f bsdfSample(                        RadicalInverse(haltonDim, haltonIndex),                        RadicalInverse(haltonDim + 1, haltonIndex));                    haltonDim += 2;                    Spectrum fr = photonBSDF.Sample_f(wo, &wi, bsdfSample, &pdf,                                                      BSDF_ALL, &flags);                    if (fr.IsBlack() || pdf == 0.f) break;                    Spectrum bnew =                        beta * fr * AbsDot(wi, isect.shading.n) / pdf;                    // Possibly terminate photon path with Russian roulette                    Float q = std::max((Float)0, 1 - bnew.y() / beta.y());                    if (RadicalInverse(haltonDim++, haltonIndex) < q) break;                    beta = bnew / (1 - q);                    photonRay = (RayDifferential)isect.SpawnRay(wi);                }                arena.Reset();            }, photonsPerIteration, 8192);            progress.Update();            photonPaths += photonsPerIteration;        }        // Update pixel values from this pass's photons        {            ProfilePhase _(Prof::SPPMStatsUpdate);            ParallelFor([&](int i) {                SPPMPixel &p = pixels[i];                if (p.M > 0) {                    // Update pixel photon count, search radius, and $\tau$ from                    // photons                    Float gamma = (Float)2 / (Float)3;                    Float Nnew = p.N + gamma * p.M;                    Float Rnew = p.radius * std::sqrt(Nnew / (p.N + p.M));                    Spectrum Phi;                    for (int j = 0; j < Spectrum::nSamples; ++j)                        Phi[j] = p.Phi[j];                    p.tau = (p.tau + p.vp.beta * Phi) * (Rnew * Rnew) /                            (p.radius * p.radius);                    p.N = Nnew;                    p.radius = Rnew;                    p.M = 0;                    for (int j = 0; j < Spectrum::nSamples; ++j)                        p.Phi[j] = (Float)0;                }                // Reset _VisiblePoint_ in pixel                p.vp.beta = 0.;                p.vp.bsdf = nullptr;            }, nPixels, 4096);        }        // Periodically store SPPM image in film and write image        if (iter + 1 == nIterations || ((iter + 1) % writeFrequency) == 0) {            int x0 = pixelBounds.pMin.x;            int x1 = pixelBounds.pMax.x;            uint64_t Np = (uint64_t)(iter + 1) * (uint64_t)photonsPerIteration;            std::unique_ptr<Spectrum[]> image(new Spectrum[pixelBounds.Area()]);            int offset = 0;            for (int y = pixelBounds.pMin.y; y < pixelBounds.pMax.y; ++y) {                for (int x = x0; x < x1; ++x) {                    // Compute radiance _L_ for SPPM pixel _pixel_                    const SPPMPixel &pixel =                        pixels[(y - pixelBounds.pMin.y) * (x1 - x0) + (x - x0)];                    Spectrum L = pixel.Ld / (iter + 1);                    L += pixel.tau / (Np * Pi * pixel.radius * pixel.radius);                    image[offset++] = L;                }            }            camera->film->SetImage(image.get());            camera->film->WriteImage();            // Write SPPM radius image, if requested            if (getenv("SPPM_RADIUS")) {                std::unique_ptr<Float[]> rimg(                    new Float[3 * pixelBounds.Area()]);                Float minrad = 1e30f, maxrad = 0;                for (int y = pixelBounds.pMin.y; y < pixelBounds.pMax.y; ++y) {                    for (int x = x0; x < x1; ++x) {                        const SPPMPixel &p =                            pixels[(y - pixelBounds.pMin.y) * (x1 - x0) +                                   (x - x0)];                        minrad = std::min(minrad, p.radius);                        maxrad = std::max(maxrad, p.radius);                    }                }                fprintf(stderr,                        "iterations: %d (%.2f s) radius range: %f - %f\n",                        iter + 1, progress.ElapsedMS() / 1000., minrad, maxrad);                int offset = 0;                for (int y = pixelBounds.pMin.y; y < pixelBounds.pMax.y; ++y) {                    for (int x = x0; x < x1; ++x) {                        const SPPMPixel &p =                            pixels[(y - pixelBounds.pMin.y) * (x1 - x0) +                                   (x - x0)];                        Float v = 1.f - (p.radius - minrad) / (maxrad - minrad);                        rimg[offset++] = v;                        rimg[offset++] = v;                        rimg[offset++] = v;                    }                }                Point2i res(pixelBounds.pMax.x - pixelBounds.pMin.x,                            pixelBounds.pMax.y - pixelBounds.pMin.y);                WriteImage("sppm_radius.png", rimg.get(), pixelBounds, res);            }        }    }    progress.Done();}

四、个人笔记

如下是原始笔记的截图。记录的是最初的理解（有的错误修正了，有的没有修正）和过程中的“吐槽”。