8.8 Image Reconstruction

As discussed in Section 5.4.3, each pixel in the Film computes an estimate of the integral of the product of a filter function with samples taken from the image function. In Section 8.1, we saw that sampling theory provides a mathematical foundation for how this filtering operation should be performed in order to achieve an antialiased result. We should, in principle:

1. Reconstruct a continuous image function from the set of image samples.
2. Prefilter that function to remove any frequencies past the Nyquist limit for the pixel spacing.
3. Sample the prefiltered function at the pixel locations to compute the final pixel values.

Because we know that we will be resampling the function at only the pixel locations, it is not necessary to construct an explicit representation of the function. Instead, we can combine the first two steps using a single filter function.

Recall that if the original function had been uniformly sampled at a frequency greater than the Nyquist frequency and reconstructed with the sinc filter, then the reconstructed function in the first step would match the original image function perfectly—quite a feat since we only have point samples. But because the image function almost always will have higher frequencies than could be accounted for by the sampling rate (due to edges, etc.), we chose to sample it nonuniformly, trading off noise for aliasing.

The theory behind ideal reconstruction depends on the samples being uniformly spaced. While a number of attempts have been made to extend the theory to nonuniform sampling, there is not yet an accepted approach to this problem. Furthermore, because the sampling rate is known to be insufficient to capture the function, perfect reconstruction is not possible. Recent research in the field of sampling theory has revisited the issue of reconstruction with the explicit acknowledgment that perfect reconstruction is not generally attainable in practice. This slight shift in perspective has led to powerful new reconstruction techniques. In particular, the goal of research in reconstruction theory has shifted from perfect reconstruction to developing reconstruction techniques that can be shown to minimize error between the reconstructed function and the original function, regardless of whether the original was band limited.

The sinc filter is not an appropriate choice here: recall that the ideal sinc filter is prone to ringing when the underlying function has frequencies beyond the Nyquist limit, meaning edges in the image have faint replicated copies of the edge in nearby pixels (the Gibbs phenomenon; see Section 8.1.5). Furthermore, the sinc filter has infinite support: it does not fall off to zero at a finite distance from its center, so all the image samples would need to be filtered for each output pixel. In practice, there is no single best filter function. Choosing the best one for a particular scene takes a mixture of quantitative evaluation and qualitative judgment. pbrt therefore provides a variety of choices.

Figure 8.48 shows comparisons of zoomed-in regions of images rendered using a variety of the filters from this section to reconstruct pixel values.

8.8.1 Filter Interface

The Filter class defines the interface for pixel reconstruction filters in pbrt. It is defined in the file base/filter.h.

All filters are 2D functions that are centered at the origin and define a radius beyond which they have a value of 0. The radii are different in the and directions but are assumed to be symmetric in each. A filter provides its radius via the Radius() method. The filter’s overall extent in each direction (its support) is twice the value of its corresponding radius (Figure 8.49).

Filter implementations must also provide a method that evaluates their filter function. This function may be called with points that are outside of the filter’s radius; in this case it is the responsibility of the implementation to detect this case and return the value 0. It is not required for the filter values returned by Evaluate() to integrate to 1 since the estimator used to compute pixel values, Equation (5.13), is self-normalizing.

Filters also must be able to return their integral. Most are able to compute this value in closed form. Thus, if calling code requires a normalized filter function, it is easy enough to find it by dividing values returned by Evaluate() by the integral.

Filters must also provide an importance sampling method, Sample, which takes a random sample u in .

The returned FilterSample structure stores both the sampled position p and a weight, which is the ratio of the value of the filter function at p to the value of the PDF used for sampling there. Because some filters are able to exactly sample from their distribution, returning this ratio directly allows them to save the trouble of evaluating those two values and instead to always return a weight of 1.

Given the specification of this interface, we can now implement the GetCameraSample() function that most integrators use to compute the CameraSamples that they pass to the Camera::GenerateRay() methods.

One subtlety in the definition of this function is that it is templated based on the type of the sampler passed to it. If a value of type Sampler is passed to this method, then it proceeds using pbrt’s usual dynamic dispatch mechanism to call the corresponding methods of the Sampler implementation. However, if a concrete sampler type (e.g., HaltonSampler) is passed to it, then the corresponding methods can be called directly (and are generally expanded inline in the function). This capability is used to improve performance in pbrt’s GPU rendering path; see Section 15.3.3.

