10.4 Image Texture

Image textures store 2D arrays of point-sampled values of a texture function. They use these samples to reconstruct a continuous image function that can be evaluated at an arbitrary left-parenthesis s comma t right-parenthesis position. These sample values are often called texels, since they are similar to pixels in an image but are used in the context of a texture. Image textures are the most widely used type of texture in computer graphics; digital photographs, scanned artwork, images created with image-editing programs, and images generated by renderers are all extremely useful sources of data for this particular texture representation (Figure 10.14).

Figure 10.14: An Example of Image Textures. Image textures are used throughout the Watercolor scene to represent spatially varying surface appearance properties. (a) Scene rendered with image textures. (b) Each image texture has been replaced with its average value. Note how much visual richness is lost. (Scene courtesy of Angelo Ferretti.)

As with most of the other types of texture, pbrt provides both Float and spectral variants. Both implementations inherit from ImageTextureBase, which provides some common functionality.

<<ImageTextureBase Definition>>= 
class ImageTextureBase { public: <<ImageTextureBase Public Methods>> 
ImageTextureBase(TextureMapping2D mapping, std::string filename, MIPMapFilterOptions filterOptions, WrapMode wrapMode, Float scale, bool invert, ColorEncoding encoding, Allocator alloc) : mapping(mapping), filename(filename), scale(scale), invert(invert) { <<Get MIPMap from texture cache if present>> 
TexInfo texInfo(filename, filterOptions, wrapMode, encoding); std::unique_lock<std::mutex> lock(textureCacheMutex); if (auto iter = textureCache.find(texInfo); iter != textureCache.end()) { mipmap = iter->second; return; } lock.unlock();
<<Create MIPMap for filename and add to texture cache>> 
mipmap = MIPMap::CreateFromFile(filename, filterOptions, wrapMode, encoding, alloc); lock.lock(); textureCache[texInfo] = mipmap;
} static void ClearCache() { textureCache.clear(); } void MultiplyScale(Float s) { scale *= s; }
protected: <<ImageTextureBase Protected Members>> 
TextureMapping2D mapping; std::string filename; Float scale; bool invert; MIPMap *mipmap;
private: <<ImageTextureBase Private Members>> 
static std::mutex textureCacheMutex; static std::map<TexInfo, MIPMap *> textureCache;
};

In the following, we will present the implementation of SpectrumImageTexture; FloatImageTexture is analogous and does not add anything new.

<<SpectrumImageTexture Definition>>= 
class SpectrumImageTexture : public ImageTextureBase { public: <<SpectrumImageTexture Public Methods>> 
SpectrumImageTexture(TextureMapping2D mapping, std::string filename, MIPMapFilterOptions filterOptions, WrapMode wrapMode, Float scale, bool invert, ColorEncoding encoding, SpectrumType spectrumType, Allocator alloc) : ImageTextureBase(mapping, filename, filterOptions, wrapMode, scale, invert, encoding, alloc), spectrumType(spectrumType) {} PBRT_CPU_GPU SampledSpectrum Evaluate(TextureEvalContext ctx, SampledWavelengths lambda) const; static SpectrumImageTexture *Create(const Transform &renderFromTexture, const TextureParameterDictionary &parameters, SpectrumType spectrumType, const FileLoc *loc, Allocator alloc); std::string ToString() const;
private: <<SpectrumImageTexture Private Members>> 
SpectrumType spectrumType;
};

10.4.1 Texture Memory Management

The caller of SpectrumImageTexture’s constructor provides a texture mapping function, the filename of an image, various parameters that control the filtering of the image map, how boundary conditions are managed, and how colors are converted to spectral samples. All the necessary initialization is handled by ImageTextureBase.

<<SpectrumImageTexture Public Methods>>= 
SpectrumImageTexture(TextureMapping2D mapping, std::string filename, MIPMapFilterOptions filterOptions, WrapMode wrapMode, Float scale, bool invert, ColorEncoding encoding, SpectrumType spectrumType, Allocator alloc) : ImageTextureBase(mapping, filename, filterOptions, wrapMode, scale, invert, encoding, alloc), spectrumType(spectrumType) {}

As was discussed in Section 4.6.6, RGB colors are transformed into spectra differently depending on whether or not they represent reflectances. The spectrumType records what type of RGB a texture represents.

<<SpectrumImageTexture Private Members>>= 
SpectrumType spectrumType;

The contents of the image file are used to create an instance of the MIPMap class that stores the texels in memory and handles the details of reconstruction and filtering to reduce aliasing.

<<ImageTextureBase Public Methods>>= 
ImageTextureBase(TextureMapping2D mapping, std::string filename, MIPMapFilterOptions filterOptions, WrapMode wrapMode, Float scale, bool invert, ColorEncoding encoding, Allocator alloc) : mapping(mapping), filename(filename), scale(scale), invert(invert) { <<Get MIPMap from texture cache if present>> 
TexInfo texInfo(filename, filterOptions, wrapMode, encoding); std::unique_lock<std::mutex> lock(textureCacheMutex); if (auto iter = textureCache.find(texInfo); iter != textureCache.end()) { mipmap = iter->second; return; } lock.unlock();
<<Create MIPMap for filename and add to texture cache>> 
mipmap = MIPMap::CreateFromFile(filename, filterOptions, wrapMode, encoding, alloc); lock.lock(); textureCache[texInfo] = mipmap;
}

