7.8 Image Reconstruction

Given carefully chosen image samples, we need to convert the samples and their computed radiance values into pixel values for display or storage. According to signal processing theory, we need to do three things to compute final values for each of the pixels in the output image:

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

Because we know that we will be resampling the function upper L overTilde at only the pixel locations, it’s 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 isn’t 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. See, for example, Unser (2000) for a survey of these developments. 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.

While the reconstruction techniques used in pbrt are not directly built on these new approaches, they serve to explain the experience of practitioners that applying perfect reconstruction techniques to samples taken for image synthesis generally does not result in the highest quality images.

Figure 7.38: 2D Image Filtering. To compute a filtered pixel value for the pixel marked with a filled circle located at left-parenthesis x comma y right-parenthesis , all of the image samples inside the box around left-parenthesis x comma y right-parenthesis with extent radius.x and radius.y need to be considered. Each of the image samples left-parenthesis x Subscript i Baseline comma y Subscript i Baseline right-parenthesis , denoted by open circles, is weighted by a 2D filter function, f left-parenthesis x minus x Subscript i Baseline comma y minus y Subscript i Baseline right-parenthesis . The weighted average of all samples is the final pixel value.

To reconstruct pixel values, we will consider the problem of interpolating the samples near a particular pixel. To compute a final value for a pixel upper I left-parenthesis x comma y right-parenthesis , interpolation results in computing a weighted average

upper I left-parenthesis x comma y right-parenthesis equals StartFraction sigma-summation Underscript i Endscripts f left-parenthesis x minus x Subscript i Baseline comma y minus y Subscript i Baseline right-parenthesis w left-parenthesis x Subscript i Baseline comma y Subscript i Baseline right-parenthesis upper L left-parenthesis x Subscript i Baseline comma y Subscript i Baseline right-parenthesis Over sigma-summation Underscript i Endscripts f left-parenthesis x minus x Subscript i Baseline comma y minus y Subscript i Baseline right-parenthesis EndFraction comma
(7.12)

where

  • upper L left-parenthesis x Subscript i Baseline comma y Subscript i Baseline right-parenthesis is the radiance value of the i th sample located at left-parenthesis x Subscript i Baseline comma y Subscript i Baseline right-parenthesis
  • w left-parenthesis x Subscript i Baseline comma y Subscript i Baseline right-parenthesis is the sample contribution weight returned by the Camera. As described in Sections 6.4.7 and 13.6.6, the manner in which these weights are computed determines which radiometric quantity the film measures.
  • f is the filter function.

Figure 7.38 shows a pixel at location left-parenthesis x comma y right-parenthesis that has a pixel filter with extent radius.x in the x direction and radius.y in the y direction. All of the samples inside the box given by the filter extent may contribute to the pixel’s value, depending on the filter function’s value for f left-parenthesis x minus x Subscript i Baseline comma y minus y Subscript i Baseline right-parenthesis .

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 (Gibbs phenomenon), meaning edges in the image have faint replicated copies of the edge in nearby pixels. Furthermore, the sinc filter has infinite support: it doesn’t fall off to zero at a finite distance from its center, so all of 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.

7.8.1 Filter Functions

All filter implementations in pbrt are derived from an abstract Filter class, which provides the interface for the f left-parenthesis x comma y right-parenthesis functions used in filtering; see Equation (7.12). The Film class (described in the Section 7.9) stores a pointer to a Filter and uses it to filter image sample contributions when accumulating them into the final image. (Figure 7.39 shows comparisons of zoomed-in regions of images rendered using a variety of the filters from this section to reconstruct pixel values.) The Filter base class is defined in the files core/filter.h and core/filter.cpp.

Figure 7.39: The pixel reconstruction filter used to convert the image samples into pixel values can have a noticeable effect on the character of the final image. Here, we see blowups of a region of the imperial crown model, filtered with (1) the box filter, (2) Gaussian, and (3) Mitchell–Netravali filter. Note that the Mitchell filter gives the sharpest image, while the Gaussian blurs it. The box is the least desirable, since it allows high-frequency aliasing to leak into the final image. (Note the stair-step pattern along bright gold edges, for example.)

<<Filter Declarations>>= 
class Filter { public: <<Filter Interface>> 
virtual ~Filter(); Filter(const Vector2f &radius) : radius(radius), invRadius(Vector2f(1 / radius.x, 1 / radius.y)) { } virtual Float Evaluate(const Point2f &p) const = 0;
<<Filter Public Data>> 
const Vector2f radius, invRadius;
};

