## 16.2 Stochastic Progressive Photon Mapping

Photon mapping is one of a family of particle-tracing algorithms, which are based on the idea of constructing paths starting from the lights and connecting vertices in these paths to the camera to deposit energy on the film. In this section, we will start by introducing a theory of particle-tracing algorithms and will discuss the conditions that must be fulfilled by a particle-tracing algorithm so that arbitrary measurements can be computed correctly using the particles created by the algorithm. We will then describe an implementation of a photon mapping integrator that uses particles to estimate illumination by interpolating lighting contributions from particles close to but not quite at the point being shaded.

### 16.2.1 Theoretical Basis for Particle Tracing

Particle-tracing algorithms in computer graphics are often explained in terms of packets of energy being shot from the light sources in the scene that deposit energy at surfaces they intersect before scattering in new directions. This is an intuitive way of thinking about particle tracing, but the intuition that it provides doesn’t make it easy to answer basic questions about how propagation and scattering affect the particles. For example, does their contribution fall off with squared distance like flux density? Or, which terms, if any, affect particles after they scatter from a surface?

In order to give a solid theoretical basis for particle tracing, we will describe it using a framework introduced by Veach (1997, Appendix 4.A), which instead interprets the stored particle histories as samples from the scene’s equilibrium radiance distribution. Under certain conditions on the distribution and weights of the particles, the particles can be used to compute estimates of nearly any measurement based on the light distribution in the scene. In this framework, it is quite easy to answer questions about the details of particle propagation like the ones earlier. After developing this theory here, the remainder of this section will demonstrate its application to photon mapping.

A particle-tracing algorithm generates a set of samples of illumination at points , on surfaces in the scene

where each sample records incident illumination from direction and has some throughput weight associated with it (Figure 16.5). As the notation already indicates, this weight will contain ratios of terms of the throughput function and the associated sampling PDFs much like the variable of the path tracer (Section 14.5.4). We would like to determine the conditions on the weights and distribution of particle positions so that we can use them to correctly compute estimates of arbitrary measurements.

Given an importance function that describes the measurement to be taken, the natural condition we would like to be fulfilled is that the particles should be distributed and weighted such that using them to compute an estimate has the same expected value as the measurement equation for the same importance function:

For example, we might want to use the particles to compute the total flux incident on a wall. Using the definition of flux,

the following importance function selects the particles that lie on the wall and arrived from the hemisphere around the normal:

If the conditions on the distribution of particle weights and positions are true for arbitrary importance functions such that Equation (16.7) holds, then the flux estimate can be computed directly as a sum of the particle weights for the particles on the wall. If we want to estimate flux over a different wall, a subset of the original wall, and so on, we only need to recompute the weighted sum with an updated importance function. The particles and weights can be reused, and we have an unbiased estimate for all of these measurements. (The estimates will be correlated, however, which is potentially a source of image artifacts.)

To see how to generate and weight particles that fulfill these conditions, consider the task of evaluating the measurement equation integral

where the vertex densities are expressed as a probability per unit area and where the importance function that describes the measurement is a black box and thus cannot be used to drive the sampling of the integral at all. We can still compute an estimate of the integral with Monte Carlo integration but must sample a set of points and from all of the surfaces in the scene, using some sampling distribution that doesn’t depend on (e.g., by uniformly sampling points by surface area).

By expanding the LTE in the integrand and applying the standard Monte Carlo estimator for samples, we can find the estimator for this measurement,

We can further expand out the term into the sum over paths and use the fact that and the fact that for a particular sample, the expected value

can be written as a finite sum of terms in just the same way that we generated a finite set of weighted path vertices for path tracing. If the sum is truncated with Russian roulette such that the probability of terminating the sum after terms is , then the th term of the th sample has contribution

Note that the path integral framework provides us with the freedom to generate a set of particles in all sorts of different ways—i.e., with different underlying vertex probability densities . Although the natural approach is to start from points on lights and incrementally sample paths using the BSDFs at the path vertices, similar to how the path-tracing integrator generates paths (starting here from the light, rather than from the camera), we could generate them with any number of different sampling strategies, as long as there is nonzero probability of generating a particle at any point where the numerator is nonzero and the particle weights are computed using the above definition.

If we only had a single measurement to make, it would be better if we used information about and could compute the estimate more intelligently, since the general particle-tracing approach described here may generate many useless samples if only covers a small subset of the points on scene objects. If we will be computing many measurements, however, the key advantage that particle tracing brings is that we can generate the samples and weights once and can then reuse them over a large number of measurements, potentially computing results much more efficiently than if the measurements were all computed from scratch.

### 16.2.2 Photon Mapping

The photon mapping algorithm is based on tracing particles into the scene and blurring their contribution to approximate the incident illumination at shading points. For consistency with other descriptions of the algorithm, we will refer to particles generated for photon mapping as photons.