A floating-point scale can be specified with each texture; it is applied to the values returned by the Evaluate() method. Further, a true value for the invert parameter causes the texture value to be subtracted from 1 before it is returned. While the same functionality can be achieved with scale and mix textures, it is easy to also provide that functionality directly in the texture here. Doing so can lead to more efficient texture evaluation on GPUs, as is discussed further in Section 15.3.9.

<<ImageTextureBase Protected Members>>= 
TextureMapping2D mapping; std::string filename; Float scale; bool invert; MIPMap *mipmap;

Each MIP map may require a meaningful amount of memory, and a complex scene may have thousands of image textures. Because an on-disk image may be reused for multiple textures in a scene, pbrt maintains a table of MIP maps that have been loaded so far so that they are only loaded into memory once even if they are used in more than one image texture.

pbrt loads textures in parallel after the scene description has been parsed; doing so reduces startup time before rendering begins. Therefore, a mutex is used here to ensure that only one thread accesses the texture cache at a time. Note that if the MIPMap is not found in the cache, the lock is released before it is read so that other threads can access the cache in the meantime.

<<Get MIPMap from texture cache if present>>= 
TexInfo texInfo(filename, filterOptions, wrapMode, encoding); std::unique_lock<std::mutex> lock(textureCacheMutex); if (auto iter = textureCache.find(texInfo); iter != textureCache.end()) { mipmap = iter->second; return; } lock.unlock();

The texture cache itself is managed with a std::map.

<<ImageTextureBase Private Members>>= 
static std::mutex textureCacheMutex; static std::map<TexInfo, MIPMap *> textureCache;

TexInfo is a simple structure that acts as a key for the texture cache std::map. It holds all the specifics that must match for a MIPMap to be reused in another image texture.

<<TexInfo Definition>>= 
struct TexInfo { <<TexInfo Public Methods>> 
TexInfo(const std::string &f, MIPMapFilterOptions filterOptions, WrapMode wm, ColorEncoding encoding) : filename(f), filterOptions(filterOptions), wrapMode(wm), encoding(encoding) {} bool operator<(const TexInfo &t) const { return std::tie(filename, filterOptions, encoding, wrapMode) < std::tie(t.filename, t.filterOptions, t.encoding, t.wrapMode); } std::string ToString() const;
std::string filename; MIPMapFilterOptions filterOptions; WrapMode wrapMode; ColorEncoding encoding; };

The TexInfo constructor, not included here, sets its member variables with provided values. Its only other method is a comparison operator, which is required by std::map.

<<TexInfo Public Methods>>= 
bool operator<(const TexInfo &t) const { return std::tie(filename, filterOptions, encoding, wrapMode) < std::tie(t.filename, t.filterOptions, t.encoding, t.wrapMode); }

If the texture has not yet been loaded, a call to CreateFromFile() yields a MIPMap for it. If the file is not found or there is an error reading it, pbrt exits with an error message, so a nullptr return value does not need to be handled here.

<<Create MIPMap for filename and add to texture cache>>= 
mipmap = MIPMap::CreateFromFile(filename, filterOptions, wrapMode, encoding, alloc); lock.lock(); textureCache[texInfo] = mipmap;

10.4.2 Image Texture Evaluation

Before describing the MIPMap implementation, we will discuss the SpectrumImageTexture Evaluate() method.

<<SpectrumImageTexture Method Definitions>>= 
SampledSpectrum SpectrumImageTexture::Evaluate( TextureEvalContext ctx, SampledWavelengths lambda) const { <<Apply texture mapping and flip t coordinate for image texture lookup>> 
TexCoord2D c = mapping.Map(ctx); c.st[1] = 1 - c.st[1];
<<Lookup filtered RGB value in MIPMap>> 
RGB rgb = scale * mipmap->Filter<RGB>(c.st, {c.dsdx, c.dtdx}, {c.dsdy, c.dtdy}); rgb = ClampZero(invert ? (RGB(1, 1, 1) - rgb) : rgb);
<<Return SampledSpectrum for RGB image texture value>> 
if (const RGBColorSpace *cs = mipmap->GetRGBColorSpace(); cs) { if (spectrumType == SpectrumType::Unbounded) return RGBUnboundedSpectrum(*cs, rgb).Sample(lambda); else if (spectrumType == SpectrumType::Albedo) return RGBAlbedoSpectrum(*cs, Clamp(rgb, 0, 1)).Sample(lambda); else return RGBIlluminantSpectrum(*cs, rgb).Sample(lambda); } DCHECK(rgb[0] == rgb[1] && rgb[1] == rgb[2]); return SampledSpectrum(rgb[0]);
}

It is easy to compute the left-parenthesis s comma t right-parenthesis texture coordinates and their derivatives for filtering with the TextureMapping2D’s Map() method. However, the t coordinate must be flipped, because pbrt’s Image class (and in turn, MIPMap, which is based on it) defines left-parenthesis 0 comma 0 right-parenthesis to be the upper left corner of the image, while image textures have left-parenthesis 0 comma 0 right-parenthesis at the lower left. (These are the typical conventions for indexing these entities in computer graphics.)

<<Apply texture mapping and flip t coordinate for image texture lookup>>= 
TexCoord2D c = mapping.Map(ctx); c.st[1] = 1 - c.st[1];

