## 6.6 Bilinear Patches

It is useful to have a shape defined by four vertices. One option would be a planar quadrilateral, though not requiring all four vertices to be coplanar is preferable, as it is less restrictive. Such a shape is the bilinear patch, which is a parametric surface defined by four vertices , , , and . Each vertex gives the position associated with a corner of the parametric domain and points on the surface are defined via bilinear interpolation:

The bilinear patch is a *doubly ruled surface*: there are two straight
lines through every point on it. (This can be seen by considering a
parametric point on the surface and then fixing either of and
and considering the function that results: it is linear.)

Not only can bilinear patches be used to represent planar quadrilaterals, but they can also represent simple curved surfaces. They are a useful target for converting higher-order parametric surfaces to simpler shapes that are amenable to direct ray intersection. Figure 6.22 shows two bilinear patches.

`pbrt` allows the specification of bilinear patch meshes for the same
reasons that triangle meshes can be specified: to allow per-vertex
attributes like position and surface normal to be shared by multiple
patches and to allow mesh-wide properties to be stored just once.
To this end, `BilinearPatchMesh` plays the equivalent role to the
`TriangleMesh`.

We will skip past the `BilinearPatchMesh` constructor, as it mirrors
the `TriangleMesh`’s, transforming the positions and normals to
rendering space and using the `BufferCache` to avoid storing redundant
buffers in memory.

The `BilinearPatch` class implements the `Shape` interface and
represents a single patch in a bilinear patch mesh.

`dpdu`and

`dpdv`to updated partial derivatives>>

`pError`>>

`SurfaceInteraction`for bilinear patch intersection>>

`p00`,

`p01`,

`p10`, and

`p11`>>

`area`>>

`p00`,

`p01`,

`p10`, and

`p11`>>

Also similar to triangles, each `BilinearPatch` stores the index of
the mesh that it is a part of as well as its own index in the mesh’s
patches.

The `GetMesh()` method makes it easy for a `BilinearPatch` to get
the pointer to its associated mesh.

There is a subtlety that comes with the use of a `vector` to store
the meshes. `pbrt`’s scene initialization code in Appendix C
does its best to parallelize its work, which includes the parallelization
of reading binary files that encode meshes from disk. A mutex is used to
protect adding meshes to this vector, though as this vector grows, it is
periodically reallocated to make more space. A consequence is that the
`BilinearPatch` constructor must not call the `GetMesh()` method to
get its `BilinearPatchMesh *`, since `GetMesh()` accesses `allMeshes`
without mutual exclusion. Thus, the mesh is passed to the constructor as a
parameter above.

The area of a parametric surface defined over is given by the integral

The partial derivatives of a bilinear patch are easily derived. They are:

However, it is not generally possible to
evaluate the area integral from Equation (6.12)
in closed form with these partial derivatives. Therefore, the
`BilinearPatch` constructor caches the
patch’s surface area in a member variable, using numerical integration to
compute its value if necessary.

Because bilinear patches are often used to represent rectangles, the constructor checks for that case and takes the product of the lengths of the sides of the rectangle to compute the area when appropriate. In the general case, the fragment <<Compute approximate area of bilinear patch>> uses a Riemann sum evaluated at points to approximate Equation (6.12). We do not include that code fragment here.

`area`>>=

`p00`,

`p01`,

`p10`, and

`p11`>>

This fragment, which loads the four vertices of a patch into local variables, will be reused in many of the following methods.

`p00`,

`p01`,

`p10`, and

`p11`>>=

In addition to the surface area computation, there will be a number of
additional cases where we will find it useful to use specialized algorithms
if a `BilinearPatch` is a rectangle. Therefore, this check is
encapsulated in the `IsRectangle()` method.

It first tests to see if any two neighboring vertices are coincident, in which case the patch is certainly not a rectangle. This check is important to perform first, since the following ones would otherwise end up trying to perform invalid operations like normalizing degenerate vectors in that case.

`p00`,

`p01`,

`p10`, and

`p11`>>

If the four vertices are not coplanar, then they do not form a rectangle. We can check this case by computing the surface normal of the plane formed by three of the vertices and then testing if the vector from one of those three to the fourth vertex is not (nearly) perpendicular to the plane normal.

Four coplanar vertices form a rectangle if they all have the same distance from the average of their positions. The implementation here computes the squared distance to save the square root operations and then tests the relative error with respect to the first squared distance. Because the test is based on relative error, it is not sensitive to the absolute size of the patch; scaling all the vertex positions does not affect it.