In order to compute reflected radiance at a point, we need to estimate the exitant radiance equation at a point in a direction , which can equivalently (and cumbersomely) be written as a measurement over all points on surfaces in the scene where a Dirac delta distribution selects only particles precisely at :

Unfortunately, because there is a delta distribution in , all of the particle histories that were generated during the particle-tracing step have zero probability of having nonzero contribution if Equation (16.7) is used to compute the estimate of the measurement value (just as we will never be able to choose a direction from a diffuse surface that intersects a point light source unless the direction is sampled based on the light source’s position).

Here is the point at which bias is introduced into the photon mapping algorithm. Under the assumption that the information about illumination at nearby points gives a reasonable approximation to illumination at the shading point, photon mapping interpolates illumination from nearby photons around the point being shaded; the delta function of position in Equation (16.9) is effectively converted to a filter function. Given Equation (16.7), the more photons there are around the point and the higher their weights, the more radiance is estimated to be incident at the point.

One factor that contributes to photon mapping’s efficiency is this reuse of photons: having gone through the trouble to compute a light-carrying path, allowing it to potentially contribute illumination at multiple points amortizes the cost of generating it. While photon mapping derives some benefit from this efficiency improvement, there’s a subtler but much more important benefit from using nearby photons at a point: some light-carrying paths are impossible to sample with unbiased algorithms based on incremental path construction (including path tracing and bidirectional path tracing), but are handled well with photon mapping. These paths can arise in very common situations.

For example, consider the task of rendering an image of a photograph with a plate of glass in front of it. Assume a pinhole camera model and a point light illuminating the scene, and also assume for simplicity that the glass is only transmissive (Figure 16.6). If we start a path from the camera that passes through the glass, then the point that it intersects the photograph is completely determined by the effect of the refraction. At this point, none of the direct lighting strategies we have available has any chance of sampling an incident direction that will reach the light: because any sampled incident direction leaving the diffuse surface will be refracted on the way out through the glass, there’s no chance that the refracted ray will hit the point light. Even with an area light source, some refracted rays may be lucky enough to hit the light, but in general, variance will be high since most will miss it.

With photon mapping, we can trace photons leaving the light, let them refract through the glass, and deposit illumination on the diffuse surface. With a sufficient number of photons, the surface will be densely covered, and the photons around a point give a good estimate of the incident illumination.

A statistical technique called *density estimation*
provides the
mathematical tools to perform this interpolation. Density estimation
constructs a PDF given a set of sample points under the assumption that the
samples are distributed according to the overall distribution of some
function of interest.
Histogramming is a straightforward example of the
idea. In 1D, the line is divided into intervals with some width, and one
can count how many samples land in each interval and normalize so that the
areas of the intervals sum to one.

*Kernel methods* are a more sophisticated density estimation
technique. They generally give better results and smoother
PDFs that don’t suffer from the discontinuities that histograms do. Given a
kernel function that integrates to 1.

the kernel estimator for samples at locations is

where is the window width (also known as the *smoothing
parameter* or *kernel bandwidth*).
Kernel methods can be thought of as placing a series of bumps at
observation points, where the sum of the bumps forms a PDF since they
individually integrate to 1 and the sum is normalized.
Figure 16.7 shows an example of density estimation in 1D,
where a smooth PDF is computed from a set of sample points.

The key question with kernel methods is how the window width is chosen.
If it is too wide, the PDF will blur out relevant detail in parts of the
domain with many samples; if it is too narrow, the PDF will be too
bumpy in the tails of the distribution where there aren’t many samples.
Nearest-neighbor techniques solve this problem by choosing adaptively
based on local density of samples. Where there are many samples, the width
is small; where there are few samples, the width is large. For
example, one approach is to pick a number and find the distance to the
th nearest sample from the point and use that distance, , for
the window width. This is the *generalized nth nearest-neighbor
estimate:*

In dimensions, this generalizes to

Substituting into the measurement equation, it can be shown that the appropriate estimator for the measurement we’d like to compute, the exitant radiance at the point in direction , is given by

where we’ve switched to using to denote the total number of emitted photons, the sum is over all of the photons, and scale factors for the photons are computed based on the density estimation Equation, (16.10). Because we know that the kernel function is zero for points farther away than the th nearest neighbor distance , implementations of this sum only need to sum over the closest neighbors.

The error introduced by this interpolation can be difficult to quantify. Tracing more photons generally increases photon density and will almost always improve the results. When the photons are more closely spaced, it isn’t necessary to use photons from as far away in the nearest-neighbor estimate. In general, the error at any given point will depend on how quickly the illumination is changing in its vicinity. One can always construct pathological cases where this error is unacceptable, but in practice it usually isn’t too bad. Because the interpolation step tends to blur out illumination, high-frequency changes in lighting are sometimes poorly reconstructed with photon mapping. If traditional methods are used for direct lighting, then this is generally less of a problem since indirect illumination tends to be low frequency.