Figure 7.40: The extent of filters in pbrt is specified in terms of their radius from the origin to its cutoff point. The support of a filter is its total non-zero extent, here equal to twice its radius.

All filters are centered at the origin left-parenthesis 0 comma 0 right-parenthesis and define a radius beyond which they have a value of 0; this width may be different in the x and y directions. The constructor takes the radius values and stores them along with their reciprocals, for use by the filter implementations. The filter’s overall extent in each direction (its support) is twice the value of its corresponding radius (Figure 7.40).

<<Filter Interface>>= 
Filter(const Vector2f &radius) : radius(radius), invRadius(Vector2f(1 / radius.x, 1 / radius.y)) { }

<<Filter Public Data>>= 
const Vector2f radius, invRadius;

The sole method that Filter implementations need to provide is Evaluate(). It takes as a parameter a 2D point that gives the position of the sample point relative to the center of the filter. The filter’s value at that point is returned. Code elsewhere in the system will never call the filter function with points outside of the filter’s extent, so filter implementations don’t need to check for this case.

<<Filter Interface>>+= 
virtual Float Evaluate(const Point2f &p) const = 0;

Box Filter

One of the most commonly used filters in graphics is the box filter (and, in fact, when filtering and reconstruction aren’t 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’s just about the worst filter possible. Recall from the discussion in Section 7.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 7.41(a) shows a graph of the box filter, and Figure 7.42 shows the result of using the box filter to reconstruct two 1D functions.

Figure 7.41: Graphs of the (a) box filter and (b) triangle filter. Although neither of these is a particularly good filter, they are both computationally efficient, easy to implement, and good baselines for evaluating other filters.

Figure 7.42: The box filter reconstructing (a) a step function and (b) a sinusoidal function with increasing frequency as x increases. This filter does well with the step function, as expected, but does an extremely poor job with the sinusoidal function.

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 x 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.

<<BoxFilter Declarations>>= 
class BoxFilter : public Filter { public: BoxFilter(const Vector2f &radius) : Filter(radius) { } Float Evaluate(const Point2f &p) const; };

Because the evaluation function won’t be called with left-parenthesis x comma y right-parenthesis values outside of the filter’s extent, it can always return 1 for the filter function’s value.

<<BoxFilter Method Definitions>>= 
Float BoxFilter::Evaluate(const Point2f &p) const { return 1.; }

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 7.41(b) for a graph of the triangle filter.

<<TriangleFilter Declarations>>= 
class TriangleFilter : public Filter { public: TriangleFilter(const Vector2f &radius) : Filter(radius) { } Float Evaluate(const Point2f &p) const; };

Evaluating the triangle filter is simple: the implementation just computes a linear function based on the width of the filter in both the x and y directions.

<<TriangleFilter Method Definitions>>= 
Float TriangleFilter::Evaluate(const Point2f &p) const { return std::max((Float)0, radius.x - std::abs(p.x)) * std::max((Float)0, radius.y - std::abs(p.y)); }

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. The Gaussian’s value at the end of its extent is subtracted from the filter value, in order to make the filter go to 0 at its limit (Figure 7.43). The Gaussian does tend to cause slight blurring of the final image compared to some of the other filters, but this blurring can actually help mask any remaining aliasing in the image.

Figure 7.43: Graphs of (a) the Gaussian filter and (b) the Mitchell filter with upper B equals one-third and upper C equals one-third , each with a width of 2. The Gaussian gives images that tend to be a bit blurry, while the negative lobes of the Mitchell filter help to accentuate and sharpen edges in final images.

<<GaussianFilter Declarations>>= 
class GaussianFilter : public Filter { public: <<GaussianFilter Public Methods>> 
GaussianFilter(const Vector2f &radius, Float alpha) : Filter(radius), alpha(alpha), expX(std::exp(-alpha * radius.x * radius.x)), expY(std::exp(-alpha * radius.y * radius.y)) { } Float Evaluate(const Point2f &p) const;
private: <<GaussianFilter Private Data>> 
const Float alpha; const Float expX, expY;
<<GaussianFilter Utility Functions>> 
Float Gaussian(Float d, Float expv) const { return std::max((Float)0, Float(std::exp(-alpha * d * d) - expv)); }
};

The 1D Gaussian filter function of radius r is

f left-parenthesis x right-parenthesis equals normal e Superscript minus alpha x squared Baseline minus normal e Superscript minus alpha r squared Baseline comma

where alpha controls the rate of falloff of the filter. Smaller values cause a slower falloff, giving a blurrier image. The second term here ensures that the Gaussian goes to 0 at the end of its extent rather than having an abrupt cliff. For efficiency, the constructor precomputes the constant term for normal e Superscript minus alpha r squared in each direction.

<<GaussianFilter Public Methods>>= 
GaussianFilter(const Vector2f &radius, Float alpha) : Filter(radius), alpha(alpha), expX(std::exp(-alpha * radius.x * radius.x)), expY(std::exp(-alpha * radius.y * radius.y)) { }

<<GaussianFilter Private Data>>= 
const Float alpha; const Float expX, expY;

Since a 2D Gaussian function is separable into the product of two 1D Gaussians, the implementation calls the Gaussian() function twice and multiplies the results.

<<GaussianFilter Method Definitions>>= 
Float GaussianFilter::Evaluate(const Point2f &p) const { return Gaussian(p.x, expX) * Gaussian(p.y, expY); }

<<GaussianFilter Utility Functions>>= 
Float Gaussian(Float d, Float expv) const { return std::max((Float)0, Float(std::exp(-alpha * d * d) - expv)); }

Mitchell Filter

Filter design is notoriously difficult, mixing mathematical analysis and perceptual experiments. Mitchell and Netravali (1988) have 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 7.43(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. Also, because the final pixel values can therefore become negative, they will eventually need to be clamped to a legal output range.

Figure 7.44 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 very good job with the sinusoidal function, up until the point where the sampling rate isn’t sufficient to capture the function’s detail.

Figure 7.44: The Mitchell–Netravali Filter Used to Reconstruct the Example Functions. It does a good job with both of these functions, (a) introducing minimal ringing with the step function and (b) accurately representing the sinusoid until aliasing from undersampling starts to dominate.

<<MitchellFilter Declarations>>= 
class MitchellFilter : public Filter { public: <<MitchellFilter Public Methods>> 
MitchellFilter(const Vector2f &radius, Float B, Float C) : Filter(radius), B(B), C(C) { } Float Evaluate(const Point2f &p) const; Float Mitchell1D(Float x) const { x = std::abs(2 * x); if (x > 1) return ((-B - 6*C) * x*x*x + (6*B + 30*C) * x*x + (-12*B - 48*C) * x + (8*B + 24*C)) * (1.f/6.f); else return ((12 - 9*B - 6*C) * x*x*x + (-18 + 12*B + 6*C) * x*x + (6 - 2*B)) * (1.f/6.f); }
private: const Float B, C; };

The Mitchell filter has two parameters called upper B and upper C . Although any values can be used for these parameters, Mitchell and Netravali recommend that they lie along the line upper B plus 2 upper C equals 1 .

<<MitchellFilter Public Methods>>= 
MitchellFilter(const Vector2f &radius, Float B, Float C) : Filter(radius), B(B), C(C) { }

The Mitchell-Netravali filter is the product of 1D filter functions in the x and y directions and is therefore separable, like the Gaussian filter. (In fact, all of the provided filters in pbrt are separable.) Nevertheless, the Filter::Evaluate() interface does not enforce this requirement, giving more flexibility in implementing new filters in the future.

<<MitchellFilter Method Definitions>>= 
Float MitchellFilter::Evaluate(const Point2f &p) const { return Mitchell1D(p.x * invRadius.x) * Mitchell1D(p.y * invRadius.y); }

The 1D function used in the Mitchell filter is an even function defined over the range left-bracket negative 2 comma 2 right-bracket . This function is made by joining a cubic polynomial defined over left-bracket 0 comma 1 right-bracket with another cubic polynomial defined over left-bracket 1 comma 2 right-bracket . This combined polynomial is also reflected around the x equals 0 plane to give the complete function. These polynomials are controlled by the upper B and upper C parameters and are chosen carefully to guarantee upper C Superscript 0 and upper C Superscript 1 continuity at x equals 0 , x equals 1 , and x equals 2 . The polynomials are

StartLayout 1st Row 1st Column Blank 2nd Column f left-parenthesis x right-parenthesis 2nd Row 1st Column Blank 2nd Column equals one-sixth StartLayout Enlarged left-brace 1st Row 1st Column left-parenthesis 12 minus 9 upper B minus 6 upper C right-parenthesis StartAbsoluteValue x EndAbsoluteValue cubed plus left-parenthesis negative 18 plus 12 upper B plus 6 upper C right-parenthesis StartAbsoluteValue x EndAbsoluteValue squared plus left-parenthesis 6 minus 2 upper B right-parenthesis 2nd Column StartAbsoluteValue x EndAbsoluteValue less-than 1 2nd Row 1st Column left-parenthesis negative upper B minus 6 upper C right-parenthesis StartAbsoluteValue x EndAbsoluteValue cubed plus left-parenthesis 6 upper B plus 30 upper C right-parenthesis StartAbsoluteValue x EndAbsoluteValue squared plus left-parenthesis minus 12 upper B minus 48 upper C right-parenthesis StartAbsoluteValue x EndAbsoluteValue 2nd Column 1 less-than-or-equal-to StartAbsoluteValue x EndAbsoluteValue less-than 2 3rd Row 1st Column plus left-parenthesis 8 upper B plus 24 upper C right-parenthesis 4th Row 1st Column 0 2nd Column otherwise period EndLayout EndLayout

<<MitchellFilter Public Methods>>+= 
Float Mitchell1D(Float x) const { x = std::abs(2 * x); if (x > 1) return ((-B - 6*C) * x*x*x + (6*B + 30*C) * x*x + (-12*B - 48*C) * x + (8*B + 24*C)) * (1.f/6.f); else return ((12 - 9*B - 6*C) * x*x*x + (-18 + 12*B + 6*C) * x*x + (6 - 2*B)) * (1.f/6.f); }

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, which is necessary for an implementation with reasonable performance. An additional parameter tau controls how many cycles the sinc function passes through before it is clamped to a value of 0. Figure 7.45 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 tau cycles:

w left-parenthesis x right-parenthesis equals StartFraction sine left-parenthesis pi x slash tau right-parenthesis Over pi x slash tau EndFraction period

Figure 7.45 also shows the filter that we will implement here, which is the product of the sinc function and the windowing function.

Figure 7.45: Graphs of the Sinc Filter. (a) The sinc function, truncated after three cycles (blue line) and the Lanczos windowing function (orange line). (b) The product of these two functions, as implemented in the LanczosSincFilter.

Figure 7.46 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 7.11). The windowed sinc filter also does extremely well at reconstructing the sinusoidal function until prealiasing begins.

Figure 7.46: Results of Using the Windowed Sinc Filter to Reconstruct the Example Functions. Here, tau equals 3 . (a) Like the infinite sinc, it suffers from ringing with the step function, although there is much less ringing in the windowed version. (b) The filter does quite well with the sinusoid, however.

<<Sinc Filter Declarations>>= 
class LanczosSincFilter : public Filter { public: <<LanczosSincFilter Public Methods>> 
LanczosSincFilter(const Vector2f &radius, Float tau) : Filter(radius), tau(tau) { } Float Evaluate(const Point2f &p) const; Float Sinc(Float x) const { x = std::abs(x); if (x < 1e-5) return 1; return std::sin(Pi * x) / (Pi * x); } Float WindowedSinc(Float x, Float radius) const { x = std::abs(x); if (x > radius) return 0; Float lanczos = Sinc(x / tau); return Sinc(x) * lanczos; }
private: const Float tau; };

<<LanczosSincFilter Public Methods>>= 
LanczosSincFilter(const Vector2f &radius, Float tau) : Filter(radius), tau(tau) { }

<<Sinc Filter Method Definitions>>= 
Float LanczosSincFilter::Evaluate(const Point2f &p) const { return WindowedSinc(p.x, radius.x) * WindowedSinc(p.y, radius.y); }

The implementation computes the value of the sinc function and then multiplies it by the value of the Lanczos windowing function.

<<LanczosSincFilter Public Methods>>+=  
Float Sinc(Float x) const { x = std::abs(x); if (x < 1e-5) return 1; return std::sin(Pi * x) / (Pi * x); }

<<LanczosSincFilter Public Methods>>+= 
Float WindowedSinc(Float x, Float radius) const { x = std::abs(x); if (x > radius) return 0; Float lanczos = Sinc(x / tau); return Sinc(x) * lanczos; }