After the filter’s Sample() method has returned a FilterSample, the image sample position can be found by adding the filter’s sampled offset to the pixel coordinate before a shift of in each dimension accounts for the mapping from discrete to continuous pixel coordinates (recall Section 8.1.4). The filter’s weight is passed along in the CameraSample so that it is available to the Film when its AddSample() method is called.

8.8.2 FilterSampler

Not all Filters are able to easily sample from the distributions of their filter functions. Therefore, pbrt provides a FilterSampler class that wraps up the details of sampling based on a tabularized representation of the filter.

Only the Filter and an allocator are provided to the constructor. We have not found it particularly useful to allow the caller to specify the rate at which the filter function is sampled to construct the table used for sampling, so instead hardcode a sampling rate of 32 times per unit filter extent in each dimension.

domain gives the bounds of the filter and f stores tabularized filter function values.

All the filters currently implemented in pbrt are symmetric about the origin, which means that they could be tabularized over a single quadrant. Further, they are all separable into the product of two 1D functions. Either of these properties could be exploited to reduce the amount of storage required for the tables used for sampling. However, to allow full flexibility with the definition of additional filter functions, the FilterSampler simply evaluates the filter at equally spaced positions over its entire domain to initialize the f array.

Given a tabularized function, it is easy to initialize the sampling distribution.

There are two important details in the implementation of its Sample() method. First, the implementation does not use Filter::Evaluate() to evaluate the filter function but instead uses the tabularized version of it in f. By using the piecewise constant approximation of the filter function, it ensures that the returned weight for a sampled point is always for a constant . If it did not do this, there would be variation in the returned weight for non-constant filter functions, due to the sampling distribution not being exactly proportional to the filter function—see Figure 8.50, which illustrates the issue.

A second important detail is that the integer coordinates of the sample returned by PiecewiseConstant2D::Sample() are used to index into f for filter function evaluation. If instead the point p was scaled up by the size of the f array in each dimension and converted to an integer, the result would occasionally differ from the integer coordinates computed during sampling by PiecewiseConstant2D due to floating-point round-off error. (Using the notation of Section 6.8.1, the issue is that with floating-point arithmetic, .) Again, variance would result, as the ratio might not be .

8.8.3 Box Filter

One of the most commonly used filters in graphics is the box filter (and, in fact, when filtering and reconstruction are not addressed explicitly, the box filter is the de facto result). The box filter equally weights all samples within a square region of the image. Although computationally efficient, it is just about the worst filter possible. Recall from the discussion in Section 8.1.2 that the box filter allows high-frequency sample data to leak into the reconstructed values. This causes postaliasing—even if the original sample values were at a high enough frequency to avoid aliasing, errors are introduced by poor filtering.

Figure 8.51(a) shows a graph of the box filter, and Figure 8.52 shows the result of using the box filter to reconstruct two 1D functions.

For the step function we used previously to illustrate the Gibbs phenomenon, the box does reasonably well. However, the results are much worse for a sinusoidal function that has increasing frequency along the axis. Not only does the box filter do a poor job of reconstructing the function when the frequency is low, giving a discontinuous result even though the original function was smooth, but it also does an extremely poor job of reconstruction as the function’s frequency approaches and passes the Nyquist limit.

For this filter and all the following ones, we will not include the rudimentary constructors and Radius() method implementations.

Evaluating the box filter requires checking that the given point is inside the box.

Sampling is also easy: the random sample u is used to linearly interpolate within the filter’s extent. Since sampling is exact and the filter function is positive, the weight is always 1.

Finally, the integral is equal to the filter’s area.

8.8.4 Triangle Filter

The triangle filter gives slightly better results than the box: the weight falls off linearly from the filter center over the square extent of the filter. See Figure 8.51(b) for a graph of the triangle filter.

Evaluating the triangle filter is simple: it is the product of two linear functions that go to 0 after the width of the filter in both the and directions. Here we have defined the filter to have a slope of , though the filter could alternatively have been defined to have a value of 1 at the origin and a slope that depends on the radius.

Because the filter is separable, its PDF is as well, and so each dimension can be sampled independently. The sampling method uses a separate SampleTent() utility function that is defined in Section A.4.1. Once again, the weight returned in the FilterSample is always 1 because the filter is positive and sampling is exact.

Finally, the triangle filter is easily integrated.

8.8.5 Gaussian Filter

Unlike the box and triangle filters, the Gaussian filter gives a reasonably good result in practice. This filter applies a Gaussian bump that is centered at the pixel and radially symmetric around it. Figure 8.53 compares plots of the Gaussian filter and the Mitchell filter (described in Section 8.8.6). The Gaussian does tend to cause slight blurring of the final image compared to some of the other filters, but this blurring can help mask any remaining aliasing. This filter is the default one used in pbrt.