With the area cached, implementation of the `Area()` method is trivial.

The bounds of a bilinear patch are given by the bounding box that bounds
its four corner vertices. As with `Triangle`s, the mesh vertices are
already in rendering space, so no further transformation is necessary here.

`p00`,

`p01`,

`p10`, and

`p11`>>

Although a planar patch has a single surface normal, the surface normal of a nonplanar patch varies across its surface.

`p00`,

`p01`,

`p10`, and

`p11`>>

`n00`at >>

`n10`,

`n01`, and

`n11`>>

If the bilinear patch is actually a triangle, the <<If patch is a
triangle, return bounds for single surface normal>> fragment evaluates its
surface normal and returns the corresponding `DirectionCone`. We have
not included that straightforward fragment here.

Otherwise, the normals are computed at the four corners of the patch. The following fragment computes the normal at the parametric position. It is particularly easy to evaluate the partial derivatives at the corners; they work out to be the differences with the adjacent vertices in and . Some care is necessary with the orientation of the normals, however. As with triangle meshes, if per-vertex shading normals were specified, they determine which side of the surface the geometric normal lies on. Otherwise, the normal may need to be flipped, depending on the user-specified orientation and the handedness of the rendering-to-object-space transformation.

`n00`at >>=

Normals at the other three vertices are computed in an equivalent manner, so the fragment that handles the rest is not included here.

A bounding cone for the normals is found by taking their average and then finding the cosine of the maximum angle that any of them makes with their average. Although this does not necessarily give an optimal bound, it usually works well in practice. (See the “Further Reading” section in Chapter 3 for more information on this topic.)

### 6.6.1 Intersection Tests

Unlike triangles (but like spheres and cylinders), a ray may intersect a bilinear patch twice, in which case the closest of the two intersections is returned. An example is shown in Figure 6.23.

As with triangles, it is useful to have a stand-alone ray–bilinear patch
intersection test function rather than only providing this functionality
through an instance of a `BilinearPatch` object. Rather than being based on
computing values along the ray and then finding the coordinates
for any found intersections, the algorithm here first determines the
parametric coordinates of any intersections. Only if any are found
within are the corresponding and values computed to find
the full intersection information.

`eps`to ensure that candidate is greater than zero>>

`u`,

`v`, and

`t`if intersection is valid>>

`tMax`and possibly return intersection>>

Going back to the definition of the bilinear surface, Equation (6.11), we can see that if we fix one of or , then we are left with an equation that defines a line. For example, with fixed, we have

with

(See Figure 6.24.)

The first step of the intersection test considers the set of all such lines defined by the patch’s vertices. For any intersection, the minimum distance from a point along the ray to a point along one of these lines will be zero. Therefore, we start with the task of trying to find values that give lines with zero distance to a point on the ray.

Given two infinite and non-parallel lines, one defined by the two points and and the other defined by and , the minimum distance between them can be found by determining the pair of parallel planes that each contain one of the lines and then finding the distance between them. (See Figure 6.25.)

To find the coefficients of those plane equations, we start by taking the cross product . This gives a vector that is perpendicular to both lines and provides the first three coefficients of the plane equation . In turn, the and coefficients can be found for each line’s plane by substituting a point along the respective line into the plane equation and solving for . Because the planes are parallel, the distance between them is then

In the case of ray–bilinear patch intersection, one line corresponds to the ray and the other to a line from the family of lines given by .

Given a ray and the bilinear patch vertices, we have , the ray’s origin, and can be set as the point along the ray . Then, and can be set as and from Equation (6.14). After taking the cross product to find the plane coefficients, finding each value, and simplifying, we can find that is a quadratic equation in . (That it is quadratic is reassuring, since a ray can intersect a bilinear patch up to two times.)

Because we only care about finding zeros of the distance function, we can neglect the denominator of Equation (6.15). After equating the difference to 0, collecting terms and simplifying, we end up with the following code to compute the quadratic coefficients.

The values where the ray intersects the patch are given by the solution to the corresponding quadratic equation. If there are no real solutions, then there is no intersection and the function returns.

The two values are handled in turn. The first step is to check whether each is between 0 and 1. If not, it does not represent a valid intersection in the patch’s parametric domain. Otherwise, the and values for the intersection point are computed.

`u`,

`v`, and

`t`if intersection is valid>>