The MIPMap’s Filter() method provides the filtered value of the image texture over the specified region; any specified scale or inversion is easily applied to the value it returns. A call to ClampZero() here ensures that no negative values are returned after inversion.

<<Lookup filtered RGB value in MIPMap>>= 
RGB rgb = scale * mipmap->Filter<RGB>(c.st, {c.dsdx, c.dtdx}, {c.dsdy, c.dtdy}); rgb = ClampZero(invert ? (RGB(1, 1, 1) - rgb) : rgb);

As discussed in Section 4.6.2, an RGB color space is necessary in order to interpret the meaning of an RGB color value. Normally, the code that reads image file formats from disk returns an RGBCcolorSpace with the read image. Most RGB image formats default to sRGB, and some allow specifying an alternative color space. (For example, OpenEXR allows specifying the primaries of an arbitrary RGB color space in the image file’s metadata.) A color space and the value of spectrumType make it possible to create the appropriate type RGB spectrum, and in turn, its Spectrum::Sample() can be called to get the SampledSpectrum that will be returned.

If the MIPMap has no associated color space, the image is assumed to have the same value in all channels and a constant value is returned for all the spectrum samples. This assumption is verified by a DCHECK() call in non-optimized builds.

<<Return SampledSpectrum for RGB image texture value>>= 
if (const RGBColorSpace *cs = mipmap->GetRGBColorSpace(); cs) { if (spectrumType == SpectrumType::Unbounded) return RGBUnboundedSpectrum(*cs, rgb).Sample(lambda); else if (spectrumType == SpectrumType::Albedo) return RGBAlbedoSpectrum(*cs, Clamp(rgb, 0, 1)).Sample(lambda); else return RGBIlluminantSpectrum(*cs, rgb).Sample(lambda); } DCHECK(rgb[0] == rgb[1] && rgb[1] == rgb[2]); return SampledSpectrum(rgb[0]);

10.4.3 MIP Maps

As always, if the image texture function has higher frequency detail than can be represented by the texture sampling rate, aliasing will be present in the final image. Any frequencies higher than the Nyquist limit must be removed by prefiltering before the function is evaluated. Figure 10.15 shows the basic problem we face: an image texture has texels that are samples of some image function at a fixed frequency. The filter region for the lookup is given by its left-parenthesis s comma t right-parenthesis center point and offsets to the estimated texture coordinate locations for the adjacent image samples. Because these offsets are estimates of the texture sampling rate, we must remove any frequencies higher than twice the distance to the adjacent samples in order to satisfy the Nyquist criterion.

Figure 10.15: Given a point at which to perform an image map lookup (denoted by the solid point in the center) and estimates of the texture-space sampling rate (denoted by adjacent solid points), it may be necessary to filter the contributions of a large number of texels in the image map (denoted by open points).

The texture sampling and reconstruction process has a few key differences from the image sampling process discussed in Chapter 8. These differences make it possible to address the antialiasing problem with more effective and less computationally expensive techniques. For example, here it is inexpensive to get the value of a sample—only an array lookup is necessary (as opposed to having to trace a number of rays to compute radiance). Further, because the texture image function is fully defined by the set of samples and there is no mystery about what its highest frequency could be, there is no uncertainty related to the function’s behavior between samples. These differences make it possible to remove detail from the texture before sampling, thus eliminating aliasing.

However, the texture sampling rate will typically change from pixel to pixel. The sampling rate is determined by scene geometry and its orientation, the texture coordinate mapping function, and the camera projection and image sampling rate. Because the texture sampling rate is not fixed, texture filtering algorithms need to be able to filter over arbitrary regions of texture samples efficiently.

The MIPMap class implements a number of methods for texture filtering with spatially varying filter widths. It can be found in the files util/mipmap.h and util/mipmap.cpp. The filtering algorithms it offers range from simple point sampling to bilinear interpolation and trilinear interpolation, which is fast and easy to implement and was widely used for texture filtering in early graphics hardware, to elliptically weighted averaging, which is more complex but returns extremely high-quality results. Figure 10.16 compares the result of texture filtering using trilinear interpolation and the EWA algorithm.

<<MIPMap Definition>>= 
class MIPMap { public: <<MIPMap Public Methods>> 
MIPMap(Image image, const RGBColorSpace *colorSpace, WrapMode wrapMode, Allocator alloc, const MIPMapFilterOptions &options); static MIPMap *CreateFromFile(const std::string &filename, const MIPMapFilterOptions &options, WrapMode wrapMode, ColorEncoding encoding, Allocator alloc); template <typename T> T Filter(Point2f st, Vector2f dstdx, Vector2f dstdy) const; std::string ToString() const; Point2i LevelResolution(int level) const { return pyramid[level].Resolution(); } int Levels() const { return int(pyramid.size()); } const RGBColorSpace *GetRGBColorSpace() const { return colorSpace; } const Image &GetLevel(int level) const { return pyramid[level]; }
private: <<MIPMap Private Methods>> 
template <typename T> T Texel(int level, Point2i st) const; template <typename T> T Bilerp(int level, Point2f st) const; template <typename T> T EWA(int level, Point2f st, Vector2f dst0, Vector2f dst1) const;
<<MIPMap Private Members>> 
pstd::vector<Image> pyramid; const RGBColorSpace *colorSpace; WrapMode wrapMode; MIPMapFilterOptions options;
};