The original formulation of photon mapping was based on a two-pass algorithm where photons are first traced from the light sources. Photon interactions are recorded on surfaces of the scene and then organized in a spatial data structure (generally a kd-tree) for use at rendering time. A second pass follows paths starting from the camera; at each path vertex, nearby photons are used to estimate indirect illumination.

While this approach is effective, it has the limitation that the number of photons that can be used is limited by available memory since all of the photons must be stored. No improvement is possible once memory is full, even if one wants a higher quality result from using more photons. (In contrast, with path tracing, for example, one can always add more samples per pixel without any incremental storage cost—only the computational expense increases.)

#### Progressive Photon Mapping

The *progressive photon mapping* algorithm addressed this issue by
restructuring the algorithm: first, a camera pass traces paths starting
from the camera. Each pixel stores a representation of all of the non-specular path
vertices found when generating paths that started in its extent. (For example,
if a camera ray
hits a diffuse surface, we might record the geometric information about the
intersection point and the diffuse reflectance. If it hit a perfectly
specular surface and then a diffuse surface, we would record the diffuse
reflectance scaled by the specular BSDF value, and so forth.) We will dub
these stored path vertices *visible points* in the
following.
A second pass traces photons from the
light sources; at each photon–surface intersection, the photon contributes
to the reflected radiance estimate for nearby visible points.

To understand how progressive photon mapping works, we consider a decomposition of the LTE into separate integrals over direct and indirect incident radiance at each vertex. (Recall the discussion of partitioning the LTE in Section 14.4.6.)

The emitted term is straightforward, and direct lighting can be handled using the usual approaches from Section 14.3. The indirect term with the integral over is handled in one of two ways. First, a ray may be sampled from the BSDF’s sampling distribution and traced to find the next vertex in the path, where the same issue of how to estimate outgoing radiance is faced—just like a path tracer. Alternatively, the current point can be saved to receive illumination from photons. The final contribution of such a point to radiance at the film plane is found by summing the product of the photons’ weights and the path throughput weight of the sequence of vertices before it in the path.

For a perfectly specular BSDF, the only reasonable choice is to trace another ray: a photon will never arrive at exactly the right direction to match the delta distribution in the BSDF. For highly glossy surfaces, it’s also advisable to trace a ray, since it will take many photons for enough to hit a narrow specular lobe to compute an accurate estimate.

For diffuse surfaces, it’s generally worth using photons, though it can be
worthwhile to trace one last bounce. This approach is called *final
gathering;* if photons are used only after a diffuse bounce, then any
errors from insufficient photon density are generally less visible, though
more camera paths may need to be traced to eliminate noise. (See
Exercise 16.8 for further discussion of final gathering.)

With this approach, no photon storage is needed, and an arbitrary number of photons can be traced; the memory limit is instead tied to the storage required for the visible points and their reflection information. For high-resolution images or images that require many samples per pixel to resolve motion blur or depth of field, memory can still be a limitation.

#### Stochastic Progressive Photon Mapping

Stochastic progressive photon mapping is a modification of progressive photon mapping that doesn’t suffer from either of these memory limits. Like progressive photon mapping, it generates a set of visible points from the camera but at a low sampling rate (e.g., following just one camera path per pixel). Next, a number of photons are shot from lights, accumulating contributions at nearby visible points. This process then repeats: the visible points are discarded, a new set is generated at different positions, another round of photons is traced, and so forth.

SPPM starts with the photon estimation equation, (16.11), and makes two adjustments. First, it uses a constant kernel function; in conjunction with the fact that the estimation is over 2D (the local tangent plane around the visible point), we have

where, as before, is the number of photons emitted from the light sources, and is the surface area of the disk-shaped kernel function.

The second adjustment, based on an approach first implemented in progressive photon mapping, is to progressively reduce the photon search radius as more photons contribute to the visible point. The general idea is that as more photons are found within the search radius, we have more evidence that a sufficient density of photons is arriving to estimate the incident radiance distribution well. By reducing the radius, we ensure that future photons that are used will be closer to the point and thus contribute to a more accurate estimate of the incident radiance distribution.

Reducing the radius requires an adjustment to how the reflected radiance estimate is computed, as now the photons in the sum in Equation (16.11) come from different radii. The following three update rules describe how to update the radius and dependent variables:

where is the number of photons that have contributed to the point after the th iteration, is the number of photons that contributed during the current iteration, is the search radius to use for the th iteration, maintains the sum of products of photons with BSDF values, and is computed during the th iteration as

The parameter, which is typically around , determines how quickly the contributions from photons from earlier iterations, with wider search radii, are faded out. (Hachisuka and Jensen’s original paper on SPPM (2009) uses the notation for this quantity; we opt for here, having used for other quantities already.)