The Gaussian function is parameterized by the position of the peak and the standard deviation :

Larger values of cause a slower falloff, which leads to a blurrier image when the Gaussian is used as a filter.

The GaussianFilter is centered at the origin, so . Further, the filter function subtracts the value of the Gaussian at the end of its extent from the filter value in order to make the filter go to 0 at its limit:

For efficiency, the constructor precomputes the constant term for in each direction.

The product of the two 1D Gaussian functions gives the overall filter value according to Equation (8.26). The calls to std::max() ensure that the value of 0 is returned for points outside of the filter’s extent.

The integral of the Gaussian is

where is the error function. GaussianIntegral() evaluates its value over a given range. The filter function’s integral can be computed by evaluating the Gaussian’s integral over the filter’s range and subtracting the integral of the offset that takes the filter to zero at the end of its extent.

It is possible to sample from the Gaussian function using a polynomial approximation to the inverse error function, though that is not sufficient in this case, given the presence of the second term of the filter function in Equation (8.26). pbrt’s GaussianFilter implementation therefore uses a FilterSampler for sampling.

8.8.6 Mitchell Filter

Filter design is notoriously difficult, mixing mathematical analysis and perceptual experiments. Mitchell and Netravali (1988) developed a family of parameterized filter functions in order to be able to explore this space in a systematic manner. After analyzing test subjects’ subjective responses to images filtered with a variety of parameter values, they developed a filter that tends to do a good job of trading off between ringing (phantom edges next to actual edges in the image) and blurring (excessively blurred results)—two common artifacts from poor reconstruction filters.

Note from the graph in Figure 8.53(b) that this filter function takes on negative values out by its edges; it has negative lobes. In practice these negative regions improve the sharpness of edges, giving crisper images (reduced blurring). If they become too large, however, ringing tends to start to enter the image. Furthermore, because the final pixel values can become negative, they will eventually need to be clamped to a legal output range.

Figure 8.54 shows this filter reconstructing the two test functions. It does extremely well with both of them: there is minimal ringing with the step function, and it does a good job with the sinusoidal function, up until the point where the sampling rate is not sufficient to capture the function’s detail.

The Mitchell filter has two parameters called and . Although any values can be used for these parameters, Mitchell and Netravali recommend that they lie along the line .

The Mitchell–Netravali filter is the product of 1D filter functions in the and directions and is therefore separable.

The 1D function used in the Mitchell filter is an even function defined over the range . This function is made by joining a cubic polynomial defined over with another cubic polynomial defined over . This combined polynomial is also reflected around the plane to give the complete function. These polynomials are controlled by the and parameters and are chosen carefully to guarantee and continuity at , , and . The polynomials are

Mitchell1D() evaluates this function. Its implementation is straightforward and is not included here.

As a cubic polynomial, sampling this filter function directly would require inverting a quartic. Therefore, the MitchellFilter uses the FilterSampler for sampling.

However, the function is easily integrated. The result is independent of the values of and .

8.8.7 Windowed Sinc Filter

Finally, the LanczosSincFilter class implements a filter based on the sinc function. In practice, the sinc filter is often multiplied by another function that goes to 0 after some distance. This gives a filter function with finite extent. An additional parameter controls how many cycles the sinc function passes through before it is clamped to a value of 0. Figure 8.55 shows a graph of three cycles of the sinc function, along with a graph of the windowing function we use, which was developed by Lanczos. The Lanczos window is just the central lobe of the sinc function, scaled to cover the cycles:

Figure 8.55 also shows the filter that we will implement here, which is the product of the sinc function and the windowing function. It is evaluated by the WindowedSinc() utility function.

Its implementation uses the Sinc() function, which in turn is implemented using the numerically robust SinXOverX() function.

Figure 8.56 shows the windowed sinc’s reconstruction results for uniform 1D samples. Thanks to the windowing, the reconstructed step function exhibits far less ringing than the reconstruction using the infinite-extent sinc function (compare to Figure 8.11). The windowed sinc filter also does extremely well at reconstructing the sinusoidal function until prealiasing begins.

The evaluation method is easily implemented in terms of the WindowedSinc() function.

There is no convenient closed-form approach for sampling from the windowed sinc function’s distribution, so a FilterSampler is used here as well.

There is no closed-form expression of the filter’s integral, so its Integral() method, not included in the text, approximates it using a Riemann sum.