一、概述
一开始出现的是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 Definitions
void 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();
}四、个人笔记
如下是原始笔记的截图。记录的是最初的理解（有的错误修正了，有的没有修正）和过程中的“吐槽”。
277
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