Note that the radius is a per-pixel property, not a per-visible-point property. Remarkably, a consistent estimate for the reflected radiance is computed even with this single radius shared over all of the series of visible points in the pixel. We won’t show this derivation here, but with the rules in Equation (16.13), as the number of photons traced , , and the reflected radiance estimates are consistent and converge to the correct values.

Figure 16.8 shows a rendered scene with path tracing and SPPM. SPPM is much more effective at handling light that passes through the glass light fixtures than the path tracing algorithm is.

*Scene courtesy “Wig42” from blendswap.com*)

### 16.2.3 SPPMIntegrator

The `SPPMIntegrator`, implemented in the files
`integrators/sppm.h` and
`integrators/sppm.cpp`, implements
the SPPM light transport algorithm.

The `SPPMIntegrator` constructor isn’t particularly interesting; it
just sets various member variables with values passed in. We therefore
won’t include it here but will discuss the various member variables that
configure the `SPPMIntegrator`’s operation as they appear in the
following.

The `SPPMIntegrator` is not a `SamplerIntegrator`, so it implements
its own `Render()` method. After some initial setup has been
performed, it runs a number of iterations of the SPPM algorithm, finding a
set of visible points and then accumulating illumination from photons at
them. Each iteration creates a new path starting from the camera in each
pixel, which helps with antialiasing geometric edges and sampling motion
blur and depth of field well.

`pixelBounds`and

`pixels`array for SPPM>>

`lightDistr`for sampling lights proportional to power>>

`nIterations`of SPPM integration>>

`tile`in image for SPPM>>

`tileBounds`for SPPM tile>>

`tileSampler`for

`pPixel`>>

`SPPMPixel`for

`pPixel`>>

`photonIndex`>>

`photonRay`from light source and initialize

`beta`>>

`grid[h]`>> }

`fr`and direction

`wi`for reflected photon>>

`bsdfSample`for outgoing photon sample>>

`VisiblePoint`in pixel>> }

`L`for SPPM pixel

`pixel`>> image[offset++] = L; } } camera->film->SetImage(image.get()); camera->film->WriteImage(); }

The `pixels` array stores a `SPPMPixel` (to be defined shortly) for each pixel in the
final image.

`pixelBounds`and

`pixels`array for SPPM>>=

A user-supplied radius, `initialSearchRadius`, is used for , the
initial search radius for photons. If the supplied radius is too large,
too many photons will contribute to visible points during early iterations
(before the radius is automatically decreased), which may be inefficient.
If it’s too small, not enough photons will be found to estimate incident
radiance well. A radius corresponding to a few pixels in the final image
is generally a good starting point.

The `SPPMPixel` structure serves three purposes. First, it stores the
current estimated average radiance visible over the extent of a pixel
(including the time the shutter is open and accounting for depth of field,
if present). Second, it stores parameters related to the photon density
estimation for the pixel (e.g., various quantities from
Equation (16.13)). Finally, it stores the geometric and
reflection information for a visible point in the pixel after the camera
pass.

Working with weighted photons allows for two very different overall sampling approaches: on the one hand, we could try to distribute photons uniformly with a weight that approximates the irradiance on surfaces, while also accounting for discontinuities and other important geometric and lighting features. However, this type of photon distribution is challenging to realize for general input; hence we will instead strive to generate photons that have the same (or similar) weights, so that it is their varying density throughout the scene that represents the variation in illumination.

Furthermore, if the photon weights on a surface have substantial variation that is not related to the irradiance, there can be unpleasant image artifacts: if one photon takes on a much larger weight than the others, a bright circular artifact will reveal the region of the scene where that photon contributes to the interpolated radiance.

Therefore, we’d like to shoot more photons from the brighter lights so that the initial weights of photons leaving all lights will be of similar magnitudes, and thus, the light to start each photon path from is chosen according to a PDF defined by the lights’ respective powers. Thus, it is a greater number of photons from the brighter lights that accounts for their greater contributions to illumination in the scene rather than the same number of photons from all, with larger weights for the more powerful lights.

`lightDistr`for sampling lights proportional to power>>=

Each iteration of the SPPM algorithm traces a new path from the camera at each pixel and then collects incident photons at each path’s endpoint. Figure 16.9 shows the effect of increasing the number of iterations with SPPM. Here, we see the caustic from light passing through the glass becoming increasingly sharp. In general, more iterations improve the sampling of visible edges, motion blur, and depth of field, but as more photons are accumulated in each pixel, the indirect lighting estimate becomes more accurate. Note how the artifacts from under-sampling here are low-frequency blotches, a quite different visual artifact from the high-frequency noise from under-sampling seen in path tracing.

`SPPMIntegrator`. As the number of iterations increases (and thus, the more photons are traced), the quality of the final result improves. Note in particular how the size of the visible circular blotches gets smaller. (1) 10 iterations, (2) 100 iterations, (3) 10,000 iterations. (

*Model courtesy Simon Wendsche.*)