Figure 10.16: Filtering the image map properly substantially improves the image. Trilinear interpolation is used for the sphere on the left and the EWA algorithm is used for the sphere on the right. Both of these approaches give a much better result than the unfiltered image map on the sphere on the left in Figure 10.1. Trilinear interpolation is not as effective as EWA at handling strongly anisotropic filter footprints, which is why the lines on the left sphere are blurrier. In regions with highly anisotropic filter footprints such as the pole of the sphere and toward the edges, EWA resolves much more detail.

If an RGB image is provided to the MIPMap constructor, its channels should be stored in R, G, B order in memory; for efficiency, the following code assumes that this is the case. All the code that currently uses MIPMaps in pbrt ensures that this is so.

<<MIPMap Method Definitions>>= 
MIPMap::MIPMap( Image image, const RGBColorSpace *colorSpace, WrapMode wrapMode, Allocator alloc, const MIPMapFilterOptions &options) : colorSpace(colorSpace), wrapMode(wrapMode), options(options) { pyramid = Image::GeneratePyramid(std::move(image), wrapMode, alloc); }

To limit the potential number of texels that need to be accessed, these filtering methods use an image pyramid of increasingly lower resolution prefiltered versions of the original image to accelerate their operation. The original image texels are at the bottom level of the pyramid, and the image at each level is half the resolution of the previous level, up to the top level, which has a single texel representing the average of all the texels in the original image. This collection of images needs at most 1 slash 3 more memory than storing the most detailed level alone and can be used to quickly find filtered values over large regions of the original image. The basic idea behind the pyramid is that if a large area of texels needs to be filtered, a reasonable approximation is to use a higher level of the pyramid and do the filtering over the same area there, accessing many fewer texels.

The MIPMap’s image pyramid is represented by a vector of Images. See Section B.5 for the implementation of Image and Section B.5.5 for its GeneratePyramid() method, which generates image pyramids.

<<MIPMap Private Members>>= 
pstd::vector<Image> pyramid; const RGBColorSpace *colorSpace; WrapMode wrapMode; MIPMapFilterOptions options;

The choice of filtering algorithm and a parameter used by the EWA method are represented by MIPMapFilterOptions.

<<MIPMapFilterOptions Definition>>= 
struct MIPMapFilterOptions { FilterFunction filter = FilterFunction::EWA; Float maxAnisotropy = 8.f; };

<<FilterFunction Definition>>= 
enum class FilterFunction { Point, Bilinear, Trilinear, EWA };

A few simple utility methods return information about the image pyramid and the MIPMap’s color space.

<<MIPMap Public Methods>>= 
Point2i LevelResolution(int level) const { return pyramid[level].Resolution(); } int Levels() const { return int(pyramid.size()); } const RGBColorSpace *GetRGBColorSpace() const { return colorSpace; } const Image &GetLevel(int level) const { return pyramid[level]; }

Given the image pyramid, we will define some utility MIPMap methods that retrieve the texel value at a specified pyramid level and discrete integer pixel coordinates. For the RGB variant, there is an implicit assumption that the image channels are laid out in R, G, B (and maybe A) order.

<<MIPMap Method Definitions>>+=  
template <> RGB MIPMap::Texel(int level, Point2i st) const { if (int nc = pyramid[level].NChannels(); nc == 3 || nc == 4) return RGB(pyramid[level].GetChannel(st, 0, wrapMode), pyramid[level].GetChannel(st, 1, wrapMode), pyramid[level].GetChannel(st, 2, wrapMode)); else { Float v = pyramid[level].GetChannel(st, 0, wrapMode); return RGB(v, v, v); } }

The Float specialization of Texel(), not included here, is analogous.

10.4.4 Image Map Filtering

The MIPMap Filter() method returns a filtered image function value at the provided left-parenthesis s comma t right-parenthesis coordinates. It takes two derivatives that give the change in left-parenthesis s comma t right-parenthesis with respect to image pixel samples.

<<MIPMap Method Definitions>>+=  
template <typename T> T MIPMap::Filter(Point2f st, Vector2f dst0, Vector2f dst1) const { if (options.filter != FilterFunction::EWA) { <<Handle non-EWA MIP Map filter>> 
Float width = 2 * std::max({std::abs(dst0[0]), std::abs(dst0[1]), std::abs(dst1[0]), std::abs(dst1[1])}); <<Compute MIP Map level for width and handle very wide filter>> 
int nLevels = Levels(); Float level = nLevels - 1 + Log2(std::max<Float>(width, 1e-8)); if (level >= Levels() - 1) return Texel<T>(Levels() - 1, Point2i(0, 0)); int iLevel = std::max(0, int(pstd::floor(level)));
if (options.filter == FilterFunction::Point) { <<Return point-sampled value at selected MIP level>> 
Point2i resolution = LevelResolution(iLevel); Point2i sti(pstd::round(st[0] * resolution[0] - 0.5f), pstd::round(st[1] * resolution[1] - 0.5f)); return Texel<T>(iLevel, sti);
} else if (options.filter == FilterFunction::Bilinear) { <<Return bilinear-filtered value at selected MIP level>> 
return Bilerp<T>(iLevel, st);
} else { <<Return trilinear-filtered value at selected MIP level>> 
if (iLevel == 0) return Bilerp<T>(0, st); else return Lerp(level - iLevel, Bilerp<T>(iLevel, st), Bilerp<T>(iLevel + 1, st));
}
} <<Compute EWA ellipse axes>> 
if (LengthSquared(dst0) < LengthSquared(dst1)) pstd::swap(dst0, dst1); Float longerVecLength = Length(dst0), shorterVecLength = Length(dst1);
<<Clamp ellipse vector ratio if too large>> 
if (shorterVecLength * options.maxAnisotropy < longerVecLength && shorterVecLength > 0) { Float scale = longerVecLength / (shorterVecLength * options.maxAnisotropy); dst1 *= scale; shorterVecLength *= scale; } if (shorterVecLength == 0) return Bilerp<T>(0, st);
<<Choose level of detail for EWA lookup and perform EWA filtering>> 
Float lod = std::max<Float>(0, Levels() - 1 + Log2(shorterVecLength)); int ilod = pstd::floor(lod); return Lerp(lod - ilod, EWA<T>(ilod, st, dst0, dst1), EWA<T>(ilod + 1, st, dst0, dst1));
}