One way to compute the and values is to find the parametric values
along the ray and the line where the distance between them is
minimized. Although this distance should be zero since we have determined
that there is an intersection between the ray and , there may be some
round-off error in the computed `u` value. Thus, formulating this
computation in terms of minimizing that distance is a reasonable
way to make the most of the values at hand.

With the ray’s origin and its direction, the parameter values where the distances are minimized are given by

and

where is shorthand for the determinant of the matrix formed from the three column vectors. We will not derive these equations here. The “Further Reading” section has more details.

We start by computing a handful of common values that are used in computing the matrix determinants and final parametric values.

The matrix determinants in the numerators can easily be computed using the
`SquareMatrix` class. Note that there are some common subexpressions
among the two of them, though we leave it to the compiler to handle them.
In a more optimized implementation, writing out the determinant
computations explicitly in order to do so manually could be worthwhile.

Due to round-off error, it is possible that the computed distance is
positive and seemingly represents a valid intersection even though the true
value of is negative and corresponds to a point behind the ray’s
origin. Testing against an epsilon value (which is discussed further
in Section 6.8.7) helps avoid reporting
incorrect intersections in such cases. Because we defer the division to
compute the final value, it is necessary to test `t1` against
`p2 * eps` here.

`u`,

`v`, and

`t`if intersection is valid>>=

The second root is handled with equivalent code, though with added logic to keep the closer of the intersections if there are two of them. That fragment is not included here.

If the final closest value is less than the given `tMax`, then an
intersection is returned.

`tMax`and possibly return intersection>>=

The coordinates and ray parametric value are sufficient to encapsulate the intersection so that the rest of its geometric properties can be computed later.

The `InteractionFromIntersection()` method computes all the
geometric information necessary to return the `SurfaceInteraction`
corresponding to a specified point on a bilinear patch, as is found
by the intersection routine.

`dpdu`and

`dpdv`to updated partial derivatives>>

`pError`>>

`SurfaceInteraction`for bilinear patch intersection>>

Given the parametric coordinates of an intersection point, it is easy to compute the corresponding point on the bilinear patch using Equation (6.11) and its partial derivatives with Equation (6.13).

`p00`,

`p01`,

`p10`, and

`p11`>>

If per-vertex texture coordinates have been specified, then they, too, are interpolated at the intersection point. Otherwise, the parametric coordinates are used for texture mapping. For the remainder of this method, we will denote the texture coordinates as to distinguish them from the patch’s parameterization. (Because this method does not use the parametric distance along the ray, this notation is unambiguous.)

Variables are also defined here to store the partial derivatives between the two sets of coordinates: , , , and . These are initialized for now to the appropriate values for when .

`dpdu`and

`dpdv`to updated partial derivatives>>

If per-vertex texture coordinates have been specified, they are bilinearly interpolated in the usual manner.

Because the partial derivatives and in the
`SurfaceInteraction` are in terms of the parameterization used
for texturing, these values must be updated if texture coordinates have
been specified.

`dpdu`and

`dpdv`to updated partial derivatives>>

The first step is to compute the updated partial derivatives and so forth. These can be found by first taking the corresponding partial derivatives of the bilinear interpolation used to compute to find and . (Note the similar form to how the partial derivatives of were computed earlier.) The desired partial derivatives can be found by taking reciprocals.

Given the partial derivatives, the chain rule can be applied to compute the updated partial derivatives of position. For example,

and similarly for .

If the provided texture coordinates specify a degenerate mapping, or may be zero. In that case,
`dpdu` and `dpdv` are left unchanged, as at least their cross
product provides a correct normal vector. A dot product checks that the
normal given by
lies in the same hemisphere as the normal given by the cross product of the
original partial derivatives of , flipping if necessary. Finally, `dpdu` and `dpdv` can be updated.

`dpdu`and

`dpdv`to updated partial derivatives>>=

The second partial derivatives of are easily found to compute the partial derivatives of the surface normal; all but are zero vectors. Thence, the partial derivatives of the normal can be computed using the regular approach. These are then adjusted to account for the parameterization in the same way that and were. The corresponding fragment follows the same form as <<Compute partial derivatives of with respect to >> and is therefore not included here.

All the necessary information for initializing the
`SurfaceInteraction` is now at hand.

`SurfaceInteraction`for bilinear patch intersection>>=

Shading geometry is set in the `SurfaceInteraction` after it is
created. Therefore, per-vertex shading normals are handled next.