The implementation here uses a `HaltonSampler` to generate camera
paths; doing so ensures that well-distributed samples are used over all of
the iterations in the aggregate.

`nIterations`of SPPM integration>>=

`tile`in image for SPPM>>

`tileBounds`for SPPM tile>>

`tileSampler`for

`pPixel`>>

`SPPMPixel`for

`pPixel`>>

`photonIndex`>>

`photonRay`from light source and initialize

`beta`>>

`grid[h]`>> }

`fr`and direction

`wi`for reflected photon>>

`bsdfSample`for outgoing photon sample>>

`VisiblePoint`in pixel>> }

`L`for SPPM pixel

`pixel`>> image[offset++] = L; } } camera->film->SetImage(image.get()); camera->film->WriteImage(); }

### 16.2.4 Accumulating Visible Points

Similar to the `SamplerIntegrator`, the `SPPMIntegrator`
decomposes the image into tiles of 16 pixels square and parallelizes the
generation of camera paths and visible points over these tiles. The number
of tiles is computed in the same manner as in the fragment
<<Compute number of tiles, `nTiles`, to use for parallel
rendering>>.

Unlike the path tracer, where the BSDF can be discarded after estimating
the direct illumination and sampling the outgoing direction at each vertex,
here we need to store the `BSDF`s for the visible points until the
photon pass for the current iteration is done. Therefore, the
`MemoryArena`s used for allocating the `BSDF`s during camera path
tracing aren’t reset at the end of loop iterations here.

Note also that we only allocate one arena per worker thread used to run the
parallel `for` loop and use the `ThreadIndex` global
to index into the vector, rather than allocating a
separate arena for each loop iteration. In this way, we avoid the overhead
of having many separate `MemoryArena`s while still ensuring that each
arena is not used by more than one processing thread (which would lead to
race conditions in the `MemoryArena` methods).

`tile`in image for SPPM>>

`tileBounds`for SPPM tile>>

`tileSampler`for

`pPixel`>>

`SPPMPixel`for

`pPixel`>>

We also need a unique `Sampler` for the thread processing the tile; as
before, `Sampler::Clone()` provides one.

`tile`in image for SPPM>>=

`tileBounds`for SPPM tile>>

`tileSampler`for

`pPixel`>>

`SPPMPixel`for

`pPixel`>>

The fragment <<Compute `tileBounds` for SPPM tile>> is very
similar to <<Compute sample bounds for tile>> from
Section 1.3.4 and therefore won’t be included here.

Recall that in `SamplerIntegrator`s, the sample vectors in each pixel are
all requested in sequence until the last one is consumed, at which point work
starts on the next pixel. In contrast, the
first SPPM iteration uses the first sample vector for
each pixel; later, the second iteration uses the second sample vector, and
so forth. (In other words, the loop nesting has been interchanged from
“for each pixel, for each sample number,” to “for each sample number,
for each pixel.”) It is for just this use case that the `Sampler`
provides the `SetSampleNumber()` method, which configures the sampler
to provide samples from the given sample vector for the pixel.

`tileSampler`for

`pPixel`>>=

We can now start the path from the camera, following the usual approach. As
with the `PathIntegrator`, the `beta` variable holds the current path
throughput weight .

Path tracing can now proceed. As with most other `Integrator`s, the
path length is limited by a predefined maximum depth.

`SPPMPixel`for

`pPixel`>>

To find the `SPPMPixel` in the `pixels` array for the current
pixel, we need to offset by the minimum pixel coordinate and convert to a
linear index.

`SPPMPixel`for

`pPixel`>>=

As described in Section 16.2.2, a regular direct
lighting calculation is performed at each vertex of the camera path. Thus,
for rays that don’t intersect any scene geometry, infinite area lights must
be allowed a chance to contribute direct lighting via the
`Light::Le()` method.

`SPPMPixel::Ld` records the weighted sum of emitted and reflected
direct illumination for all camera path vertices for the pixel (in other
words, the first two terms of Equation (16.12)). Note
that these terms are also evaluated at vertices found by sampling the BSDF
and tracing a ray for the third term; `Ld` doesn’t just store direct
illumination at the first vertex. Because this sum of outgoing radiance
includes contributions of all of the samples in a pixel, this value must be
divided by `SPPMIntegrator::nIterations` to get the average direct
illumination for the final pixel radiance estimate.

More commonly, the ray intersects a surface and <<Process SPPM camera ray intersection>> executes

First, we need the BSDF at the intersection point. Recall from
Section 11.3 that a `nullptr`-valued `BSDF *` means that
an intersection should be ignored, as the surface intersected is only in the
scene to delineate the boundary of participating media. The
`SPPMIntegrator` does not account for participating media; thus we simply
skip over the intersection and restart the current loop iteration.

As in other integrators, emitted radiance from the surface is only included for the first intersection from the camera or after a specular bounce, where no direct lighting calculation was possible.