The EWA filtering technique to be described shortly uses both derivatives of left-parenthesis s comma t right-parenthesis to compute an anisotropic filter—one that filters by different amounts in the different dimensions. The other three use an isotropic filter that filters both equally. The isotropic filters are more computationally efficient than the anisotropic filter, though they do not give results that are as good. For them, only a single value is needed to specify the width of the filter. The width here is conservatively chosen to avoid aliasing in both the s and t directions, though this choice means that textures viewed at an oblique angle will appear blurry, since the required sampling rate in one direction will be very different from the sampling rate along the other in this case.

<<Handle non-EWA MIP Map filter>>= 
Float width = 2 * std::max({std::abs(dst0[0]), std::abs(dst0[1]), std::abs(dst1[0]), std::abs(dst1[1])}); <<Compute MIP Map level for width and handle very wide filter>> 
int nLevels = Levels(); Float level = nLevels - 1 + Log2(std::max<Float>(width, 1e-8)); if (level >= Levels() - 1) return Texel<T>(Levels() - 1, Point2i(0, 0)); int iLevel = std::max(0, int(pstd::floor(level)));
if (options.filter == FilterFunction::Point) { <<Return point-sampled value at selected MIP level>> 
Point2i resolution = LevelResolution(iLevel); Point2i sti(pstd::round(st[0] * resolution[0] - 0.5f), pstd::round(st[1] * resolution[1] - 0.5f)); return Texel<T>(iLevel, sti);
} else if (options.filter == FilterFunction::Bilinear) { <<Return bilinear-filtered value at selected MIP level>> 
return Bilerp<T>(iLevel, st);
} else { <<Return trilinear-filtered value at selected MIP level>> 
if (iLevel == 0) return Bilerp<T>(0, st); else return Lerp(level - iLevel, Bilerp<T>(iLevel, st), Bilerp<T>(iLevel + 1, st));
}

Because filtering over many texels for wide filter widths would be inefficient, this method chooses a MIP map level from the pyramid such that the filter region at that level would cover four texels at that level. Figure 10.17 illustrates this idea.

Figure 10.17: Choosing a MIP Map Level for the Triangle Filter. The MIPMap chooses a level such that the filter covers four texels.

Since the resolutions of the levels of the pyramid are all powers of two, the resolution of level l is 2 Superscript normal n normal upper L normal e normal v normal e normal l normal s minus 1 minus l . Therefore, to find the level with a texel spacing width w requires solving

StartFraction 1 Over w EndFraction equals 2 Superscript normal n normal upper L normal e normal v normal e normal l normal s minus 1 minus l

for l . In general, this will be a floating-point value between two MIP map levels. Values of l greater than the number of pyramid levels correspond to a filter width wider than the image, in which case the single pixel at the top level is returned.

<<Compute MIP Map level for width and handle very wide filter>>= 
int nLevels = Levels(); Float level = nLevels - 1 + Log2(std::max<Float>(width, 1e-8)); if (level >= Levels() - 1) return Texel<T>(Levels() - 1, Point2i(0, 0)); int iLevel = std::max(0, int(pstd::floor(level)));

For a point-sampled texture lookup, it is only necessary to convert the continuous texture coordinates over left-bracket 0 comma 1 right-bracket to discrete coordinates over the image resolution and to retrieve the appropriate texel value via the MIPMap’s Texel() method.

<<Return point-sampled value at selected MIP level>>= 
Point2i resolution = LevelResolution(iLevel); Point2i sti(pstd::round(st[0] * resolution[0] - 0.5f), pstd::round(st[1] * resolution[1] - 0.5f)); return Texel<T>(iLevel, sti);

Bilinear filtering, which is equivalent to filtering using a triangle filter, is easily implemented via a call to Bilerp().

<<Return bilinear-filtered value at selected MIP level>>= 
return Bilerp<T>(iLevel, st);

Bilinear interpolation is provided in a separate method so that it can also be used for trilinear filtering.

<<MIPMap Method Definitions>>+=  
template <> RGB MIPMap::Bilerp(int level, Point2f st) const { if (int nc = pyramid[level].NChannels(); nc == 3 || nc == 4) return RGB(pyramid[level].BilerpChannel(st, 0, wrapMode), pyramid[level].BilerpChannel(st, 1, wrapMode), pyramid[level].BilerpChannel(st, 2, wrapMode)); else { Float v = pyramid[level].BilerpChannel(st, 0, wrapMode); return RGB(v, v, v); } }