The usual bilinear interpolation is performed and if the resulting normal
is non-degenerate, the shading geometry is provided to the
`SurfaceInteraction`.

The partial derivatives of the shading normal are computed in the same
manner as the partial derivatives of were found, including the
adjustment for the parameterization given by per-vertex texture
coordinates, if provided. Because shading geometry is specified via
shading and vectors, here we find the rotation matrix that
takes the geometric normal to the shading normal and apply it to
`dpdu` and `dpdv`. The cross product of the resulting vectors
then gives the shading normal.

Given the intersection and `InteractionFromIntersection()` methods,
both of the
`BilinearPatch::Intersect()`
and `IntersectP()` methods are
easy to implement. Since they both follow what should be by now a familiar
form, we have not included them here.

### 6.6.2 Sampling

The sampling routines for bilinear patches select between sampling algorithms depending on the characteristics of the patch. For area sampling, both rectangular patches and patches that have an emission distribution defined by an image map are given special treatment. When sampling by solid angle from a reference point, rectangular patches are projected on to the sphere and sampled as spherical rectangles. For both cases, general-purpose sampling routines are used otherwise.

The area sampling method first samples parametric coordinates, from which the rest of the necessary geometric values are derived.

`p00`,

`p01`,

`p10`, and

`p11`>>

`w`array with differential area at bilinear patch corners>>

`pError`for sampled bilinear patch >>

`ShapeSample`for sampled bilinear patch point>>

While all the `Shape` implementations we have implemented so far
can be used as area light sources, none of their sampling routines have
accounted for the fact that `pbrt`’s `DiffuseAreaLight` allows
specifying an image that is used to represent spatially varying emission
over the shape’s surface. Because such emission profiles are most
frequently used with rectangular light sources, the `BilinearPatch` has
the capability of sampling in according to the emission function.
Figure 6.26 demonstrates the value of
doing so.

*(Bunny model courtesy of the Stanford Computer Graphics Laboratory.)*

Otherwise, if the patch is not a rectangle, an approximation to uniform
area sampling is used. If it is a rectangle, then uniform area sampling is
trivial and the provided sample value is used directly for . In all
of these cases, the `pdf` value is with respect to the
parametric domain over .

`w`array with differential area at bilinear patch corners>>

For patches without an emissive image to sample from, we would like to uniformly sample over their surface area, as we have done with all the shapes so far. Unfortunately, an exact equal-area sampling algorithm cannot be derived for bilinear patches. This is a consequence of the fact that it is not possible to integrate the expression that gives the area of a bilinear patch, Equation (6.12). It is easy to uniformly sample in parametric space, but doing so can give a poor distribution of samples, especially for bilinear patches where two vertices are close together and the others are far apart. Figure 6.27 shows such a patch as well as the distribution of points that results from uniform parametric sampling. While the nonuniform distribution of points can be accounted for in the PDF such that Monte Carlo estimates using such samples still converge to the correct result, the resulting estimators will generally have higher error than if a more uniform distribution is used.

An exact equal-area sampling algorithm would sample points with probability proportional to its differential surface area . Lacking the ability to sample directly from this distribution, we will approximate it with a bilinear function where the value of the function at each corner is given by the patch’s differential surface area there. Sampling a location from that distribution generally works well to approximate exact equal-area sampling; see Figure 6.28.

`w`array with differential area at bilinear patch corners>>

It is especially easy to compute the partial derivatives at the patch corners; they are just differences with the adjacent vertices.

`w`array with differential area at bilinear patch corners>>=

Given a position on the patch, the corresponding position, partial
derivatives, and surface normal can all be computed, following the same
approach as was implemented in `InteractionFromIntersection()`. The
fragment <<Compute , , and for sampled
>> is therefore not included here, and texture coordinates for
the sampled point are computed using a fragment defined earlier.

`pError`for sampled bilinear patch >>

Only the geometric normal is needed for sampled points. It is easily found
via the cross product of partial derivatives of the position. The
<<Flip normal at sampled if necessary>> fragment negates the
normal if necessary, depending on the mesh properties and shading normals,
if present. It follows the same form as earlier fragments that orient the
geometric normal based on the shading normal, if present, and otherwise the
`reverseOrientation` and `transformSwapsHandedness` properties of
the mesh.