The implementation here creates a visible point at the first diffuse surface found or if the path is about to hit its maximum length and we have a glossy surface. As explained earlier, a visible point at a perfectly specular surface will never respond to incident photons, and a glossy surface may have high variance if the specular lobe is tight and not enough photons have arrived to represent the incident radiance distribution well.

The `VisiblePoint` structure records a point found along a camera path
at which we’ll look for nearby photons during the photon shooting pass. It
stores enough information to compute the reflected radiance due to an
incident photon, which in turn is scaled by the path throughput weight of the
visible point to find the overall contribution to the original image
sample.

Sampling a ray leaving a vertex follows the same form as <<Sample
BSDF to get new path direction>> in the `PathIntegrator` and is
therefore not included here.

### 16.2.5 Visible Point Grid Construction

During the photon pass, whenever a photon intersects a surface, we need to efficiently find the visible points where the distance from the visible point to the photon’s intersection is less than the current search radius for the visible point’s pixel. The implementation here uses a uniform grid over the bounding box of all of the visible points; Exercise 16.7 at the end of the chapter has you implement other data structures in place of the grid and investigate trade-offs.

The grid is usually sparsely filled; many voxels have no visible points
inside their extent. (For starters, any voxel in the grid with no surfaces
in its volume will never have a visible point in it.) Therefore, rather
than allocating storage for all of the voxels, the grid is represented by a
hash table where a hash function transforms 3D voxel coordinates to an
index into the `grid` array.

In the following, we’ll construct the grid in parallel using multiple
threads. These threads will use atomic operations to update the grid, so a
`std::atomic` is used for the grid voxel element type here.

Each grid cell stores a linked list, where each node points to an
`SPPMPixel` whose visible point’s search volume overlaps the grid cell. A visible
point may overlap multiple grid cells, so it’s desirable to keep the node representation
compact by only storing a pointer to the `SPPMPixel` rather than making a
copy for each cell it overlaps.

If there’s no visible point for the pixel for the current iteration, the pixel will have a path throughput weight (and no attempt should be made to place a visible point in the grid for that pixel). This case can happen if the path from the camera leaves the scene or is terminated before intersecting a diffuse surface.

Otherwise, the implementation here computes the bounding box centered at the visible point with extent , the current photon search radius for the pixel. In turn, when we later have a photon intersection point, we will only need to consider the visible points for the grid cell that the photon is in to find the visible points that the photon may contribute to (Figure 16.10). Because different visible points will have different search radii depending on how many photons have contributed to their pixel so far, it would otherwise be unwieldy to find the potentially relevant visible points for a photon if the visible points were stored without accounting for the volume of space they will accept photons from.

Given the overall bound of the grid, we need to decide how large the voxels should be and thus how finely to subdivide space. On one hand, if the voxels are too large, the photon shooting pass will be inefficient, as each photon will have to check many visible points to see if it contributes to each one. If they’re too small, then each visible point will overlap many voxels, and the memory needed to represent the grid will be excessive.

Here, we compute an initial resolution such that the voxel width in the largest grid dimension is roughly equal to the largest current search radius of all of the visible points. This limits the maximum number of voxels any visible point can overlap. In turn, resolutions for the other two dimensions of the grid are set so that voxels have roughly the same width in all dimensions.

The visible points can now be added to the grid. Because both the grid and
the `BSDF`s for visible points from the camera path must be kept
around until the end of the photon tracing pass, we can reuse the per-thread
`MemoryArena`s that were created earlier for the `BSDF`s to
allocate `SPPMPixelListNode`s.

Each visible point is added to all of the grid cells that its bounding box overlaps.

`ToGrid()` returns the coordinates of the voxel in the grid that the
given point lies in. The Boolean return value indicates whether the point
`p` is inside the grid’s bounds; if it isn’t, the returned coordinates
`pi` are clamped to be inside the range of valid coordinates.

`hash()` hashes the coordinates of the voxel,
returning an index into the `grid` array defined earlier; it is a
straightforward hash function, and its implementation isn’t included
here.
A new `SPPMPixelListNode` is allocated and the task now
is to add this list node to the head of the linked list in `grid[h]`.

`node`to the start of

`grid[h]`’s linked list>>

Given the grid index, atomic operations can be used to allow multiple
threads to add visible points to the grid concurrently without needing to
hold any locks. (See Appendix A.6.2 for further discussion
of atomic operations.) In the absence of concurrency, we’d just want to set
`node->next` to point to the head of the list in `grid[h]` and
assign `node` to `grid[h]`. That approach will not work if
multiple threads are updating the lists concurrently; it’s possible that the
first assignment will become stale: another thread might modify `grid[h]`
between the current thread initializing `node->next` and then
trying to assign `node` to `grid[h]`.