As shown by Figure 10.17, applying a triangle filter to the four texels around the sample point will either filter over too small a region or too large a region (except for very carefully selected filter widths). Therefore, the Trilinear filtering option applies the triangle filter at both of these levels and blends between them according to how close level is to each of them. This helps hide the transitions from one MIP map level to the next at nearby pixels in the final image. While applying a triangle filter to four texels at two levels in this manner does not generally give exactly the same result as applying a triangle filter to the original pixels, the difference is not too bad in practice, and the efficiency of this approach is worth this penalty. In any case, the following elliptically weighted average filtering approach should be used when texture quality is important.

<<Return trilinear-filtered value at selected MIP level>>= 
if (iLevel == 0) return Bilerp<T>(0, st); else return Lerp(level - iLevel, Bilerp<T>(iLevel, st), Bilerp<T>(iLevel + 1, st));

The elliptically weighted average (EWA) algorithm fits an ellipse to the two differential vectors in texture space and then filters the texture with a Gaussian filter function (Figure 10.18). It is widely regarded as one of the best texture filtering algorithms in graphics and has been carefully derived from the basic principles of sampling theory. Unlike the triangle filter, it can filter over arbitrarily oriented regions of the texture, with different filter extents in different directions. The quality of its results is improved by it being an anisotropic filter, since it can adapt to different sampling rates along the two image axes.

Figure 10.18: The EWA filter applies a Gaussian filter to the texels in an elliptical area around the evaluation point. The extent of the ellipse is such that its edge passes through the positions of the adjacent texture samples as estimated by the texture coordinate partial derivatives.

We will not show the full derivation of this filter here, although we do note that it is distinguished by being a unified resampling filter: it simultaneously computes the result of a Gaussian filtered texture function convolved with a Gaussian reconstruction filter in image space. This is in contrast to many other texture filtering methods that ignore the effect of the image-space filter or equivalently assume that it is a box. Even if a Gaussian is not being used for filtering the samples for the image being rendered, taking some account of the spatial variation of the image filter improves the results, assuming that the filter being used is somewhat similar in shape to the Gaussian, as the Mitchell and windowed sinc filters are.

The screen-space partial derivatives of the texture coordinates define the ellipse. The lookup method starts out by determining which of the two axes is the longer of the two, swapping them if needed so that dst0 is the longer vector. The length of the shorter vector will be used to select a MIP map level.

<<Compute EWA ellipse axes>>= 
if (LengthSquared(dst0) < LengthSquared(dst1)) pstd::swap(dst0, dst1); Float longerVecLength = Length(dst0), shorterVecLength = Length(dst1);

Next the ratio of the length of the longer vector to the length of the shorter one is considered. A large ratio indicates a very long and skinny ellipse. Because this method filters texels from a MIP map level chosen based on the length of the shorter differential vector, a large ratio means that a large number of texels need to be filtered. To avoid this expense (and to ensure that any EWA lookup takes a bounded amount of time), the length of the shorter vector may be increased to limit this ratio. The result may be an increase in blurring, although this effect usually is not noticeable in practice.

<<Clamp ellipse vector ratio if too large>>= 
if (shorterVecLength * options.maxAnisotropy < longerVecLength && shorterVecLength > 0) { Float scale = longerVecLength / (shorterVecLength * options.maxAnisotropy); dst1 *= scale; shorterVecLength *= scale; } if (shorterVecLength == 0) return Bilerp<T>(0, st);

Like the triangle filter, the EWA filter uses the image pyramid to reduce the number of texels to be filtered for a particular texture lookup, choosing a MIP map level based on the length of the shorter vector. Given the limited ratio from the clamping above, the total number of texels used is thus bounded. Given the length of the shorter vector, the computation to find the appropriate pyramid level is the same as was used for the triangle filter. Similarly, the implementation here blends between the filtered results at the two levels around the computed level of detail, again to reduce artifacts from transitions from one level to another.

<<Choose level of detail for EWA lookup and perform EWA filtering>>= 
Float lod = std::max<Float>(0, Levels() - 1 + Log2(shorterVecLength)); int ilod = pstd::floor(lod); return Lerp(lod - ilod, EWA<T>(ilod, st, dst0, dst1), EWA<T>(ilod + 1, st, dst0, dst1));

The MIPMap::EWA() method actually applies the filter at a particular level.