The PDF value stored in `pdf` gives the probability of
sampling the position `uv` in parametric space. In order to return a
PDF defined with respect to the patch’s surface area, it is necessary to
account for the corresponding change of variables, which results in an
additional factor of .

`ShapeSample`for sampled bilinear patch point>>=

The PDF for sampling a given point on a bilinear patch is found by first computing the probability density for sampling its position in parametric space and then transforming that density to be with respect to the patch’s surface area.

`p00`,

`p01`,

`p10`, and

`p11`>>

`intr.uv`>>

`w`array with differential area at bilinear patch corners>>

If coordinates have been specified at the vertices of a bilinear
patch, then the member variable `Interaction::uv` stores interpolated
texture coordinates. In the following, we will need the parametric
coordinates over the patch’s domain, which can be found via a
call to `InvertBilinear()`.

Regardless of which sampling technique is used for a bilinear patch, finding the PDF for a sample is straightforward.

`intr.uv`>>=

`w`array with differential area at bilinear patch corners>>

The partial derivatives are computed from as they have been before
and so we omit the corresponding fragment. Given `dpdu` and
`dpdv`, the same scaling factor as was used in `Sample()` is used
to transform the PDF.

The solid angle sampling method handles general bilinear patches with the
same approach that was used for `Cylinder`s and `Disk`s: a sample
is taken with respect to surface area on the surface using the first
sampling method and is returned with a probability density expressed in
terms of solid angle. Rectangular patches are handled using a specialized
sampling technique that gives better results.

`p00`,

`p01`,

`p10`, and

`p11`>>

`wi`>>

`ss`to solid angle measure>>

`u`to account for incident factor>>

If the patch is not a rectangle, is very small, or has a sampling distribution associated with it to match textured emission, then the regular area sampling method is used and the PDF with respect to area is converted to be with respect to solid angle. Otherwise, a specialized solid angle sampling technique is used.

`wi`>>

`ss`to solid angle measure>>

`u`to account for incident factor>>

Rectangular patches are projected on the sphere to form a spherical rectangle that can then be sampled directly. Doing so gives similar benefits to sampling spherical triangles, as was implemented in Section 6.5.4.

`u`to account for incident factor>>

Also as with triangles, we incorporate an approximation to the
factor at the reference point by warping the uniform sample
`u` using a bilinear approximation to the function in
the sampling space.

`u`to account for incident factor>>=

The spherical rectangle sampling algorithm maintains the relationship between each corner of the sampling space and the corresponding corner of the patch in its parametric space. Therefore, each bilinear sampling weight is set using the factor to the corresponding vertex of the patch.

In addition to the sample `u`, `SampleSphericalRectangle()`
takes a reference point, the corner of the rectangle, and
two vectors that define its edges. It returns a point on the patch and the
sampling PDF, which is one over the solid angle that it subtends.

The implementation of the `SampleSphericalRectangle()` function is
another interesting exercise in spherical trigonometry, like
`SampleSphericalTriangle()` was. However, in the interests of space,
we will not include discussion of its implementation here; the “Further
Reading” section has a pointer to the paper that introduced the approach
and describes its derivation.

A rectangle has the same surface normal across its entire surface, which
can be found by taking the cross product of the two edge vectors. The
parametric coordinates for the point `p` are then found by
computing the normalized projection of `p` onto each of the edges.

If the bilinear patch has per-vertex texture coordinates associated with it, the interpolated texture coordinates at the sampled point are easily computed.

The associated `PDF()` method follows the usual form, determining
which sampling method would be used for the patch and then computing its
solid angle PDF.

`p00`,

`p01`,

`p10`, and

`p11`>>

`ctx`>>

In all cases, the `SurfaceInteraction` corresponding to the
intersection of the ray from `ctx` in the direction `wi` will be
needed, so the method starts by performing a ray–patch intersection test.

`ctx`>>=

If one of the area sampling approaches was used, then a call to the other
`PDF()` method provides the PDF with respect to surface area, which is
converted to be with respect to solid angle before it is returned.

Otherwise, the spherical rectangle sampling technique was used. The
uniform solid angle PDF is 1 over the solid angle that the rectangle
subtends. If the reference point is on a surface and the approximation to
the factor would be applied in the `Sample()` method,
then the PDF for sampling `u` must be included as well.

The `InvertSphericalRectangleSample()` function, not included here,
returns the sample value `u` that maps to the given point `pRect`
on the rectangle when the `SampleSphericalRectangle()` is called with
the given reference point `pRef`.