The `compare_exchange_weak()` method of `std::atomic` addresses
this issue: the first parameter is the value that we expect `grid[h]`
to have, and the second gives the value that we’d like to set it to. It
performs the assignment only if our expectation was correct. Thus, in the
common case, `node->next` and `grid[h]` have the same pointer
value, the assignment occurs, and `true` is returned. The node has
been added to the list.

If the pointer stored in `grid[h]` has been changed by another thread,
then `compare_exchange_weak()` actually updates the first parameter,
`node->next`, with the current value of `grid[h]` before
returning `false`. We are thus all set to try the atomic compare and
exchange again, as `node->next` points to the new head of the list.
This process continues until the assignment is successful.

`node`to the start of

`grid[h]`’s linked list>>=

### 16.2.6 Accumulating Photon Contributions

Given the grid of visible points, the `SPPMIntegrator` can now follow
photon paths through the scene. The total of `photonsPerIteration`
photons to be traced for the current iteration are traced using multiple
threads. A separate `MemoryArena` is available for each worker thread;
this arena is reset after each photon path, so a new pool of per-thread
arenas is allocated rather than reusing the one used for BSDFs and grid linked
list nodes.

`photonIndex`>>

`photonRay`from light source and initialize

`beta`>>

`grid[h]`>> }

`fr`and direction

`wi`for reflected photon>>

`bsdfSample`for outgoing photon sample>>

There’s a balance to strike in choosing the number of photons to trace for each SPPM iteration: too many, and pixels’ radii won’t have a chance to decrease as more photons arrive and too many too-far-away photons are used. Too few, and the overhead of finding visible points and making the grid of them won’t be amortized over enough photons. In practice, a few hundred thousand to a few million per iteration generally works well (see Figure 16.11).

A Halton sequence provides a set of well-distributed sample points for all
of the photon paths over all of the iterations. `haltonIndex` records
the index of the Halton sequence (corresponding to in
Equation (7.7)) for the current photon; it can also be
seen as a global index of photons traced. In other words, it starts at 0
for the first photon and passes through all subsequent integer values
for following photons. It’s important to use a 64-bit integer for this value, since an
unsigned 32-bit integer would overflow after roughly 4 billion photons; many more
photons may be needed for high-quality images.

The dimension of the sample, corresponding to the th prime number in
Equation (7.7), is maintained in `haltonDim`.

`photonIndex`>>=

`photonRay`from light source and initialize

`beta`>>

`grid[h]`>> }

`fr`and direction

`wi`for reflected photon>>

`bsdfSample`for outgoing photon sample>>

Which light to start the path from is determined by sampling from the PDF based on light power computed previously. The first dimension of the Halton sequence is used for the sample value.

The next five dimensions of the sample from the Halton sequence are used for the sample values used to generate the ray leaving the light source.

After the light has been chosen, its `Sample_Le()` method is used to
sample an outgoing ray. Given a light, a ray from the light source is
sampled and its value is initialized based on
Equation (16.8):

where is the probability for sampling this particular light and is the product of the area and directional densities for sampling this particular ray leaving the light. Intersecting this ray against the scene geometry to obtain also samples part of the geometric term except for a cosine factor that must be explicitly integrated into the particle weight .

`photonRay`from light source and initialize

`beta`>>=

Now the integrator can start following the path through the scene, updating after each scattering event. The photon makes no contribution at the first intersection found after it left the light source, since that intersection represents direct illumination, which was already accounted for when tracing paths starting from the camera. For subsequent intersections, illumination is contributed to nearby visible points.

`grid[h]`>> }

`fr`and direction

`wi`for reflected photon>>

`bsdfSample`for outgoing photon sample>>

Given a photon intersection, `ToGrid()`’s return value indicates if
it’s within the extent of the grid. If it isn’t, then by construction,
none of the visible points is interested in this photon’s contribution.
Otherwise, the visible points in the grid cell all need to be checked to
see if the photon is within their radius.

`grid[h]`>> }

Recall that `grid` stores `std::atomic` pointers to
`SPPMPixelListNode`s. Normally, reading from a `std::atomic` value
means that the compiler must be careful to not reorder instructions that read
or write memory around the read of the value of `grid[h]`; this
constraint is necessary so that lock-free algorithms will work as expected.
In this case, the grid has been constructed and no other threads are
concurrently modifying it. Therefore, it’s worthwhile to use the
`std::atomic.load()` method and letting it know that the “relaxed”
memory model, which doesn’t have these constraints, can be used to read the
initial grid pointer. This approach has a significant performance benefit:
for a simple scene of a few hundred triangles (where not too much time is
spent tracing rays), the photon pass runs in 20% less time using this
memory model on a 2015-era CPU.

`grid[h]`>>=

`pixel`and for nearby photon>> }

Given a photon contribution, we need to update the sum for the pixel’s
scattered radiance estimate from Equation (16.14). The total
number of contributing photons in this pass is stored in `M`, and the
sum of the product of BSDF values with particle weights is stored in
`Phi`.