<<MIPMap Method Definitions>>+= 
template <typename T> T MIPMap::EWA(int level, Point2f st, Vector2f dst0, Vector2f dst1) const { if (level >= Levels()) return Texel<T>(Levels() - 1, {0, 0}); <<Convert EWA coordinates to appropriate scale for level>> 
Point2i levelRes = LevelResolution(level); st[0] = st[0] * levelRes[0] - 0.5f; st[1] = st[1] * levelRes[1] - 0.5f; dst0[0] *= levelRes[0]; dst0[1] *= levelRes[1]; dst1[0] *= levelRes[0]; dst1[1] *= levelRes[1];
<<Find ellipse coefficients that bound EWA filter region>> 
Float A = Sqr(dst0[1]) + Sqr(dst1[1]) + 1; Float B = -2 * (dst0[0] * dst0[1] + dst1[0] * dst1[1]); Float C = Sqr(dst0[0]) + Sqr(dst1[0]) + 1; Float invF = 1 / (A * C - Sqr(B) * 0.25f); A *= invF; B *= invF; C *= invF;
<<Compute the ellipse’s left-parenthesis s comma t right-parenthesis bounding box in texture space>> 
Float det = -Sqr(B) + 4 * A * C; Float invDet = 1 / det; Float uSqrt = SafeSqrt(det * C), vSqrt = SafeSqrt(A * det); int s0 = pstd::ceil(st[0] - 2 * invDet * uSqrt); int s1 = pstd::floor(st[0] + 2 * invDet * uSqrt); int t0 = pstd::ceil(st[1] - 2 * invDet * vSqrt); int t1 = pstd::floor(st[1] + 2 * invDet * vSqrt);
<<Scan over ellipse bound and evaluate quadratic equation to filter image>> 
T sum{}; Float sumWts = 0; for (int it = t0; it <= t1; ++it) { Float tt = it - st[1]; for (int is = s0; is <= s1; ++is) { Float ss = is - st[0]; <<Compute squared radius and filter texel if it is inside the ellipse>> 
Float r2 = A * Sqr(ss) + B * ss * tt + C * Sqr(tt); if (r2 < 1) { int index = std::min<int>(r2 * MIPFilterLUTSize, MIPFilterLUTSize - 1); Float weight = MIPFilterLUT[index]; sum += weight * Texel<T>(level, {is, it}); sumWts += weight; }
} } return sum / sumWts;
}

This method first converts from texture coordinates in left-bracket 0 comma 1 right-bracket to coordinates and differentials in terms of the resolution of the chosen MIP map level. It also subtracts 0.5 from the continuous position coordinate to align the sample point with the discrete texel coordinates, as was done in MIPMap::Bilerp().

<<Convert EWA coordinates to appropriate scale for level>>= 
Point2i levelRes = LevelResolution(level); st[0] = st[0] * levelRes[0] - 0.5f; st[1] = st[1] * levelRes[1] - 0.5f; dst0[0] *= levelRes[0]; dst0[1] *= levelRes[1]; dst1[0] *= levelRes[0]; dst1[1] *= levelRes[1];

It next computes the coefficients of the implicit equation for the ellipse centered at the origin that is defined by the vectors (ds0,dt0) and (ds1,dt1). Placing the ellipse at the origin rather than at left-parenthesis s comma t right-parenthesis simplifies the implicit equation and the computation of its coefficients and can be easily corrected for when the equation is evaluated later. The general form of the implicit equation for all points left-parenthesis s comma t right-parenthesis inside such an ellipse is

e left-parenthesis s comma t right-parenthesis equals upper A s squared plus upper B s t plus upper C t squared less-than upper F comma

although it is more computationally efficient to divide through by upper F and express this as

e left-parenthesis s comma t right-parenthesis equals StartFraction upper A Over upper F EndFraction s squared plus StartFraction upper B Over upper F EndFraction s t plus StartFraction upper C Over upper F EndFraction t squared equals upper A prime s squared plus upper B prime s t plus upper C prime t squared less-than 1 period

We will not derive the equations that give the values of the coefficients, although the interested reader can easily verify their correctness.

<<Find ellipse coefficients that bound EWA filter region>>= 
Float A = Sqr(dst0[1]) + Sqr(dst1[1]) + 1; Float B = -2 * (dst0[0] * dst0[1] + dst1[0] * dst1[1]); Float C = Sqr(dst0[0]) + Sqr(dst1[0]) + 1; Float invF = 1 / (A * C - Sqr(B) * 0.25f); A *= invF; B *= invF; C *= invF;

The next step is to find the axis-aligned bounding box in discrete integer texel coordinates of the texels that are potentially inside the ellipse. The EWA algorithm loops over all of these candidate texels, filtering the contributions of those that are in fact inside the ellipse. The bounding box is found by determining the minimum and maximum values that the ellipse takes in the s and t directions. These extrema can be calculated by finding the partial derivatives partial-differential e slash partial-differential s and partial-differential e slash partial-differential t , finding their solutions for s equals 0 and t equals 0 , and adding the offset to the ellipse center. For brevity, we will not include the derivation for these expressions here.

<<Compute the ellipse’s left-parenthesis s comma t right-parenthesis bounding box in texture space>>= 
Float det = -Sqr(B) + 4 * A * C; Float invDet = 1 / det; Float uSqrt = SafeSqrt(det * C), vSqrt = SafeSqrt(A * det); int s0 = pstd::ceil(st[0] - 2 * invDet * uSqrt); int s1 = pstd::floor(st[0] + 2 * invDet * uSqrt); int t0 = pstd::ceil(st[1] - 2 * invDet * vSqrt); int t1 = pstd::floor(st[1] + 2 * invDet * vSqrt);

Figure 10.19: Finding the r squared Ellipse Value for the EWA Filter Table Lookup.

Now that the bounding box is known, the EWA algorithm loops over the texels, transforming each one to the coordinate system where the texture lookup point left-parenthesis s comma t right-parenthesis is at the origin with a translation. It then evaluates the ellipse equation to see if the texel is inside the ellipse (Figure 10.19) and computes the filter weight for the texel if so. The final filtered value returned is a weighted sum over texels upper T left-parenthesis s prime comma t prime right-parenthesis inside the ellipse, where f is the Gaussian filter function:

StartFraction sigma-summation f left-parenthesis s prime minus s comma t prime minus t right-parenthesis upper T left-parenthesis s prime comma t Superscript prime Baseline right-parenthesis Over sigma-summation f left-parenthesis s prime minus s comma t prime minus t right-parenthesis EndFraction period

<<Scan over ellipse bound and evaluate quadratic equation to filter image>>= 
T sum{}; Float sumWts = 0; for (int it = t0; it <= t1; ++it) { Float tt = it - st[1]; for (int is = s0; is <= s1; ++is) { Float ss = is - st[0]; <<Compute squared radius and filter texel if it is inside the ellipse>> 
Float r2 = A * Sqr(ss) + B * ss * tt + C * Sqr(tt); if (r2 < 1) { int index = std::min<int>(r2 * MIPFilterLUTSize, MIPFilterLUTSize - 1); Float weight = MIPFilterLUT[index]; sum += weight * Texel<T>(level, {is, it}); sumWts += weight; }
} } return sum / sumWts;

A nice feature of the implicit equation e left-parenthesis s comma t right-parenthesis is that its value at a particular texel is the squared ratio of the distance from the center of the ellipse to the texel to the distance from the center of the ellipse to the ellipse boundary along the line through that texel (Figure 10.19). This value can be used to index into a precomputed lookup table of Gaussian filter function values.

<<Compute squared radius and filter texel if it is inside the ellipse>>= 
Float r2 = A * Sqr(ss) + B * ss * tt + C * Sqr(tt); if (r2 < 1) { int index = std::min<int>(r2 * MIPFilterLUTSize, MIPFilterLUTSize - 1); Float weight = MIPFilterLUT[index]; sum += weight * Texel<T>(level, {is, it}); sumWts += weight; }

The lookup table is precomputed and available as a constant array. Similar to the GaussianFilter used for image reconstruction, the filter function is offset so that it goes to zero at the end of its extent rather than having an abrupt step. It is

normal e Superscript minus alpha r squared Baseline minus normal e Superscript negative alpha Baseline period

alpha equals 2 was used for the table in pbrt. Because the table is indexed with squared distances from the filter center r squared , each entry stores a value normal e Superscript minus alpha r , rather than normal e Superscript minus alpha r squared .

<<MIPMap EWA Lookup Table Definition>>= 
static constexpr int MIPFilterLUTSize = 128; static PBRT_CONST Float MIPFilterLUT[MIPFilterLUTSize] = { <<MIPMap EWA Lookup Table Values>> 
0.864664733f, 0.849040031f, 0.83365953f, 0.818519294f, 0.80361563f, 0.788944781f, 0.774503231f, 0.760287285f, 0.746293485f, 0.732518315f, 0.718958378f, 0.705610275f, 0.692470789f, 0.679536581f, 0.666804492f, 0.654271305f, 0.641933978f, 0.629789352f, 0.617834508f, 0.606066525f, 0.594482362f, 0.583079159f, 0.571854174f, 0.560804546f, 0.549927592f, 0.539220572f, 0.528680861f, 0.518305838f, 0.50809288f, 0.498039544f, 0.488143265f, 0.478401601f, 0.468812168f, 0.45937258f, 0.450080454f, 0.440933526f, 0.431929469f, 0.423066139f, 0.414341331f, 0.405752778f, 0.397298455f, 0.388976216f, 0.380784035f, 0.372719884f, 0.364781618f, 0.356967449f, 0.34927541f, 0.341703475f, 0.334249914f, 0.32691282f, 0.319690347f, 0.312580705f, 0.305582166f, 0.298692942f, 0.291911423f, 0.285235822f, 0.278664529f, 0.272195935f, 0.265828371f, 0.259560347f, 0.253390193f, 0.247316495f, 0.241337672f, 0.235452279f, 0.229658857f, 0.223955944f, 0.21834214f, 0.212816045f, 0.207376286f, 0.202021524f, 0.196750447f, 0.191561714f, 0.186454013f, 0.181426153f, 0.176476851f, 0.171604887f, 0.166809067f, 0.162088141f, 0.157441005f, 0.152866468f, 0.148363426f, 0.143930718f, 0.139567271f, 0.135272011f, 0.131043866f, 0.126881793f, 0.122784719f, 0.11875169f, 0.114781633f, 0.11087364f, 0.107026696f, 0.103239879f, 0.0995122194f, 0.0958427936f, 0.0922307223f, 0.0886750817f, 0.0851749927f, 0.0817295909f, 0.0783380121f, 0.0749994367f, 0.0717130303f, 0.0684779733f, 0.0652934611f, 0.0621587038f, 0.0590728968f, 0.0560353249f, 0.0530452281f, 0.0501018465f, 0.0472044498f, 0.0443523228f, 0.0415447652f, 0.0387810767f, 0.0360605568f, 0.0333825648f, 0.0307464004f, 0.0281514227f, 0.0255970061f, 0.0230824798f, 0.0206072628f, 0.0181707144f, 0.0157722086f, 0.013411209f, 0.0110870898f, 0.0087992847f, 0.0065472275f, 0.00433036685f, 0.0021481365f, 0.f
};