`pixel`and for nearby photon>>=

Each pixel’s and values are stored using atomic variables, which
in turn allows multiple threads to safely concurrently update their
values. Because `pbrt`’s `Spectrum` class doesn’t allow atomic updates,
`Phi` is instead represented with an array of `AtomicFloat`
coefficients for each spectral sample. This representation will require
some manual copying of values to and from `Phi` to
`Spectrum`-typed variables in the following, but we think that this
small awkwardness is preferable to the complexity of, for example, a new
`AtomicSpectrum` type.

Having recorded the photon’s contribution, the integrator needs to choose a new outgoing direction from the intersection point and update the value to account for the effect of scattering. Equation (16.7) shows how to incrementally update the particle weight after a scattering event: given some weight that represents the weight for the th intersection of the th particle history, after a scattering event where a new vertex has been sampled, the weight should be set to be

As with the path-tracing integrator, there are a number of reasons to choose the next vertex in the path by sampling the BSDF’s distribution at the intersection point to get a direction and tracing a ray in that direction rather than directly sampling by area on the scene surfaces. Therefore, we again apply the Jacobian to account for this change in measure, all of the terms in except for a single cancel out, and the expression is

`fr`and direction

`wi`for reflected photon>>

`bsdfSample`for outgoing photon sample>>

As before, a `nullptr`-valued `BSDF *` indicates an intersection that
should be ignored.

Sampling the BSDF to find the scattered photon direction follows the usual model.

`fr`and direction

`wi`for reflected photon>>=

`bsdfSample`for outgoing photon sample>>

The next two dimensions of the Halton sample vector are used for the BSDF sample.

`bsdfSample`for outgoing photon sample>>=

The photon scattering step should be implemented carefully in order to keep the photon weights as similar to each other as possible. A method that gives a distribution of photons where all have exactly equal weights was suggested by Jensen (2001, Section 5.2). First, the reflectance is computed at the intersection point. A random decision is then made whether or not to continue the photon’s path with probability proportional to this reflectance. If the photon continues, its scattered direction is found by sampling from the BSDF’s distribution, but it continues with its weight unchanged except for adjusting the spectral distribution based on the surface’s color. Thus, a surface that reflects very little light will reflect few of the photons that reach it, but those that are scattered will continue on with unchanged contributions and so forth.

This particular approach isn’t possible in `pbrt` due to a subtle
implementation detail (a similar issue held for light source sampling
previously as well): in `pbrt`, the `BxDF` interfaces are written so
that the distribution used for importance sampling BSDFs doesn’t
necessarily have to perfectly match the actual distribution of the function
being sampled. It is all the better if it does, but for many complex
BSDFs exactly sampling from its distribution is difficult or impossible.

Therefore, here we will use an approach that generally leads to a similar result but offers more flexibility: at each intersection point, an outgoing direction is sampled with the BSDF’s sampling distribution, and the photon’s updated weight is computed using Equation (16.16). Then the ratio of the luminance of to the luminance of the photon’s old weight is used to set the probability of continuing the path after applying Russian roulette.

The termination probability is thus set so that if the photon’s weight is significantly decreased at the scattering point, the termination probability will be high and if the photon’s weight is essentially unchanged, the termination probability is low. In particular, the termination probability is chosen in a way such that if the photon continues, after its weight has been adjusted for the possibility of termination, its luminance will be the same as it was before scattering. It is easy to verify this property from the fragment below. (This property actually doesn’t hold for the case where , as can happen when the ratio of the BSDF’s value and the PDF is greater than 1.)

After all of the photons for the iteration have been traced, the estimate of the incident radiance visible in each pixel area can now be updated based on contributions from photons in the current pass.

`VisiblePoint`in pixel>> }

Equation (16.13) gives the rules to update the search radius and other quantities related to the photon estimate.

Note that the number of photons that have contributed to the pixel `N`
is actually stored as a `Float`. This quantity must be treated as
continuously valued, not a discrete integer, for the progressive radiance
estimate to converge to the correct value in the limit.

Before the next SPPM iteration begins, it’s necessary to zero out the
visible point in the pixel so that we don’t attempt to re-use this one if
no visible point and `BDSF *` is found in the next iteration.

`VisiblePoint`in pixel>>=

Most of the <<Periodically store SPPM image in film and write
image>> fragment is a straightforward matter of allocating an image of
`Spectrum` values to pass to `Film::SetImage()` and then
initializing the pixels in the image and before calling `Film::WriteImage()`.
We won’t include that boilerplate here, in order to focus on the last step
of the SPPM algorithm, which combines the direct and indirect radiance
estimates.

As described earlier, the direct lighting estimate needs to be divided by the number of pixel samples (which in turn is how many iterations have completed at this point) to get its average value. The indirect photon term is computed using Equation (16.11)—the two values then just need to be added together.

`L`for SPPM pixel

`pixel`>>=