1.7.1 Class Design

A key part of any graphics program is to have good classes or routines for geometric entities such as vectors and matrices, as well as graphics entities such as RGB colors and images. These routines should be made as clean and efficient as possible. A universal design question is whether locations and displacements should be separate classes because they have different operations; e.g., a location multiplied by one-half makes no geometric sense while one-half of a displacement does (Goldman, 1985; DeRose, 1989). There is little agreement on this question, which can spur hours of heated debate among graphics practitioners, but for the sake of example, let’s assume we will not make the distinction.

This implies that some basic classes to be written include

• vector2. A 2D vector class that stores an x- and y-component. It should store these components in a length-2 array so that an indexing operator can be well supported. You should also include operations for vector addition, vector subtraction, dot product, cross product, scalar multiplication, and scalar division.

• vector3. A 3D vector class analogous to vector2.

• hvector. A homogeneous vector with four components (see Chapter 8).

• rgb. An RGB color that stores three components. You should also include operations for RGB addition, RGB subtraction, RGB multiplication, scalar multiplication, and scalar division.

• transform. A 4 × 4 matrix for transformations. You should include a matrix multiply and member functions to apply to locations, directions, and surface normal vectors. As shown in Chapter 7, these are all different.

• image. A 2D array of RGB pixels with an output operation.

In addition, you might or might not want to add classes for intervals, orthonormal bases, and coordinate frames.

1.7.2 Float vs. Double

Modern architecture suggests that keeping memory use down and maintaining coherent memory access are the keys to efficiency. This suggests using single-precision data. However, avoiding numerical problems suggests using double-precision arithmetic. The tradeoffs depend on the program, but it is nice to have a default in your class definitions.

1.7.3 Debugging Graphics Programs

If you ask around, you may find that as programmers become more experienced, they use traditional debuggers less and less. One reason for this is that using such debuggers is more awkward for complex programs than for simple programs. Another reason is that the most difficult errors are conceptual ones where the wrong thing is being implemented, and it is easy to waste large amounts of time stepping through variable values without detecting such cases. We have found several debugging strategies to be particularly useful in graphics.

if x = 126 and y = 247 then
   print “blarg!”

2.1.1 Inverse Mappings

If we have a function f : A ⟼ B, there may exist an inverse function f^–1: B ⟼ A, which is defined by the rule f^–1(b) = a where b = f (a) . This definition only works if every b ∈ B is an image of some point under f (i.e., the range equals the target) and if there is only one such point (i.e., there is only one a for which f (a) = b). Such mappings or functions are called bijections. A bijection maps every a ∈ A to a unique b ∈ B, and for every b ∈ B, there is exactly one a ∈ A such that f (a) = b (Figure 2.1). A bijection between a group of riders and horses indicates that everybody rides a single horse, and every horse is ridden. The two functions would be rider (horse) and horse (rider). These are inverse functions of each other. Functions that are not bijections have no inverse (Figure 2.2).

An example of a bijection is f : ℝ ⟼ ℝ, with f (x) = x³. The inverse is . This example shows that the standard notation can be somewhat awkward because x is used as a dummy variable in both f and f^–1. It is sometimes more intuitive to use different dummy variables, with y = f (x) and x = f^–1(y) . This yields the more intuitive y = x³ and . An example of a function that does not have an inverse is sqr : ℝ ⟼ ℝ, where sqr(x) = x². This is true for two reasons: first x² = (–x)², and second no members of the domain map to the negative portions of the target. Note that we can define an inverse if we restrict the domain and range to R⁺. Then, is a valid inverse.

2.1.2 Intervals

Often, we would like to specify that a function deals with real numbers that are restricted in value. One such constraint is to specify an interval. An example of an interval is the real numbers between zero and one, not including zero or one. We denote this (0, 1) . Because it does not include its endpoints, this is referred to as an open interval. The corresponding closed interval, which does contain its endpoints, is denoted with square brackets: [0, 1]. This notation can be mixed; i.e., [0, 1) includes zero but not one. When writing an interval [a, b], we assume that a ≤ b. The three common ways to represent an interval are shown in Figure 2.3. The Cartesian products of intervals are often used. For example, to indicate that a point x is in the unit cube in 3D, we say x ∈ [0, 1]³.

Intervals are particularly useful in conjunction with set operations: intersection, union, and difference. For example, the intersection of two intervals is the set of points they have in common. The symbol ∩ is used for intersection. For example, [3, 5)∩[4, 6] = [4, 5) . For unions, the symbol ∪ is used to denote points in either interval. For example, [3, 5) ∪ [4, 6] = [3, 6]. Unlike the first two operators, the difference operator produces different results depending on argument order. The minus sign is used for the difference operator, which returns the points in the left interval that are not also in the right. For example, [3, 5) – [4, 6] = [3, 4) and [4, 6] – [3, 5) = [5, 6]. These operations are particularly easy to visualize using interval diagrams (Figure 2.4).

2.1.3 Logarithms

Although not as prevalent today as they were before calculators, logarithms are often useful in problems where equations with exponential terms arise. By definition, every logarithm has a base a. The “log base a”of x is written_a x and is defined as “the exponent to which a must be raised to get x,” i.e.,

Note that the logarithm base a and the function that raises a to a power are inverses of each other. This basic definition has several consequences:

When we apply calculus to logarithms, the special number e = 2.718... often turns up. The logarithm with base e is called the natural logarithm. We adopt the common shorthand ln to denote it:

Note that the “≡” symbol can be read “is equivalent by definition.” Like π, the special number e arises in a remarkable number of contexts. Many fields use a particular base in addition to e for manipulations and omit the base in their notation, i.e., log x. For example, astronomers often use base 10 and theoretical computer scientists often use base 2. Because computer graphics borrows technology from many fields, we will avoid this shorthand.

The derivatives of logarithms and exponents illuminate why the natural logarithm is “natural”:

The constant multipliers above are unity only for a = e.

2.3.1 Angles

Although we take angles somewhat for granted, we should return to their definition so we can extend the idea of the angle onto the sphere. An angle is formed between two half-lines (infinite rays stemming from an origin) or directions, and some convention must be used to decide between the two possibilities for the angle created between them as shown in Figure 2.6. An angle is defined by the length of the arc segment it cuts out on the unit circle. A common convention is that the smaller arc length is used, and the sign of the angle is determined by the order in which the two half-lines are specified. Using that convention, all angles are in the range [–π, π].

Each of these angles is the length of the arc of the unit circle that is “cut” by the two directions. Because the perimeter of the unit circle is 2π, the two possible angles sum to 2π. The unit of these arc lengths is radians. Another common unit is degrees, where the perimeter of the circle is 360°. Thus, an angle that is π radians is 180°, usually denoted 180°. The conversion between degrees and radians is

2.3.2 Trigonometric Functions

Given a right triangle with sides of length a, o, and h, where h is the length of the longest side (which is always opposite the right angle), or hypotenuse, an important relation is described by the Pythagorean theorem:

You can see that this is true from Figure 2.7, where the big square has area (a+o)², the four triangles have the combined area 2ao, and the center square has area h².

Because the triangles and inner square subdivide the larger square evenly, we have 2ao + h² = (a + o)², which is easily manipulated to the form above.

We define sine and cosine of ϕ, as well as the other ratio-based trigonometric expressions:

These definitions allow us to set up polar coordinates, where a point is coded as a distance from the origin and a signed angle relative to the positive x-axis (Figure 2.8). Note the convention that angles are in the range ϕ ∈ (–π, π], and that the positive angles are counterclockwise from the positive x-axis. This convention that counterclockwise maps to positive numbers is arbitrary, but it is used in many contexts in graphics so it is worth committing to memory.

Trigonometric functions are periodic and can take any angle as an argument. For example, sin(A) = sin(A +2π) . This means the functions are not invertible when considered with the domain R. This problem is avoided by restricting the range of standard inverse functions, and this is done in a standard way in almost all modern math libraries (e.g., Plauger (1991)). The domains and ranges are

The last function, atan2(s, c) is often very useful. It takes an s value proportional to sin A and a c value that scales cos A by the same factor and returns A. The factor is assumed to be positive. One way to think of this is that it returns the angle of a 2D Cartesian point (s, c) in polar coordinates (Figure 2.9).

2.3.3 Useful Identities

This section lists without derivation a variety of useful trigonometric identities.

Shifting identities:

Pythagorean identities:

Half-angle identities:

Half-angle identities:

Product identities:

The following identities are for arbitrary triangles with side lengths a, b, and c, each with an angle opposite it given by A, B, C, respectively (Figure 2.10),

The area of a triangle can also be computed in terms of these side lengths:

2.3.4 Solid Angles and Spherical Trigonometry

Traditional trigonometry in this section deals with triangles on the plane. Triangles can be defined on non-planar surfaces as well, and one that arises in many fields, astronomy, for example, is triangles on the unit-radius sphere. These spherical triangles have sides that are segments of the great circles (unit-radius circles) on the sphere. The study of these triangles is a field called spherical trigonometry and is not used that commonly in graphics, but sometimes, it is critical when it does arise. We wont discuss the details of it here, but want the reader to be aware that area exists for when those problems do arise, and there are a lot of useful rules such as a spherical law of cosines and a spherical law of sines. For an example of the machinery of spherical trigonometry being used, see the paper on sampling triangle lights (which project to a spherical triangle) (Arvo, 1995b).

Of more central importance to computer graphics are solid angles. While angles allow us to quantify things like “what is the separation of those two poles in my visual field,” solid angles let us quantify things like “how much of my visual field does that airplane cover.” For traditional angles, we project the posts onto the unit circle and measure arc length between them on the unit circle. We work with angles often enough that many of us can forget this definition because it is all so intuitive to us now. Solid angles are just as simple, but they may seem more confusing because most of us learn about them as adults. For solid angles, we project the visible directions that “see” the airplane and project it onto the unit sphere and measure the area. This area is the solid angle in the same way the arc length is the angle. While angles are measured in radians and sum to 2π (the total length of a unit circle), solid angles are measured in steradians and sum to 4π (the total area of a unit sphere).

2.4.1 Vector Operations

Vectors have most of the usual arithmetic operations that we associate with real numbers. Two vectors are equal if and only if they have the same length and direction. Two vectors are added according to the parallelogram rule. This rule states that the sum of two vectors is found by placing the tail of either vector against the head of the other (Figure 2.12). The sum vector is the vector that “completes the triangle” started by the two vectors. The parallelogram is formed by taking the sum in either order. This emphasizes that vector addition is commutative:

Note that the parallelogram rule just formalizes our intuition about displacements. Think of walking along one vector, tail to head, and then walking along the other. The net displacement is just the parallelogram diagonal. You can also create a unary minus for a vector: –a (Figure 2.13) is a vector with the same length as a but opposite direction. This allows us to also define subtraction:

You can visualize vector subtraction with a parallelogram (Figure 2.14). We can write

Vectors can also be multiplied. In fact, there are several kinds of products involving vectors. First, we can scale the vector by multiplying it by a real number k.

This just multiplies the vector’s length without changing its direction. For example, 3.5a is a vector in the same direction as a, but it is 3.5 times as long as a. We discuss two products involving two vectors, the dot product and the cross product, later in this section, and a product involving three vectors, the determinant, in Chapter 6.

2.4.2 Cartesian Coordinates of a Vector

A 2D vector can be written as a combination of any two nonzero vectors which are not parallel. This property of the two vectors is called linear independence. Two linearly independent vectors form a 2D basis, and the vectors are thus referred to as basis vectors. For example, a vector c may be expressed as a combination of two basis vectors a and b (Figure 2.15):

Note that the weights a_c and b_c are unique. Bases are especially useful if the two vectors are orthogonal; i.e., they are at right angles to each other. It is even more useful if they are also unit vectors in which case they are orthonormal. If we assume two such “special” vectors x and y are known to us, then we can use them to represent all other vectors in a Cartesian coordinate system, where each vector is represented as two real numbers. For example, a vector a might be represented as

where x_a and y_a are the real Cartesian coordinates of the 2D vector a (Figure 2.16). Note that this is not really any different conceptually from Equation (2.3), where the basis vectors were not orthonormal. But there are several advantages to a Cartesian coordinate system. For instance, by the Pythagorean theorem, the length of a is

It is also simple to compute dot products, cross products, and coordinates of vectors in Cartesian systems, as we’ll see in the following sections.

By convention, we write the coordinates of a either as an ordered pair (x_a,y_a) or a column matrix:

The form we use will depend on typographic convenience. We will also occasionally write the vector as a row matrix, which we will indicate as a^T:

We can also represent 3D, 4D, etc., vectors in Cartesian coordinates. For the 3D case, we use a basis vector z that is orthogonal to both x and y.

2.4.3 Dot Product

The simplest way to multiply two vectors is the dot product. The dot product of a and b is denoted a · b and is often called the scalar product because it returns a scalar. The dot product returns a value related to its arguments’ lengths and the angle ϕ between them (Figure 2.17):

The most common use of the dot product in graphics programs is to compute the cosine of the angle between two vectors.

The dot product can also be used to find the projection of one vector onto another. This is the length a→b of a vector a that is projected at right angles onto a vector b (Figure 2.18):

The dot product obeys the familiar associative and distributive properties we have in real arithmetic:

If 2D vectors a and b are expressed in Cartesian coordinates, we can take advantage of x · x = y · y = 1 and x · y = 0 to derive that their dot product is

Similarly in 3Dwe can find

2.4.4 Cross Product

The cross product a × b is usually used only for three-dimensional vectors; generalized cross products are discussed in references given in the chapter notes. The cross product returns a 3D vector that is perpendicular to the two arguments of the cross product. The length of the resulting vector is related to sin ϕ:

The magnitude ||a × b|| is equal to the area of the parallelogram formed by vectors a and b. In addition, a × b is perpendicular to both a and b (Figure 2.19). Note that there are only two possible directions for such a vector. By definition, the vectors in the direction of the x-, y- and z-axes are given by

and we set as a convention that x × y must be in the plus or minus z direction. The choice is somewhat arbitrary, but it is standard to assume that

All possible permutations of the three Cartesian unit vectors are

Because of the sin ϕ property, we also know that a vector cross itself is the zero vector, so x × x = 0 and so on. Note that the cross product is not commutative, i.e., x × y = y × x. The careful observer will note that the above discussion does not allow us to draw an unambiguous picture of how the Cartesian axes relate. More specifically, if we put x and y on a sidewalk, with x pointing east and y pointing north, then does z point up to the sky or into the ground? The usual convention is to have z point to the sky. This is known as a right-handed coordinate system. This name comes from the memory scheme of “grabbing” x with your right palm and fingers and rotating it toward y. The vector z should align with your thumb. This is illustrated in Figure 2.20.

The cross product has the nice property that

and

However, a consequence of the right-hand rule is

In Cartesian coordinates, we can use an explicit expansion to compute the cross product:

So, in coordinate form,

2.4.5 Orthonormal Bases and Coordinate Frames

Managing coordinate systems is one of the core tasks of almost any graphics program; the key to this is managing orthonormal bases. Any set of two 2D vectors u and v form an orthonormal basis provided that they are orthogonal (at right angles) and are each of unit length. Thus,

and

In 3D, three vectors u, v,and w form an orthonormal basis if

and

This orthonormal basis is right-handed provided

and otherwise, it is left-handed.

Note that the Cartesian canonical orthonormal basis is just one of infinitely many possible orthonormal bases. What makes it special is that it and its implicit origin location are used for low-level representation within a program. Thus, the vectors x, y, and z are never explicitly stored and neither is the canonical origin location o. The global model is typically stored in this canonical coordinate system, and it is thus often called the global coordinate system. However, if we want to use another coordinate system with origin p and orthonormal basis vectors u, v, and w, then we do store those vectors explicitly. Such a system is called a frame of reference or coordinate frame. For example, in a flight simulator, we might want to maintain a coordinate system with the origin at the nose of the plane, and the orthonormal basis aligned with the airplane. Simultaneously, we would have the master canonical coordinate system (Figure 2.21). The coordinate system associated with a particular object, such as the plane, is usually called a local coordinate system.

At a low level, the local frame is stored in canonical coordinates. For example, if u has coordinates (x_u,y_u,z_u) ,

A location implicitly includes an offset from the canonical origin:

where (x_p,y_p,z_p) are the coordinates of p.

Note that if we store a vector a with respect to the u-v-w frame, we store a triple (u_a,v_a,w_a) which we can interpret geometrically as

To get the canonical coordinates of a vector a stored in the u-v-w coordinate system, simply recall that u, v,and w are themselves stored in terms of Cartesian coordinates, so the expression u_au + v_av + w_aw is already in Cartesian coordinates if evaluated explicitly. To get the u-v-w coordinates of a vector b stored in the canonical coordinate system, we can use dot products:

This works because we know that for some u_b, v_b,and w_b,

and the dot product isolates the u_b coordinate:

This works because u, v,and w are orthonormal.

Using matrices to manage changes of coordinate systems is discussed in Sections 7.2.1 and 7.5.

2.4.6 Constructing a Basis from a Single Vector

Often we need an orthonormal basis that is aligned with a given vector. That is, given a vector a, we want an orthonormal u, v, and w such that w points in the same direction as a (Hughes & Möller, 1999), but we don’t particularly care what u and v are. One vector isn’t enough to uniquely determine the answer; we just need a robust procedure that will find any one of the possible bases.

This can be done using cross products as follows. First, make w a unit vector in the direction of a:

Then, choose any vector t not collinear with w, and use the cross product to build a unit vector u perpendicular to w:

If t is collinear with w, the denominator will vanish, and if they are nearly collinear, the results will have low precision. A simple procedure to find a vector sufficiently different from w is to start with t equal to w and change the smallest magnitude component of t to 1. For example, if then . Once w and u are in hand, completing the basis is simple:

An example of a situation where this construction is used is surface shading, where a basis aligned to the surface normal is needed but the rotation around the normal is often unimportant.

For serious production code, recently researchers at Pixar have developed a rather remarkable method for constructing a vector from two vectors that is impressive in its compactness and efficiency (Duff et al., 2017). They provide battle-tested code, and readers are encouraged to use it as there are not “gotchas” that have emerged as it used throughout the industry.

2.4.7 Constructing a Basis from Two Vectors

The procedure in the previous section can also be used in situations where the rotation of the basis around the given vector is important. A common example is building a basis for a camera: it’s important to have one vector aligned in the direction the camera is looking, but the orientation of the camera around that vector is not arbitrary, and it needs to be specified somehow. Once the orientation is pinned down, the basis is completely determined.

A common way to fully specify a frame is by providing two vectors a (which specifies w) and b (which specifies v). If the two vectors are known to be perpendicular, it is a simple matter to construct the third vector by u = b × a.

To be sure that the resulting basis really is orthonormal, even if the input vectors weren’t quite, a procedure much like the single-vector procedure is advisable:

In fact, this procedure works just fine when a and b are not perpendicular. In this case, w will be constructed exactly in the direction of a,and v is chosen to be the closest vector to b among all vectors perpendicular to w.

This procedure won’t work if a and b are collinear. In this case, b is of no help in choosing which of the directions perpendicular to a we should use: it is perpendicular to all of them.

In the example of specifying camera positions (Section 4.3), we want to construct a frame that has w parallel to the direction the camera is looking, and v should point out the top of the camera. To orient the camera upright, we build the basis around the view direction, using the straight-up direction as the reference vector to establish the camera’s orientation around the view direction. Setting v as close as possible to straight up exactly matches the intuitive notion of “holding the camera straight.”

2.4.8 Squaring Up a Basis

Occasionally, you may find problems caused in your computations by a basis that is supposed to be orthonormal but where error has crept in—due to rounding error in computation, or to the basis having been stored in a file with low precision, for instance.

The procedure of the previous section can be used; simply constructing the basis anew using the existing w and v vectors will produce a new basis that is orthonormal and is close to the old one.

This approach is good for many applications, but it is not the best available. It does produce accurately orthogonal vectors, and for nearly orthogonal starting bases, the result will not stray far from the starting point. However, it is asymmetric: it “favors” w over v and v over u (whose starting value is thrown away). It chooses a basis close to the starting basis but has no guarantee of choosing the closest orthonormal basis. When this is not good enough, the SVD (Section 6.4.1) can be used to compute an orthonormal basis that is guaranteed to be closest to the original basis.

2.5.1 Averages and weighted averages

Integrals compute the total of things. Lengths, areas, volumes, etc. But they are often used to compute averages. For example, we can compute the total volume of a region by integrating the elevation over a region (like a country).

float volume = integrate(elevation(), country)

But we could also compute the average elevation:

float averageElevation = integrate(elevation(),
      country) / integrate(1, country)

This is basically “divide the volume by the area.” This can be abstracted as

Float averageElevation = average(elevation, country)

We can also take a weighted average. Here, we add a weighting function to emphasize some points in the average more than others. For example, if we want to emphasize a parts of the region by the temperature (this is pretty arbitrary, and we will see more graphics relevant examples in the next section):

float weightedAverageElevation =
      integrate(temperature()*elevation(),
      country) / integrate(temperature(), country)

It’s a good idea to keep an eye out for this form; often integrals contain a weighted average without explicitly pointing that out and it can sometimes help intuition.

2.5.2 Integrals over solid angle

One example of a type of integral we see a lot is one of these forms or something related:

float shade = integrate(cos()*f*(),
   unit-hemisphere)

Note that since integrate(cos(), unit-hemisphere) = pi, the weighted average version is just

float shade = integrate((1/pi)*cos()*f*(),
       unit-hemisphere)

The more traditional form of this integral is

Or with spherical coordinates as we might use to solve such integrals algebraically:

The sine term if an area-correction factor for spherical coordinates. Note that in graphics, we will rarely need to write that all out and will use simpler forms without explicit coordinates as we numerically solve the integrals

The particular integral above is the shade of a perfectly reflective matte (dif-fuse) surface, and it is also a weighted average of all incident colors. This structure can be great for intuition; the color of a surface is usually related to a weighted average of incident colors.

The integrals over solid angle are almost always the same but use a wide variety of notations. Key is to recognize this is just notations and map the notations you see to one you are most comfortable with. This is much like reading pseudocode!

2.7.1 2D Implicit Curves

Intuitively, a curve is a set of points that can be drawn on a piece of paper without lifting the pen. A common way to describe a curve is using an implicit equation. An implicit equation in two dimensions has the form

The function f (x, y) returns a real value. Points (x, y) where this value is zero are on the curve, and points where the value is nonzero are not on the curve. For example, let’s say that f (x, y) is

where (x_c,y_c) is a 2D point and r is a nonzero real number. If we take f (x, y) = 0, the points where this equality holds are on the circle with center (x_c,y_c) and radius r. The reason that this is called an “implicit” equation is that the points (x, y) on the curve cannot be immediately calculated from the equation and instead must be determined by solving the equation. Thus, the points on the curve are not generated by the equation explicitly, but they are buried somewhere implicitly in the equation.

It is interesting to note that f does have values for all (x, y) . We can think of f as a terrain, with sea level at f = 0 (Figure 2.24). The shore is the implicit curve. The value of f is the altitude. Another thing to note is that the curve partitions space into regions where f > 0, f < 0, and f = 0. So you evaluate f to decide whether a point is “inside” a curve. Note that f (x, y) = c is a curve for any constant c,and c = 0 is just used as a convention. For example, if f (x, y) = x² + y² – 1, varying c just gives a variety of circles centered at the origin (Figure 2.25).

We can compress our notation using vectors. If we have c = (x_c,y_c) and p = (x, y) , then our circle with center c and radius r is defined by those position vectors that satisfy

This equation, if expanded algebraically, will yield Equation (2.9), but it is easier to see that this is an equation for a circle by “reading” the equation geometrically. It reads, “points p on the circle have the following property: the vector from c to p when dotted with itself has value r².” Because a vector dotted with itself is just its own length squared, we could also read the equation as, “points p on the circle have the following property: the vector from c to p has squared length r².”

Even better, is to observe that the squared length is just the squared distance from c to p, which suggests the equivalent form

and, of course, this suggests

The above could be read “the points p on the circle are those a distance r from the center point c,” which is as good a definition of circle as any. This illustrates that the vector form of an equation often suggests more geometry and intuition than the equivalent full-blown Cartesian form with x and y. For this reason, it is usually advisable to use vector forms when possible. In addition, you can support a vector class in your code; the code is cleaner when vector forms are used. The vector-oriented equations are also less error prone in implementation: once you implement and debug vector types in your code, the cut-and-paste errors involving x, y,and z will go away. It takes a little while to get used to vectors in these equations, but once you get the hang of it, the payoff is large.

2.7.2 The 2D Gradient

If we think of the function f (x, y) as a height field with height = f (x, y) , the gradient vector points in the direction of maximum upslope, i.e., straight uphill. The gradient vector ∇f(x, y) is given by

The gradient vector evaluated at a point on the implicit curve f (x, y) = 0 is perpendicular to the tangent vector of the curve at that point. This perpendicular vector is usually called the normal vector to the curve. In addition, since the gradient points uphill, it indicates the direction of the f (x, y) > 0 region.

In the context of height fields, the geometric meaning of partial derivatives and gradients is more visible than usual. Suppose that near the point (a, b) , f (x, y) is a plane (Figure 2.26). There is a specific uphill and downhill direction. At right angles to this direction is a direction that is level with respect to the plane. Any intersection between the plane and the f (x, y) = 0 plane will be in the direction that is level. Thus, the uphill/downhill directions will be perpendicular to the line of intersection f (x, y) = 0. To see why the partial derivative has something to do with this, we need to visualize its geometric meaning. Recall that the conventional derivative of a 1D function y = g(x) is

This measures the slope of the tangent line to g (Figure 2.27).

The partial derivative is a generalization of the 1D derivative. For a 2D function f (x, y) , we can’t take the same limit for x as in Equation (2.10), because f can change in many ways for a given change in x. However, if we hold y constant, we can define an analog of the derivative, called the partial derivative (Figure 2.28):

Why is it that the partial derivatives with respect to x and y are the components of the gradient vector? Again, there is more obvious insight in the geometry than in the algebra. In Figure 2.29, we see the vector a travels along a path where f does not change. Note that this is again at a small enough scale that the surface height (x, y) = f (x, y) can be considered locally planar. From the figure, we see that the vector a = (Δx, Δy) .

Because the uphill direction is perpendicular to a, we know the dot product is equal to zero:

We also know that the change in f in the direction (x_a,y_a) equals zero:

Given any vectors (x, y) and (x,y) that are perpendicular, we know that the angle between them is 90 degrees, and thus, their dot product equals zero (recall that the dot product is proportional to the cosine of the angle between the two vectors). Thus, we have xx + yy = 0. Given (x, y) , it is easy to construct valid vectors whose dot product with (x, y) equals zero, the two most obvious being (y,–x) and (–y, x) ; you can verify that these vectors give the desired zero dot product with (x, y). A generalization of this observation is that (x, y) is perpendicular to k(y,–x) where k is any nonzero constant. This implies that

Combining Equations (2.11) and (2.12) gives

where k is any nonzero constant. By definition, “uphill” implies a positive change in f , so we would like k > 0,and k = 1 is a perfectly good convention.

As an example of the gradient, consider the implicit circle x² + y² – 1 = 0 with gradient vector (2x, 2y) , indicating that the outside of the circle is the positive region for the function f (x, y) = x² + y² – 1. Note that the length of the gradient vector can be different depending on the multiplier in the implicit equation. For example, the unit circle can be described by Ax² + Ay² – A = 0 for any nonzero A. The gradient for this curve is (2Ax, 2Ay) . This will be normal (perpendicular) to the circle, but will have a length determined by A. For A > 0, the normal will point outward from the circle, and for A < 0, it will point inward. This switch from outward to inward is as it should be, since the positive region switches inside the circle. In terms of the height-field view, h = Ax² + Ay² – A, and the circle is at zero altitude. For A > 0, the circle encloses a depression, and for A < 0, the circle encloses a bump. As A becomes more negative, the bump increases in height, but the h = 0 circle doesn’t change. The direction of maximum uphill doesn’t change, but the slope increases. The length of the gradient reflects this change in degree of the slope. So intuitively, you can think of the gradient’s direction as pointing uphill and its magnitude as measuring how uphill the slope is.

Implicit 2D Lines

The familiar “slope-intercept” form of the line is

This can be converted easily to implicit form (Figure 2.30):

Here, m is the “slope” (ratio of rise to run), and b is the y value where the line crosses the y-axis, usually called the y-intercept. The line also partitions the 2D plane, but here “inside” and “outside” might be more intuitively called “over” and “under.”

Because we can multiply an implicit equation by any constant without changing the points where it is zero, kf (x, y) = 0 is the same curve for any nonzero k. This allows several implicit forms for the same line, for example,

One reason the slope-intercept form is sometimes awkward is that it can’t represent some lines such as x = 0 because m would have to be infinite. For this reason, a more general form is often useful:

for real numbers A, B, C.

Suppose we know two points on the line, (x₀, y₀) and (x₁,y₁) . What A, B, and C describe the line through these two points? Because these points lie on the line, they must both satisfy Equation (2.15):

Unfortunately, we have two equations and three unknowns: A, B, and C. This problem arises because of the arbitrary multiplier we can have with an implicit equation. We could set C = 1 for convenience:

but we have a similar problem to the infinite slope case in slope-intercept form: lines through the origin would need to have A(0) + B(0) + 1 = 0, which is a contradiction. For example, the equation for a 45–° line through the origin can be written x – y = 0, or equally well y – x = 0, or even 17y – 17x = 0, but it cannot be written in the form Ax + By +1 = 0.

Whenever we have such pesky algebraic problems, we try to solve the problems using geometric intuition as a guide. One tool we have, as discussed in Section 2.7.2, is the gradient. For the line Ax + By + C = 0, the gradient vector is (A, B) . This vector is perpendicular to the line (Figure 2.31), and points to the side of the line where Ax + By + C is positive. Given two points on the line (x₀,y₀) and (x₁,y₁) , we know that the vector between them points in the same direction as the line. This vector is just (x₁ – x₀,y₁ – y₀), and because it is parallel to the line, it must also be perpendicular to the gradient vector (A, B) . Recall that there are an infinite number of (A, B, C) that describe the line because of the arbitrary scaling property of implicits. We want any one of the valid (A, B, C) .

We can start with any (A, B) perpendicular to (x₁–x₀,y₁–y₀). Such a vector is just (A, B) = (y₀ –y₁, x₁ – x₀) by the same reasoning as in Section 2.7.2. This means that the equation of the line through (x₀,y₀) and (x₁,y₁) is

Now we just need to find C. Because (x₀,y₀) and (x₁,y₁) are on the line, they must satisfy Equation (2.16). We can plug either value in and solve for C. Doing this for (x₀,y₀) yields C = x₀y₁ – x₁y₀, and thus, the full equation for the line is

Again, this is one of infinitely many valid implicit equations for the line through two points, but this form has no division operation and thus no numerically degenerate cases for points with finite Cartesian coordinates. A nice thing about Equation (2.17) is that we can always convert to the slope-intercept form (when it exists) by moving the non-y terms to the right-hand side of the equation and dividing by the multiplier of the y term:

An interesting property of the implicit line equation is that it can be used to find the signed distance from a point to the line. The value of Ax + By + C is proportional to the distance from the line (Figure 2.32). As shown in Figure 2.33, the distance from a point to the line is the length of the vector k(A, B) ,which is

For the point (x, y)+ k(A, B) ,thevalueof f (x, y) = Ax + By + C is

The simplification in that equation is a result of the fact that we know (x, y) is on the line, so Ax + By + C = 0. From Equations (2.18) and (2.19), we can see that the signed distance from line Ax + By + C = 0 to a point (a, b) is

Here, “signed distance” means that its magnitude (absolute value) is the geometric distance, but on one side of the line, distances are positive and on the other, they are negative. You can choose between the equally valid representations f (x, y) = 0 and –f (x, y) = 0 if your problem has some reason to prefer a particular side being positive. Note that if (A, B) is a unit vector, then f (a, b) is the signed distance. We can multiply Equation (2.17) by a constant that ensures that (A, B) is a unit vector:

Note that evaluating f (x, y) in Equation (2.20) directly gives the signed distance, but it does require a square root to set up the equation. Implicit lines will turn out to be very useful for triangle rasterization (Section 9.1.2). Other forms for 2D lines are discussed in Chapter 13.

Implicit Quadric Curves

In the previous section, we saw that a linear function f (x, y) gives rise to an implicit line f (x, y) = 0. If f is instead a quadratic function of x and y, with the general form

the resulting implicit curve is called a quadric. Two-dimensional quadric curves include ellipses and hyperbolas, as well as the special cases of parabolas, circles, and lines.

Examples of quadric curves include the circle with center (x_c, y_c) and radius r,

and axis-aligned ellipses of the form

where (x_c,y_c) is the center of the ellipse, and a and b are the minor and major semi-axes (Figure 2.34).

2.7.3 3D Implicit Surfaces

Just as implicit equations can be used to define curves in 2D, they can be used to define surfaces in 3D. As in 2D, implicit equations implicitly define a set of points that are on the surface:

Any point (x, y, z) that is on the surface results in zero when given as an argument to f . Any point not on the surface results in some number other than zero. You can check whether a point is on the surface by evaluating f , or you can check which side of the surface the point lies on by looking at the sign of f , but you cannot always explicitly construct points on the surface. Using vector notation, we will write such functions of p = (x, y, z) as

2.7.4 Surface Normal to an Implicit Surface

A surface normal (which is needed for lighting computations, among other things) is a vector perpendicular to the surface. Each point on the surface may have a different normal vector. In the same way that the gradient provides a normal to an implicit curve in 2D, the surface normal at a point p on an implicit surface is given by the gradient of the implicit function

The reasoning is the same as for the 2D case: the gradient points in the direction of fastest increase in f , which is perpendicular to all directions tangent to the surface, in which f remains constant. The gradient vector points toward the side of the surface where f (p) > 0, which we may think of as “into” the surface or “out from” the surface in a given context. If the particular form of f creates inward-facing gradients, and outward-facing gradients are desired, the surface –f (p) = 0 is the same as surface f (p) = 0 but has directionally reversed gradients, i.e., –∇f (p) = ∇(–f (p)) .

2.7.5 Implicit Planes

As an example, consider the infinite plane through point a with surface normal n. The implicit equation to describe this plane is given by

Note that a and n are known quantities. The point p is any unknown point that satisfies the equation. In geometric terms this equation says “the vector from a to p is perpendicular to the plane normal.” If p were not in the plane, then (p – a) would not make a right angle with n (Figure 2.35).

Sometimes, we want the implicit equation for a plane through points a, b, and c. The normal to this plane can be found by taking the cross product of any two vectors in the plane. One such cross product is

This allows us to write the implicit plane equation:

A geometric way to read this equation is that the volume of the parallelepiped defined by p – a, b – a, and c – a is zero; i.e., they are coplanar. This can only be true if p is in the same plane as a, b, and c. The full-blown Cartesian representation for this is given by the determinant (this is discussed in more detail in Section 6.3):

The determinant can be expanded (see Section 6.3 for the mechanics of expanding determinants) to the bloated form with many terms.

Equations (2.22) and (2.23) are equivalent, and comparing them is instructive. Equation (2.22) is easy to interpret geometrically and will yield efficient code. In addition, it is relatively easy to avoid a typographic error that compiles into incorrect code if it takes advantage of debugged cross and dot product code. Equation (2.23) is also easy to interpret geometrically and will be efficient provided an efficient 3 × 3 determinant function is implemented. It is also easy to implement without a typo if a function determinant (a, b, c) is available. It will be especially easy for others to read your code if you rename the determinant function volume. So both Equations (2.22) and (2.23) map well into code. The full expansion of either equation into x-, y-, and z-components is likely to generate typos. Such typos are likely to compile and, thus, to be especially pesky. This is an excellent example of clean math generating clean code and bloated math generating bloated code.

2.7.6 2D Parametric Curves

A parametric curve is controlled by a single parameter that can be considered a sort of index that moves continuously along the curve. Such curves have the form

Here, (x, y) is a point on the curve, and t is the parameter that influences the curve. For a given t, there will be some point determined by the functions g and h. For continuous g and h, a small change in t will yield a small change in x and y. Thus, as t continuously changes, points are swept out in a continuous curve. This is a nice feature because we can use the parameter t to explicitly construct points on the curve. Often, we can write a parametric curve in vector form,

where f is a vector-valued function, . Such vector functions can generate very clean code, so they should be used when possible.

We can think of the curve with a position as a function of time. The curve can go anywhere and could loop and cross itself. We can also think of the curve as having a velocity at any point. For example, the point p(t) is traveling slowly near t = –2 and quickly between t = 2 and t = 3. This type of “moving point” vocabulary is often used when discussing parametric curves even when the curve is not describing a moving point.

2D Parametric Lines

A parametric line in 2D that passes through points p₀ = (x₀, y₀) and p₁ = (x₁,y₁) can be written as

Because the formulas for x and y have such similar structure, we can use the vector form for p = (x, y) (Figure 2.36):

You can read this in geometric form as “start at point p₀ and go some distance toward p₁ determined by the parameter t.” A nice feature of this form is that p(0) = p₀ and p(1) = p₁. Since the point changes linearly with t, the value of t between p₀ and p₁ measures the fractional distance between the points. Points with t < 0 are to the “far” side of p₀, and points with t > 1 are to the “far” side of p₁.

Parametric lines can also be described as just a point o and a vector d:

When the vector d has unit length, the line is arc-length parameterized. This means t is an exact measure of distance along the line. Any parametric curve can be arc-length parameterized, which is obviously a very convenient form, but not all can be converted analytically.

2D Parametric Circles

A circle with center (x_c,y_c) and radius r has a parametric form:

To ensure that there is a unique parameter ϕ for every point on the curve, we can restrict its domain: ϕ ∈ [0, 2π) or ϕ ∈ (–π, π] or any other half-open interval of length 2π.

An axis-aligned ellipse can be constructed by scaling the x and y parametric equations separately:

2.7.7 3D Parametric Curves

A 3D parametric curve operates much like a 2D parametric curve:

For example, a spiral around the z-axis is written as

As with 2D curves, the functions f , g, and h are defined on a domain D ⊂ R if we want to control where the curve starts and ends. In vector form, we can write

In this chapter, we only discuss 3D parametric lines in detail. General 3D parametric curves are discussed more extensively in Chapter 15.

2.7.8 3D Parametric Surfaces

The parametric approach can be used to define surfaces in 3D space in much the same way we define curves, except that there are two parameters to address the two-dimensional area of the surface. These surfaces have the form

or, in vector form,

With implicit surfaces, the derivative of the function f gave us the surface normal. With parametric surfaces, the derivatives of p also give information about the surface geometry.

Consider the function q(t) = p(t, v₀) . This function defines a parametric curve obtained by varying u while holding v fixed at the value v₀. This curve, called an isoparametric curve (or sometimes “isoparm” for short), lies in the surface. The derivative of q gives a vector tangent to the curve, and since the curve lies in the surface, the vector q also lies in the surface. Since it was obtained by varying one argument of p, the vector q is the partial derivative of p with respect to u, which we’ll denote p_u. A similar argument shows that the partial derivative p_v gives the tangent to the isoparametric curves for constant u, which is a second tangent vector to the surface.

The derivative of p, then, gives two tangent vectors at any point on the surface. The normal to the surface may be found by taking the cross product of these vectors: since both are tangent to the surface, their cross product, which is perpendicular to both tangents, is normal to the surface. The right-hand rule for cross products provides a way to decide which side is the front, or outside, of the surface; we will use the convention that the vector

points toward the outside of the surface.

2.7.9 Summary of Curves and Surfaces

Implicit curves in 2D or surfaces in 3D are defined by scalar-valued functions of two or three variables, f : ℝ² → ℝ or f : ℝ³ → ℝ, and the surface consists of all points where the function is zero:

Parametric curves in 2D or 3D are defined by vector-valued functions of one variable, p : D ⊂ ℝ → ℝ² or p : D ⊂ ℝ → ℝ³, and the curve is swept out as t varies over all of D:

Parametric surfaces in 3D are defined by vector-valued functions of two variables, p : D ⊂ ℝ² → ℝ³, and the surface consists of the images of all points (u, v) in the domain:

For implicit curves and surfaces, the normal vector is given by the derivative of f (the gradient), and the tangent vector (for a curve) or vectors (for a surface) can be derived from the normal by constructing a basis.

For parametric curves and surfaces, the derivative of p gives the tangent vector (for a curve) or vectors (for a surface), and the normal vector can be derived from the tangents by constructing a basis.

2.9.1 2D Triangles

If we have a 2D triangle defined by 2D points a, b, and c, we can first find its area:

The derivation of this formula can be found in Section 6.3. This area will have a positive sign if the points a, b,and c are in counterclockwise order and a negative sign, otherwise.

Often in graphics, we wish to assign a property, such as color, at each triangle vertex and smoothly interpolate the value of that property across the triangle. There are a variety of ways to do this, but the simplest is to use barycentric coordinates. One way to think of barycentric coordinates is as a nonorthogonal coordinate system as was discussed briefly in Section 2.4.2. Such a coordinate system is shown in Figure 2.38, where the coordinate origin is a and the vectors from a to b and c are the basis vectors. With that origin and those basis vectors, any point p can be written as

Note that we can reorder the terms in Equation (2.28) to get

Often people define a new variable α to improve the symmetry of the equations:

which yields the equation

with the constraint that

Barycentric coordinates seem like an abstract and unintuitive construct at first, but they turn out to be powerful and convenient. You may find it useful to think of how street addresses would work in a city where there are two sets of parallel streets, but where those sets are not at right angles. The natural system would essentially be barycentric coordinates, and you would quickly get used to them. Barycentric coordinates are defined for all points on the plane. A particularly nice feature of barycentric coordinates is that a point p is inside the triangle formed by a, b,and c if and only if

If one of the coordinates is zero and the other two are between zero and one, then you are on an edge. If two of the coordinates are zero, then the other is one, and you are at a vertex. Another nice property of barycentric coordinates is that Equation (2.29) in effect mixes the coordinates of the three vertices in a smooth way. The same mixing coefficients (α, β, γ) can be used to mix other properties, such as color, as we will see in the next chapter.

Given a point p, how do we compute its barycentric coordinates? One way is to write Equation (2.28) as a linear system with unknowns β and γ,solve,andset α = 1 – β – γ . That linear system is

Although it is straightforward to solve Equation (2.31) algebraically, it is often fruitful to compute a direct geometric solution.

One geometric property of barycentric coordinates is that they are the signed scaled distance from the lines through the triangle sides, as is shown for β in Figure 2.39. Recall from Section 2.7.2 that evaluating the equation f (x, y) for the line f (x, y) = 0 returns the scaled signed distance from (x, y) to the line. Also recall that if f (x, y) = 0 is the equation for a particular line, so is kf (x, y) = 0 for any nonzero k. Changing k scales the distance and controls which side of the line has positive signed distance, and which negative. We would like to choose k such that, for example, kf (x, y) = β. Since k is only one unknown, we can force this with one constraint, namely, that at point b, we know β = 1. So if the line f_ac(x, y) = 0 goes through both a and c, then we can compute β for a point (x, y) as follows:

and we can compute γ and α in a similar fashion. For efficiency, it is usually wise to compute only two of the barycentric coordinates directly and to compute the third using Equation (2.30).

To find this “ideal” form for the line through p₀ and p₁, we can first use the technique of Section 2.7.2 to find some valid implicit lines through the vertices. Equation (2.17) gives us

Note that f_ab(x_c,y_c) probably does not equal one, so it is probably not the ideal form we seek. By dividing through by f_ab(x_c,y_c) ,weget

The presence of the division might worry us because it introduces the possibility of divide-by-zero, but this cannot occur for triangles with areas that are not near zero. There are analogous formulas for α and β, but typically only one is needed:

Another way to compute barycentric coordinates is to compute the areas A_a, A_b, and A_c, of subtriangles as shown in Figure 2.40. Barycentric coordinates obey

the rule

where A is the area of the triangle. Note that A = A_a + A_b + A_c, so it can be computed with two additions rather than a full area formula. This rule still holds for points outside the triangle if the areas are allowed to be signed. The reason for this is shown in Figure 2.41. Note that these are signed areas and will be computed correctly as long as the same signed area computation is used for both A and the subtriangles A_a, A_b,and A_c.

2.9.2 3D Triangles

One wonderful thing about barycentric coordinates is that they extend almost transparently to 3D. If we assume the points a, b, and c are 3D, then we can still use the representation

Now, as we vary β and γ, we sweep out a plane.

The normal vector to a triangle can be found by taking the cross product of any two vectors in the plane of the triangle (Figure 2.42). It is easiest to use two of the three edges as these vectors, for example,

Note that this normal vector is not necessarily of unit length, and it obeys the right-hand rule of cross products.

The area of the triangle can be found by taking the length of the cross product:

Note that this is not a signed area, so it cannot be used directly to evaluate barycentric coordinates. However, we can observe that a triangle with a “clockwise” vertex order will have a normal vector that points in the opposite direction to the normal of a triangle in the same plane with a “counterclockwise” vertex order. Recall that

where ϕ is the angle between the vectors. If a and b are parallel, then cos ϕ = ±1, and this gives a test of whether the vectors point in the same or opposite directions. This, along with Equations (2.33)–(2.35), suggest the formulas

where n is Equation (2.34) evaluated with vertices a, b, and c; n_a is Equation (2.34) evaluated with vertices b, c,and p, and so on, i.e.,

2.12.1 Importance Sampling

When a function we want to take a random average of has a wide variation in its high and low values, it can be to our advantage to concentrate samples in some areas and then correct for the nonuniformity with weights. The probability density functions give us the right tool for that: if we know the PDF of a sample, that is a direct measure of how “oversampled” that region is. If we use nonuniform samples, then we can get thus

integrate = average_of_nonuniform_samples(f()/p(),
            domain).

A neat thing about this formula is it also works for uniform random samples. In that case, the PDF p() = 1/ integrate(1, domain) so the “size” of the domain is encoded in the PDF.

For any given Monte Carlo importance sampling problem, there is a pretty formulaic approach we follow, at least to get started:

1. Identify what is the function f () and the domain of integration (e.g., points on the unit sphere or points on a triangle).

2. Pick a method for generating random samples x_i on that domain, and make sure there is a way to evaluate the PDF p(x_i) for each sample.

3. Average the ratio f (x_i)/p(x_i) for many x_i. This is our estimate of the integral.

A neat thing is that any p() can be used and you will converge to the right answer (with the caveat that where f () is nonzero p() must be nonzero). Which p() you use merely influences how fast your estimate converges. So we usually start with a constant p() for debugging our code.

3.1.1 Displays

Current displays, including televisions and digital cinematic projectors as well as displays and projectors for computers, are nearly universally based on fixed arrays of pixels. They can be separated into emissive displays, which use pixels that directly emit controllable amounts of light, and transmissive displays, in which the pixels themselves don’t emit light but instead vary the amount of light that they allow to pass through them. Transmissive displays require a light source to illuminate them: in a direct-viewed display, this is a backlight behind the array; in a projector, it is a lamp that emits light that is projected onto the screen after passing through the array. An emissive display is its own light source.

Light-emitting diode (LED) displays are an example of the emissive type. Each pixel is composed of one or more LEDs, which are semiconductor devices (based on inorganic or organic semiconductors) that emit light with intensity depending on the electrical current passing through them (see Figure 3.1).

The pixels in a color display are divided into three independently controlled subpixels—one red, one green, and one blue—each with its own LED made using different materials so that they emit light of different colors (Figure 3.2). When the display is viewed from a distance, the eye can’t separate the individual subpixels, and the perceived color is a mixture of red, green, and blue.

Liquid crystal displays (LCDs) are an example of the transmissive type. A liquid crystal is a material whose molecular structure enables it to rotate the polarization of light that passes through it, and the degree of rotation can be adjusted by an applied voltage. An LCD pixel (Figure 3.3) has a layer of polarizing film behind it, so that it is illuminated by polarized light—let’s assume it is polarized horizontally.

A second layer of polarizing film in front of the pixel is oriented to transmit only vertically polarized light. If the applied voltage is set so that the liquid crystal layer in between does not change the polarization, all light is blocked and the pixel is in the “off” (minimum intensity) state. If the voltage is set so that the liquid crystal rotates the polarization by 90°, then all the light that entered through the back of the pixel will escape through the front, and the pixel is fully “on”—it has its maximum intensity. Intermediate voltages will partly rotate the polarization so that the front polarizer partly blocks the light, resulting in intensities between the minimum and maximum (Figure 3.4). Like color LED displays, color LCDs have red, green, and blue subpixels within each pixel, which are three independent pixels with red, green, and blue color filters over them.

Any type of display with a fixed pixel grid, including these and other technologies, has a fundamentally fixed resolution determined by the size of the grid. For displays and images, resolution simply means the dimensions of the pixel grid: if a desktop monitor has a resolution of 1920 × 1200 pixels, this means that it has 2,304,000 pixels arranged in 1920 columns and 1200 rows.

An image of a different resolution, to fill the screen, must be converted into a 1920 × 1200 image using the methods of Chapter 10.

3.1.2 Hardcopy Devices

The process of recording images permanently on paper has very different constraints from showing images transiently on a display. In printing, pigments are distributed on paper or another medium so that when light reflects from the paper it forms the desired image. Printers are raster devices like displays, but many printers can only print binary images—pigment is either deposited or not at each grid position, with no intermediate amounts possible.

An ink-jet printer (Figure 3.5) is an example of a device that forms a raster image by scanning. An ink-jet print head contains liquid ink carrying pigment, which can be sprayed in very small drops under electronic control. The head moves across the paper, and drops are emitted as it passes grid positions that should receive ink; no ink is emitted in areas intended to remain blank. After each sweep, the paper is advanced slightly, and then, the next row of the grid is laid down. Color prints are made by using several print heads, each spraying ink with a different pigment, so that each grid position can receive any combination of different colored drops. Because all drops are the same, an ink-jet printer prints binary images: at each grid point, there is a drop or no drop; there are no intermediate shades.

An ink-jet printer has no physical array of pixels; the resolution is determined by how small the drops can be made and how far the paper is advanced after each sweep. Many ink-jet printers have multiple nozzles in the print head, enabling several sweeps to be made in one pass, but it is the paper advance, not the nozzle spacing, that ultimately determines the spacing of the rows.

The thermal dye transfer process is an example of a continuous tone printing process, meaning that varying amounts of dye can be deposited at each pixel—it is not all-or-nothing like an ink-jet printer (Figure 3.6). A donor ribbon containing colored dye is pressed between the paper, or dye receiver, and a print head containing a linear array of heating elements, one for each column of pixels in the image. As the paper and ribbon move past the head, the heating elements switch on and off to heat the ribbon in areas where dye is desired, causing the dye to diffuse from the ribbon to the paper. This process is repeated for each of several dye colors. Since higher temperatures cause more dye to be transferred, the amount of each dye deposited at each grid position can be controlled, allowing a continuous range of colors to be produced. The number of heating elements in the print head establishes a fixed resolution in the direction across the page, but the resolution along the page is determined by the rate of heating and cooling compared to the speed of the paper.

Unlike displays, the resolution of printers is described in terms of the pixel density instead of the total count of pixels. So a thermal dye transfer printer that has elements spaced 300 per inch across its print head has a resolution of 300 pixels per inch (ppi) across the page. If the resolution along the page is chosen to be the same, we can simply say the printer’s resolution is 300 ppi. An ink-jet printer that places dots on a grid with 1200 grid points per inch is described as having a resolution of 1200 dots per inch (dpi). Because the ink-jet printer is a binary device, it requires a much finer grid for at least two reasons. Because edges are abrupt black/white boundaries, very high resolution is required to avoid stair-stepping, or aliasing, from appearing (see Section 9.3). When continuous-tone images are printed, the high resolution is required to simulate intermediate colors by printing varying-density dot patterns called halftones.

3.1.3 Input Devices

Raster images have to come from somewhere, and any image that wasn’t computed by some algorithm has to have been measured by some raster input device, most often a camera or scanner. Even in rendering images of 3D scenes, photographs are used constantly as texture maps (see Chapter 11). A raster input device has to make a light measurement for each pixel, and (like output devices) they are usually based on arrays of sensors.

A digital camera is an example of a 2D array input device. The image sensor in a camera is a semiconductor device with a grid of light-sensitive pixels. Two common types of arrays are known as CCDs (charge-coupled devices) and CMOS (complimentary metal–oxide–semiconductor) image sensors. The camera’s lens projects an image of the scene to be photographed onto the sensor, and then, each pixel measures the light energy falling on it, ultimately resulting in a number that goes into the output image (Figure 3.7). In much the same way as color displays use red, green, and blue subpixels, most color cameras work by using a color-filter array or mosaic to allow each pixel to see only red, green, or blue light, leaving the image processing software to fill in the missing values in a process known as demosaicking (Figure 3.8).

Other cameras use three separate arrays, or three separate layers in the array, to measure independent red, green, and blue values at each pixel, producing a usable color image without further processing. The resolution of a camera is determined by the fixed number of pixels in the array and is usually quoted using the total count of pixels: a camera with an array of 3000 columns and 2000 rows produces an image of resolution 3000 × 2000, which has 6 million pixels, and is called a 6 megapixel (MP) camera. It’s important to remember that a mosaic sensor does not measure a complete color image, so a camera that measures the same number of pixels but with independent red, green, and blue measurements records more information about the image than one with a mosaic sensor.

A flatbed scanner also measures red, green, and blue values for each of a grid of pixels, but like a thermal dye transfer printer, it uses a 1D array that sweeps across the page being scanned, making many measurements per second (Figure 3.9). The resolution across the page is fixed by the size of the array, and the resolution along the page is determined by the frequency of measurements compared to the speed at which the scan head moves. A color scanner has a 3 × n_x array, where n_x is the number of pixels across the page, with the three rows covered by red, green, and blue filters. With an appropriate delay between the times at which the three colors are measured, this allows three independent color measurements at each grid point. As with continuous-tone printers, the resolution of scanners is reported in pixels per inch (ppi).

With this concrete information about where our images come from and where they will go, we’ll now discuss images more abstractly, in the way we’ll use them in graphics algorithms.

3.2.1 Pixel Values

So far we have described the values of pixels in terms of real numbers, representing intensity (possibly separately for red, green, and blue) at a point in the image. This suggests that images should be arrays of floating-point numbers, with either one (for grayscale, or black and white, images) or three (for RGB color images) 32-bit floating-point numbers stored per pixel. This format is sometimes used, when its precision and range of values are needed, but images have a lot of pixels and memory and bandwidth for storing and transmitting images are invariably scarce. Just one ten-megapixel photograph would consume about 115 MB of RAM in this format.

Less range is required for images that are meant to be displayed directly. While the range of possible light intensities is unbounded in principle, any given device has a decidedly finite maximum intensity, so in many contexts, it is perfectly sufficient for pixels to have a bounded range, usually taken to be [0, 1] for simplicity. For instance, the possible values in an 8-bit image are 0, 1/255, 2/255, ..., 254/255, 1. Images stored with floating-point numbers, allowing a wide range of values, are often called high dynamic range (HDR) images to distinguish them from fixed-range, or low dynamic range (LDR) images that are stored with integers. See Chapter 20 for an in-depth discussion of techniques and applications for high dynamic range images.

Here are some pixel formats with typical applications:

• 1-bit grayscale—text and other images where intermediate grays are not desired (high resolution required);

• 8-bit RGB fixed-range color (24 bits total per pixel)—web and email applications, consumer photographs;

• 8- or 10-bit fixed-range RGB (24–30 bits/pixel)—digital interfaces to computer displays;

• 12- to 14-bit fixed-range RGB (36–42 bits/pixel)—raw camera images for professional photography;

• 16-bit fixed-range RGB (48 bits/pixel)—professional photography and printing; intermediate format for image processing of fixed-range images;

• 16-bit fixed-range grayscale (16 bits/pixel)—radiology and medical imaging;

• 16-bit “half-precision” floating-point RGB—HDR images; intermediate format for real-time rendering;

• 32-bit floating-point RGB—general-purpose intermediate format for software rendering and processing of HDR images.

Reducing the number of bits used to store each pixel leads to two distinctive types of artifacts, or artificially introduced flaws, in images. First, encoding images with fixed-range values produces clipping when pixels that would otherwise be brighter than the maximum value are set, or clipped, to the maximum representable value. For instance, a photograph of a sunny scene may include reflections that are much brighter than white surfaces; these will be clipped (even if they were measured by the camera) when the image is converted to a fixed range to be displayed. Second, encoding images with limited precision leads to quantization artifacts, or banding, when the need to round pixel values to the nearest representable value introduces visible jumps in intensity or color. Banding can be particularly insidious in animation and video, where the bands may not be objectionable in still images, but become very visible when they move back and forth.

3.2.2 Monitor Intensities and Gamma

All modern monitors take digital input for the “value” of a pixel and convert this to an intensity level. Real monitors have some nonzero intensity when they are off because the screen reflects some light. For our purposes, we can consider this “black” and the monitor fully on as “white.” We assume a numeric description of pixel color that ranges from zero to one. Black is zero, white is one, and a gray halfway between black and white is 0.5. Note that here “halfway” refers to the physical amount of light coming from the pixel, rather than the appearance. The human perception of intensity is nonlinear and will not be part of the present discussion; see Chapter 19 for more.

There are two key issues that must be understood to produce correct images on monitors. The first is that monitors are nonlinear with respect to input. For example, if you give a monitor 0, 0.5, and 1.0 as inputs for three pixels, the intensities displayed might be 0, 0.25, and 1.0 (off, one-quarter fully on, and fully on). As an approximate characterization of this nonlinearity, monitors are commonly characterized by a γ (“gamma”) value. This value is the degree of freedom in the formula

where a is the input pixel value between zero and one. For example, if a monitor has a gamma of 2.0, and we input a value of a = 0.5, the displayed intensity will be one-fourth the maximum possible intensity because 0.5² = 0.25. Note that a = 0 maps to zero intensity and a = 1 maps to the maximum intensity regardless of the value of γ. Describing a display’s nonlinearity using γ is only an approximation; we do not need a great deal of accuracy in estimating the γ of a device. A nice visual way to gauge the nonlinearity is to find what value of a gives an intensity halfway between black and white. This a will be

If we can find that a, we can deduce γ by taking logarithms on both sides:

We can find this a by a standard technique where we display a checkerboard pattern of black and white pixels next to a square of gray pixels with input a (Figure 3.11), then ask the user to adjust a (with a slider, for instance) until the two sides match in average brightness. When you look at this image from a distance (or without glasses if you are nearsighted), the two sides of the image will look about the same when a is producing an intensity halfway between black and white. This is because the blurred checkerboard is mixing even numbers of white and black pixels so the overall effect is a uniform color halfway between white and black.

Once we know γ, we can gamma correct our input so that a value of a = 0.5 is displayed with intensity halfway between black and white. This is done with the transformation

When this formula is plugged into Equation (3.1), we get

Another important characteristic of real displays is that they take quantized input values. So while we can manipulate intensities in the floating-point range [0, 1], the detailed input to a monitor is a fixed-size integer. The most common range for this integer is 0–255 which can be held in 8 bits of storage. This means that the possible values for a are not any number in [0, 1] but instead

This means the possible displayed intensity values are approximately

where M is the maximum intensity. In applications where the exact intensities need to be controlled, we would have to actually measure the 256 possible intensities, and these intensities might be different at different points on the screen, especially for CRTs. They might also vary with viewing angle. Fortunately, few applications require such accurate calibration.

3.4.1 Image Storage

Most RGB image formats use eight bits for each of the red, green, and blue channels. This results in approximately three megabytes of raw information for a single million-pixel image. To reduce the storage requirement, most image formats allow for some kind of compression. At a high level, such compression is either lossless or lossy. No information is discarded in lossless compression, while some information is lost unrecoverably in a lossy system. Popular image storage formats include

• jpeg. This lossy format compresses image blocks based on thresholds in the human visual system. This format works well for natural images.

• tiff. This format is most commonly used to hold binary images or losslessly compressed 8- or 16-bit RGB although many other options exist.

• ppm. This very simple lossless, uncompressed format is most often used for 8-bit RGB images although many options exist.

• png. This is a set of lossless formats with a good set of open source management tools.

Because of compression and variants, writing input/output routines for images can be involved. Fortunately, one can usually rely on library routines to read and write standard file formats. For quick-and-dirty applications, where simplicity is valued above efficiency, a simple choice is to use raw ppm files, which can often be written simply by dumping the array that stores the image in memory to a file, prepending the appropriate header.

4.3.1 Orthographic Views

For an orthographic view, all the rays will have the direction –w. Even though a parallel view doesn’t have a viewpoint per se, we can still use the origin of the camera frame to define the plane where the rays start, so that it’s possible for objects to be behind the camera.

The viewing rays should start on the plane defined by the point e and the vectors u and v; the only remaining information required is where on the plane the image is supposed to be. We’ll define the image dimensions with four numbers, for the four sides of the image: l and r are the positions of the left and right edges of the image, as measured from e along the u direction; and b and t are the positions of the bottom and top edges of the image, as measured from e along the v direction. Usually, l < 0 < r and b < 0 < t. (SeeFigure4.9a.)

In Section 3.2, we discussed pixel coordinates in an image. To fit an image with n_x × n_y pixels into a rectangle of size (r – l)×(t–b) , the pixels are spaced a distance (r – l)/n_x apart horizontally and (t – b)/n_y apart vertically, with a half-pixel space around the edge to center the pixel grid within the image rectangle. This means that the pixel at position (i, j) in the raster image has the position

where (u, v) are the coordinates of the pixel’s position on the image plane, measured with respect to the origin e and the basis {u, v}.

In an orthographic view, we can simply use the pixel’s image-plane position as the ray’s starting point, and we already know the ray’s direction is the view direction. The procedure for generating orthographic viewing rays is then

compute u and v using (4.1)

ray.o ← e + u u + v v

ray.d ←–w

It’s very simple to make an oblique parallel view: just allow the image plane normal w to be specified separately from the view direction d. The procedure is then exactly the same, but with d substituted for –w. Of course, w is still used to construct u and v.

4.3.2 Perspective Views

For a perspective view, all the rays have the same origin, at the viewpoint; it is the directions that are different for each pixel. The image plane is no longer positioned at e, but rather some distance d in front of e; this distance is the image plane distance, often loosely called the focal length, because choosing d plays the same role as choosing focal length in a real camera. The direction of each ray is defined by the viewpoint and the position of the pixel on the image plane. This situation is illustrated in Figure 4.9, and the resulting procedure is similar to the orthographic one:

compute u and v using (4.1)

ray.o ← e

ray.d ←– d w + u u + v v

As with parallel projection, oblique perspective views can be achieved by specifying the image plane normal separately from the projection direction.

4.4.1 Ray-Sphere Intersection

Given a ray p(t) = e + td and an implicit surface f (p) = 0 (see Section 2.7.3), we’d like to know where they intersect. Intersection points occur when points on the ray satisfy the implicit equation, so the values of t we seek are those that solve the equation

A sphere with center c = (x_c ,y_c ,z_c) and radius R can be represented by the implicit equation

We can write this same equation in vector form:

Any point p that satisfies this equation is on the sphere. If we plug points on the ray p(t) = e + td into this equation, we get an equation in terms of t that is satisfied by the values of t that yield points on the sphere:

Rearranging terms yields

Here, everything is known except the parameter t, so this is a classic quadratic equation in t, meaning it has the form

The solution to this equation is discussed in Section 2.2. The term under the square root sign in the quadratic solution, B² – 4AC, is called the discriminant and tells us how many real solutions there are. If the discriminant is negative, its square root is imaginary and the line and sphere do not intersect. If the discriminant is positive, there are two solutions: one solution where the ray enters the sphere and one where it leaves. If the discriminant is zero, the ray grazes the sphere, touching it at exactly one point. Plugging in the actual terms for the sphere and canceling a factor of two, we get

In an actual implementation, you should first check the value of the discriminant before computing other terms. To correctly find the closest intersection in the interval [t₀,t₁], there are three cases: if the smaller of the two solutions is in the interval, it is the first hit; otherwise, if the larger solution is in the interval, it is the first hit; otherwise, there is no hit.

As discussed in Section 2.7.4, the normal vector at point p is given by the gradient n = 2(p – c) . The unit normal is (p – c)/R.

4.4.2 Ray-Triangle Intersection

There are many algorithms for computing ray-triangle intersections. We will present the form that uses barycentric coordinates for the parametric plane containing the triangle, because it requires no long-term storage other than the vertices of the triangle (Snyder & Barr, 1987).

To intersect a ray with a parametric surface, we set up a system of equations where the Cartesian coordinates all match:

Here, we have three equations and three unknowns (t, u, and v). In the case where the surface is a parametric plane, the parametric equation is linear and can be written in vector form as discussed in Section 2.9.2. If the vertices of the triangle are a, b,and c, then the intersection will occur when

for some t, β,and γ. Solving this equation tells us both t, which locates the intersection point along the ray, and (β, γ) , which locates the intersection point relative to the triangle. The intersection p will be at e+td as shown in Figure 4.10. Again from Section 2.9.2, we know the intersection is inside the triangle if and only if β > 0, γ > 0, and β + γ < 1. Otherwise, the ray has hit the plane outside the triangle, so it misses the triangle. If there are no solutions, either the triangle is degenerate or the ray is parallel to the plane containing the triangle.

To solve for t, β, and γ in Equation (4.2), we expand it from its vector form into the three equations for the three coordinates:

This can be rewritten as a standard linear system:

The fastest classic method to solve this 3 × 3 linear system is Cramer’s rule. This gives us the solutions

where the matrix A is

and |A| denotes the determinant of A. The 3 × 3 determinants have common subterms that can be exploited for efficiency in implementation. Looking at the linear systems with dummy variables

Cramer’s rule gives us

where

We can reduce the number of operations by reusing numbers such as “ei-minus-hf.”

The algorithm for the ray-triangle intersection for which we need the linear solution can have some conditions for early termination. Thus, the function should look something like:

boolean raytri (Ray r, vector3 a, vector3 b, vector3 c,
                interval [t0,t1])
compute t 
if (t < t0) or (t > t1) then
   return false
compute γ 
if (γ < 0) or (γ > 1) then
  return false
compute β 
if (β < 0) or (β > 1 – γ) then
return false
return true

4.4.3 Ray intersection in software

In a ray tracing program, it is a good idea to use an object-oriented design that has a class called something like Surface with derived classes Triangle, Sphere, etc. Anything that a ray can intersect, including groups of surfaces or efficiency structures (Section 12.3) should be a subclass of Surface. The ray-tracing program would then have one reference to a Surface for the whole model, and new types of objects and efficiency structures can be added transparently.

The key interface of the Surface class is a method to intersect a ray (Kirk & Arvo, 1988).

class Surface
    HitRecord hit(Ray r, real t0, real t1)

Here, (t₀,t₁) is the interval on the ray where hits will be returned, and HitRecord is a class that contains all the data about the surface intersection that will be needed:

class HitRecord
    Surface s  | surface that was hit
    real t | coordinate of hit point along the ray
    Vec3 n | surface normal at the hit point
    .
    .
    .

The surface that was hit, the t value, and the surface normal are the minimum required, but other data such as texture coordinates or tangent vectors may be stored as well. Depending on the language, the hit record might not be literally returned from the function but rather passed by reference and filled in. A miss can be indicated by a hit that has t = ∞.

4.4.4 Intersecting a Group of Objects

Of course, most interesting scenes consist of more than one object, and when we intersect a ray with the scene, we must find only the closest intersection to the camera along the ray. A simple way to implement this is to think of a group of objects as itself being another type of object. To intersect a ray with a group, you simply intersect the ray with the objects in the group and return the intersection with the smallest t value. The following code tests for hits in the interval t ∈ [t₀,t₁]:

class Group, subclass of Surface
    list-of-Surface surfaces | list of all surfaces in the group
    HitRecord hit(Ray ray, real t0, real t1)
       HitRecord closest-hit(∞) | initialize to indicate miss
       for surf in surfaces do
           rec = surf.hit(ray, t0, t1)
           if rec.t < ∞ then
               closest-hit = rec
               t1 = t
       return closest-hit

Note that this code shrinks the intersection interval [t₀,t₁] so that the call to surf.hit will only hit surfaces that are closer than the closest one seen so far.

Once ray-scene intersection works, we can render an image like Figure 4.11, but nicer results depend on including more visual cues, as we describe next.

4.5.1 Light Sources

To support shading, a ray tracing program always has a list of light sources. For the Chapter 5 shading model, we need three types of lights: point lights, which emit light from a point in space, directional lights, which illuminate the scene from a single direction, and ambient lights, which provide constant illumination to fill in the shadows. In fancier systems, other types of lights are supported, such as area lights (which are basically scene geometry that emits light) or environment lights (which use an image to represent light coming from far-away sources like the sky).

Computing shading from a point or directional light source requires certain geometric information, and in a ray tracer, after a viewing ray has been determined to hit the surface, we have all we need to determine these four vectors:

• The shading point x can be computed by evaluating the viewing ray at the t value of the intersection.

• The surface normal n depends on the type of surface (sphere, triangle, etc.), and every surface needs to be able to compute its normal at the point where a ray intersects it.

• The light direction l is computed from the light source position or direction as part of shading.

• The viewing direction v is simply opposite the direction of the viewing ray (v = –d/ ||d||).

The shading from an ambient source is much simpler: there is no l since light comes from everywhere; the shading does not depend on v; and for the simple models of Chapter 5, it doesn’t even depend on x or n.

Computing shading in a scene containing several lights is simply a matter of adding up the contributions of the lights. In a basic ray tracer, you can simply loop over all the light sources, computing shading from each one, and accumulate the results into the pixel color.

4.5.2 Shading in software

A ray tracing program usually contains objects representing light sources and materials. Light sources can be instances of subclasses of a Light class, and they must include enough information to fully describe the light source. Since shading also requires parameters describing the material of the surface, another class that is useful is Material, which encapsulates everything needed to evaluate the shading model.

Different systems take different approaches to breaking up the shading calculations between lights and materials. An approach that aligns with the presentation in this chapter is to make lights responsible for the overall illumination computation and materials responsible for computing BRDF values. With this setup, the interfaces of these classes might look like:

class Light
    Color illuminate(Ray ray, HitRecord hrec)

class Material
    Color evaluate(Vec3 l, Vec3 v, Vec3 n)

Each surface would then store a reference to its material, and in this way, point light illumination might be implemented as follows:

class PointLight, subclass of Light
    Color I
    Vec3 p
    Color illuminate(Ray ray, HitRecord hrec)
       Vec3 x = ray.evaluate(hrec.t)
       real r = p – x
       Vec3 l = (p – x)/r
       Vec3 n = hrec.normal
       Color E = max(0, n · l) I/r2
       Color k = hrec.surface.material.evaluate(l, v, n)
       return kE

These computations assume the class Color carries the RGB components of a color and supports componentwise multiplication. This arrangement is also amenable to treating ambient lighting as a light source, by making the ambient coefficient a property of the material:

class AmbientLight, subclass of Light
    Color Ia
    Color illuminate(Ray ray, HitRecord hrec)
       Color ka = hrec.surface.material.ka
       return ka Ia

The complete calculation for shading a ray, including the intersection and handling several lights, can look like this:

function shade-ray(Ray ray, realt0, realt1)
HitRecord rec = scene.hit(ray,t0,t1)
if rec.t < ∞ then
   Color c = 0
   for light in scene.lights do
     c = c + light.illuminate(ray, rec)
   return c
else return background-color

This setup keeps materials and lights reasonably separate and allows you to later add new kinds of materials and lights transparently. Textures add some complexity to the architecture of a ray tracer; see Section 11.2.5.

By itself, shading makes images of 3D objects more realistic and understandable, but it doesn’t show their interactions with other objects. For instance, the spheres in Figure 4.12 appear to float above the floor they are resting on.

4.5.3 Shadows

Once you have basic shading in your ray tracer, shadows for point and directional lights can be added very easily. If we imagine ourselves at a point x on a surface being shaded, the point is in shadow if we “look” towards the light source and see an object between us and the light source. If there are no objects in between, then the light is not blocked.

This is shown in Figure 4.13, where the ray x + tl does not hit any objects and thus the point x is not in shadow. On the other hand, the point x is in shadow because the ray x + tl does hit an object. The rays that determine in or out of shadow are called shadow rays to distinguish them from viewing rays.

To get the algorithm for shading, we add an if statement to the code that adds shading from a light source to first determine whether the light is shadowed. In a naive implementation, the shadow ray will check for t ∈ [0,r], but because of numerical imprecision, this can result in an intersection with the surface on which p lies. Instead, the usual adjustment to avoid that problem is to test for t ∈ [ , r] where is some small positive constant (Figure 4.14).

A shadow test can be added to the method PointLight.illuminate shown above by tracing a shadow ray and adding a conditional:

HitRecord srec = scene.hit(Ray(x, l), ,r) if srec.t < ∞ then proceed with normal illumination calculation else return 0 | shading point is in shadow

The shadow test for directional lights is similar but uses t₁ = ∞ rather than r. Note that the illumination computation for each light requires a separate shadow ray, and there is no shadow test in computing ambient shading.

Shadows serve an important visual role in showing the relationships between nearby objects, as shown in Figure 4.15.

4.5.4 Mirror Reflection

It is straightforward to add ideal specular reflection, or mirror reflection,toaray-tracing program. The key observation is shown in Figure 4.16 where a viewer looking from direction e sees what is in direction r as seen from the surface. The vector r is the reflection of the vector –d across the surface normal n, which can be computed using the projection of d onto the direction of the surface normal:

In the real world, some energy is lost when the light reflects from the surface, and this loss can be different for different colors. For example, gold reflects yellow more efficiently than blue, so it shifts the colors of the objects it reflects. This can be implemented by adding a recursive call in shade-ray that adds one more contribution after all the lights are accounted for:

where k_m (for “mirror reflection”) is the specular RGB color. We need to make sure to pass t₀ = for the same reason as we did with shadow rays; we don’t want the reflection ray to hit the object that generates it.

The problem with the recursive call above is that it may never terminate. For example, if a ray starts inside a room, it will bounce forever. This can be fixed by adding a maximum recursion depth. The code will be more efficient if a reflection ray is generated only if k_m is not zero.

Using a constant mirror reflection coefficient k_m gives a particular look characteristic of simple ray tracers (Figure 4.17); in the real world, this coefficient varies substantially depending on the incident angle. For better models, see Chapter 14.

5.1.1 Point source illumination

A point light source is described by its position, which is a point in 3D space, and its intensity, which describes the amount of light it produces. A point source can be isotropic, meaning the intensity is the same in all directions; this is normally the default, but many systems provide “spot lights” that only send light in some directions, which can be handy for controlling light in a virtual scene in the same way that a real spot light is useful for controlling light on a stage.

For an isotropic point source, it’s easy to reason about how much light falls on a surface a certain distance away. Suppose we have a point source that emits one Watt of radiant power isotropically, and we place this source at the center of a hollow sphere with one meter radius (Figure 5.1). All the power from the light falls on the inside surface of the sphere, and it’s distributed uniformly over the whole surface area of 4π m², so the density of radiant power per unit area is 1/(4π) Watts per square meter. This density is known as irradiance and is the right quantity to describe how much light is falling on a surface for the purposes of light reflection.

In the general case of a source that has power P and a receiving sphere of radius r, we find the irradiance E is

The quantity I = P/(4π) is the intensity of the source; it is a property of the source itself that is independent of what surface it’s illuminating. The r^–2 factor, often called the inverse square term, describes how irradiance depends on the distance r between the source and the surface.

One other important consideration in computing irradiance is the angle of incidence—the angle between the surface normal and the direction the light is traveling. Consider a small surface that is illuminated by a point source that is far away compared to the size of the surface. The light that falls on the surface is all travelling approximately parallel. If we tilt the surface to an angle of 60° as shown in Figure 5.2, the surface intercepts only half the light that it did when it was facing the source. In general, when rotated by an angle θ it intercepts an amount of light (radiant power) proportional to cos θ, and since the area stays the same, the irradiance (which, remember, is radiant power per unit area) is proportional to the same factor. This rule, that the irradiance on a surface falls off as the cosine of the incident angle, is known as Lambert’s cosine law because it was described by Johann Heinrich Lambert in his 1760 book Photometria.

Putting this together with the formula, we just derived for irradiance on a surface facing exactly toward the source, we get the general formula for irradiance due to a point source,

The term cos θ/r² can be called the geometry factor for a point source; it depends on the geometric relationship between source and receiving surface, but not on the specific properties of either one.

In practice, the angle θ is not normally computed, because given a unit vector n that is normal to the surface and a unit vector l that points toward the light (Figure 5.3), the cosine factor can be computed using the dot product

which is simpler and more efficient than computations with trigonometric functions.

5.1.2 Directional illumination

A directional source is a limiting case of a very bright, far-away point source. As the source gets farther and farther away, the ratio I/r² in Equation 5.2 varies less and less over the scene, and for a directional source, we replace this with a constant, H:

This constant can be called the normal irradiance since it is equal to the irradiance when the light is positioned along the surface normal. A directional source is characterized by the direction toward the source (rather than by a position) and by the normal irradiance H (rather than by an intensity). The illumination from a directional source is uniform and does not fall off with distance in the way that point source illumination does.

5.2.1 Lambertian reflection

The very simplest kind of reflection is a surface that reflects light equally to all directions, regardless of where it came from, so that the reflected light L_r seen by the observer is simply a constant multiple of the irradiance:

A surface that behaves this way is known as an ideal diffuse surface and appears the same brightness from all directions; its color is view independent and is completely described by its reflectance, R, which is the fraction of the irradiance it reflects. The coefficient relating reflected to incident light is R/π (the reason for the factor of π will have to wait for Chapter 14):

The reflectance can be different for different colors of light, and for simple modeling of color, it suffices to just keep three different reflectances, one each for red, green, and blue, so this shading equation is carried out separately for the three color channels.

Ideal diffuse shading, often called Lambertian shading because Lambert’s cosine law is the main effect it models, provides a flat, chalky appearance by itself. Physically, it models light that bounces around inside the material so that it “for-gets” where it came from and emerges randomly in all directions. It is an effective model for paper, flat paint, dirt, tree bark, stone, and other rough materials that don’t have a distinct and smooth enough top surface to produce noticeable shiny reflections.

5.2.2 Specular reflection

Many materials have some degree of shininess to them—for example, metals, plastics, gloss or semi-gloss paints, or many leaves of plants. When you look at these materials, you see reflections that move around when you move your viewpoint; you could describe their color as being view-dependent in contrast to the view-independent color of a Lambertian surface. The view-dependent part of the reflection generally happens at the top surface of the material and is known as specular reflection.

The simplest kind of specular reflection happens at perfectly smooth surfaces like a mirror or the surface of water: light reflects in a mirrorlike way so that light coming from a point source goes in exactly one direction. This is known as ideal specular reflection and generally needs to be handled as a special case. But many surfaces are not perfectly smooth, and they exhibit a more general kind of reflection known in computer graphics as glossy reflection. There are many models for glossy reflection, and better ones are discussed in Chapter 14, but a simple and well-known model was originally proposed by Phong (1975) and later updated by Blinn (1976) and others to the form most commonly used today, known as the Modified Blinn–Phong model.

Since specular reflection is view dependent, it is a function of the view vector v that points from the shading point toward the viewer, as well as the normal vector n and light direction l. The idea is to produce reflection that is at its brightest when v and l are symmetrically positioned across the surface normal, which is when mirror reflection would occur; the reflection then decreases smoothly as the vectors move away from a mirror configuration.

We can tell how close we are to a mirror configuration using the idea of a half vector, which is the vector halfway between the viewing and illumination directions and is perpendicular to the surface exactly when l and v are in a mirror reflection configuration (Figure 5.5). If the half vector is near the surface normal, the specular component should be bright; if it is far away, it should be dim. We measure the nearness of h and n by computing their dot product (remember they are unit vectors, so n · h reaches its maximum of 1 when the vectors are equal) and then take the result to a power p > 1 to make it decrease faster:

(n · h)^p

The Phong exponent, p, controls the apparent shininess of the surface: higher values make the reflection fall off faster away from the mirror direction, leading to a shinier appearance. The half vector itself is easy to compute: since v and l are the same length, their sum is a vector that bisects the angle between them, which only needs to be normalized to produce h:

To incorporate the Blinn–Phong idea into a shading computation, we add a specular component to Lambertian shading; the Lambertian part is then the diffuse component. We simply generalize the factor k from (5.3) that relates reflected light to incident irradiance to include not just the contribution of diffuse reflection but also a separate term that adds in specular reflection:

where the scale factor k_s is the specular coefficient (separate for red, green, and blue) and controls how bright the specular component is, and we have added a clamping operation to avoid surprises for corner cases in which n faces away from h.

The expression that generalizes the factor k is called the bidirectional reflectance distribution function or BRDF, because it describes how the reflectance varies as a function of the two directions l and v. The BRDF of a Lambertian surface is constant, but the BRDF of a surface that has specular reflection is not. The shading calculation then boils down to computing the irradiance (describing how much light is available to reflect) and the BRDF (describing how the surface reflects it), and then multiplying them. The BRDF is discussed more completely in Chapter 14.

5.2.3 Calculating shading

When implementing surface shading, the code needs to have access to information about the light source, the surface, and the viewing direction. Writing clean code that supports both point and directional lights is easiest to do by separating the calculation of irradiance from the calculation of reflected light. Irradiance depends only on the light source and the surface geometry, and once it’s known, calculating the reflected light only depends on the surface properties and the viewing geometry.

Basic shading calculations can be done in exactly the same way in ray tracing and rasterization systems; it’s really only how the inputs are computed that varies. To compute irradiance, we need

• The shading point x, a 3D point on a surface

• The surface normal n perpendicular to the surface at x

• The light source position p for a point light or its direction l for a directional light

• The light source intensity I for a point light or its normal irradiance H for a directional light (these are RGB colors).

For a point light, we need to compute the distance and the light direction, which are both simple to get from the vector p – x:

and for both types of lights, the cosine factor is best computed using a dot product; as long as n and l are unit vectors,

In practice, it’s a good idea when computing irradiance to clamp the dot product at zero to make sure that even if in some cases, you find the light direction is facing away from the surface normal, you won’t get negative shading. This leads to what we view as the official equations for computing irradiance:

Once the irradiance is known, it needs to be multiplied by the BDRF value, and the ingredients for calculating that value are

• The light direction l, a unit vector pointing from x toward the light (already computed as part of the irradiance calculation)

• The viewing direction v, a unit vector pointing from x toward the viewer

• The parameters describing the properties of the surface material. For this chapter’s model, this includes R, k_s, and p.

How you get these quantities differs substantially between ray tracing and rasterization systems, but the actual shading calculation itself is the same. Don’t forget that v, l, and n all must be unit vectors; failing to normalize these vectors is a very common error in shading computations.

6.2.1 Matrix Arithmetic

A matrix times a constant results in a matrix where each element has been multiplied by that constant, e.g.,

For matrix multiplication, we “multiply” rows of the first matrix with columns of the second matrix:

For matrix multiplication, we “multiply” rows of the first matrix with columns of the second matrix:

So the element p_ij of the resulting product is

Taking a product of two matrices is only possible if the number of columns of the left matrix is the same as the number of rows of the right matrix. For example,

Matrix multiplication is not commutative in most instances:

Also, if AB = AC, it does not necessarily follow that B = C. Fortunately, matrix multiplication is associative and distributive:

6.2.2 Operations on Matrices

We would like a matrix analog of the inverse of a real number. We know the inverse of a real number x is 1/x and that the product of x and its inverse is 1. We need a matrix I that we can think of as a “matrix one.” This exists only for square matrices and is known as the identity matrix; it consists of ones down the diagonal and zeroes elsewhere. For example, the four by four identity matrix is

The inverse matrix A^–1 of a matrix A is the matrix that ensures AA^– 1 = I. For example,

Note that the inverse of A^-1 is A. So AA^-1 = A^-1A = I. The inverse of a product of two matrices is the product of the inverses, but with the order reversed:

We will return to the question of computing inverses in Section 6.3.

The transpose A^T of a matrix A has the same numbers, but the rows are switched with the columns. If we label the entries of A^T as a_ij,then

For example,

The transpose of a product of two matrices obeys a rule similar to Equation (6.4):

The determinant of a square matrix is simply the determinant of the columns of the matrix, considered as a set of vectors. The determinant has several nice relationships to the matrix operations just discussed, which we list here for reference:

6.2.3 Vector Operations in Matrix Form

In graphics, we use a square matrix to transform a vector represented as a matrix. For example, if you have a 2D vector a = (x_a, y_a) and want to rotate it by 90 degrees about the origin to form vector a = (–y_a, x_a) , you can use a product of a 2 × 2 matrix and a 2 × 1 matrix, called a column vector. The operation in matrix form is

We can get the same result by using the transpose of this matrix and multiplying on the left (“premultiplying”) with a row vector:

These days, postmultiplication using column vectors is fairly standard, but in many older books and systems, you will run across row vectors and premultiplication. The only difference is that the transform matrix must be replaced with its transpose.

We also can use matrix formalism to encode operations on just vectors. If we consider the result of the dot product as a 1 × 1 matrix, it can be written

For example, if we take two 3D vectors we get

A related vector product is the outer product between two vectors, which can be expressed as a matrix multiplication with a column vector on the left and a row vector on the right: ab^T. The result is a matrix consisting of products of all pairs of an entry of a with an entry of b. For 3D vectors, we have

It is often useful to think of matrix multiplication in terms of vector operations. To illustrate using the three-dimensional case, we can think of a 3 × 3 matrix as a collection of three 3D vectors in two ways: either it is made up of three column vectors side-by-side or it is made up of three row vectors stacked up. For instance, the result of a matrix-vector multiplication y = Ax can be interpreted as a vector whose entries are the dot products of x with the rows of A. Naming these row vectors r_i, we have

Alternatively, we can think of the same product as a sum of the three columns c_i of A, weighted by the entries of x:

Using the same ideas, one can understand a matrix–matrix product AB as an array containing the pairwise dot products of all rows of A with all columns of B (cf. (6.2)); as a collection of products of the matrix A with all the column vectors of B, arranged left to right; as a collection of products of all the row vectors of A with the matrix B, stacked top to bottom; or as the sum of the pairwise outer products of all columns of A with all rows of B. (See Exercise 8.)

These interpretations of matrix multiplication can often lead to valuable geometric interpretations of operations that may otherwise seem very abstract.

6.2.4 Special Types of Matrices

The identity matrix is an example of a diagonal matrix, where all nonzero elements occur along the diagonal. The diagonal consists of those elements whose column index equals the row index counting from the upper left.

The identity matrix also has the property that it is the same as its transpose. Such matrices are called symmetric.

The identity matrix is also an orthogonal matrix, because each of its columns considered as a vector has length 1 and the columns are orthogonal to one another. The same is true of the rows (see Exercise 2). The determinant of any orthogonal matrix is either +1 or –1.

A very useful property of orthogonal matrices is that they are nearly their own inverses. Multiplying an orthogonal matrix by its transpose results in the identity,

This is easy to see because the entries of R^TR are dot products between the columns of R. Off-diagonal entries are dot products between orthogonal vectors, and the diagonal entries are dot products of the (unit-length) columns with themselves.

6.3.1 Computing Inverses

Determinants give us a tool to compute the inverse of a matrix. It is a very inefficient method for large matrices, but often in graphics, our matrices are small. A key to developing this method is that the determinant of a matrix with two identical rows is zero. This should be clear because the volume of the n-dimensional parallelepiped is zero if two of its sides are the same. Suppose we have a 4 × 4 A and we wish to find its inverse A^–1. The inverse is

Note that this is just the transpose of the matrix where elements of A are replaced by their respective cofactors multiplied by the leading constant (1 or –1). This matrix is called the adjoint of A. The adjoint is the transpose of the cofactor matrix of A. We can see why this is an inverse. Look at the product AA^–1 which we expect to be the identity. If we multiply the first row of A by the first column of the adjoint matrix we need to get |A| (remember the leading constant above divides by |A|:

This is true because the elements in the first row of A are multiplied exactly by their cofactors in the first column of the adjoint matrix which is exactly the determinant. The other values along the diagonal of the resulting matrix are |A| for analogous reasons. The zeros follow a similar logic:

Note that this product is a determinant of some matrix:

The matrix in fact is

Because the first two rows are identical, the matrix is singular, and thus, its determinant is zero.

The argument above does not apply just to four by four matrices; using that size just simplifies typography. For any matrix, the inverse is the adjoint matrix divided by the determinant of the matrix being inverted. The adjoint is the transpose of the cofactor matrix, which is just the matrix whose elements have been replaced by their cofactors.

6.3.2 Linear Systems

We often encounter linear systems in graphics with “n equations and n unknowns,” usually for n = 2 or n = 3. For example,

Here, x, y, and z are the “unknowns” for which we wish to solve. We can write this in matrix form:

A common shorthand for such systems is Ax = b where it is assumed that A is a square matrix with known constants, x is an unknown column vector (with elements x, y,and z in our example), and b is a column matrix of known constants.

There are many ways to solve such systems, and the appropriate method depends on the properties and dimensions of the matrix A. Because in graphics we so frequently work with systems of size n ≤ 4, we’ll discuss here a method appropriate for these systems, known as Cramer’s rule, which we saw earlier, from a 2D geometric viewpoint, in the example on page 108. Here, we show this algebraically. The solution to the above equation is

The rule here is to take a ratio of determinants, where the denominator is |A| and the numerator is the determinant of a matrix created by replacing a column of A with the column vector b. The column replaced corresponds to the position of the unknown in vector x. For example, y is the second unknown and the second column is replaced. Note that if |A| = 0, the division is undefined and there is no solution. This is just another version of the rule that if A is singular (zero determinant), then there is no unique solution to the equations.

6.4.1 Singular Value Decomposition

We saw in the last section that any symmetric matrix can be diagonalized, or decomposed into a convenient product of orthogonal and diagonal matrices. However, most matrices we encounter in graphics are not symmetric, and the eigenvalue decomposition for nonsymmetric matrices is not nearly so convenient or illuminating, and in general involves complex-valued eigenvalues and eigenvectors even for real-valued inputs.

There is another generalization of the symmetric eigenvalue decomposition to nonsymmetric (and even non-square) matrices; it is the singular value decomposition (SVD). The main difference between the eigenvalue decomposition of a symmetric matrix and the SVD of a nonsymmetric matrix is that the orthogonal matrices on the left and right sides are not required to be the same in the SVD:

Here, U and V are two, potentially different, orthogonal matrices, whose columns are known as the left and right singular vectors of A,and S is a diagonal matrix whose entries are known as the singular values of A. When A is symmetric and has all nonnegative eigenvalues, the SVD and the eigenvalue decomposition are the same.

There is another relationship between singular values and eigenvalues that can be used to compute the SVD (though this is not the way an industrial-strength SVD implementation works). First, we define M = AA^T. We assume that we can perform a SVD on M:

The substitution is based on the fact that (BC)^T = C^TB^T, that the transpose of an orthogonal matrix is its inverse, and that the transpose of a diagonal matrix is the matrix itself. The beauty of this new form is that M is symmetric and US²U^T is its eigenvalue decomposition, where S² contains the (all nonnegative) eigenvalues. Thus, we find that the singular values of a matrix are the square roots of the eigenvalues of the product of the matrix with its transpose, and the left singular vectors are the eigenvectors of that product. A similar argument allows V, the matrix of right singular vectors, to be computed from A^TA.

This form used the standard symbol σ_i for the ith singular value. Again, for a symmetric matrix, the eigenvalues and the singular values are the same (σ_i = λ_i). We will examine the geometry of SVD further in Section 7.1.6.

7.1.1 Scaling

The most basic transform is a scale along the coordinate axes. This transform can change length and possibly direction:

Note what this matrix does to a vector with Cartesian components (x, y) :

So, just by looking at the matrix of an axis-aligned scale, we can read off the two scale factors.

7.1.2 Shearing

A shear is something that pushes things sideways, producing something like a deck of cards across which you push your hand; the bottom card stays put and cards move more the closer they are to the top of the deck. The horizontal and vertical shear matrices are

An analogous transform vertically is (see Figure 7.4)

In both cases, the square outline of the sheared clock becomes a parallelogram, and the circular face of the sheared clock becomes an ellipse.

Another way to think of a shear is in terms of rotation of only the vertical (or horizontal) axis. The shear transform that takes a vertical axis and tilts it clockwise by an angle ϕ is

Similarly, the shear matrix which rotates the horizontal axis counterclockwise by angle ϕ is

7.1.3 Rotation

Suppose we want to rotate a vector a by an angle ϕ counterclockwise to get vector b (Figure 7.5). If a makes an angle α with the x-axis, and its length is , then we know that

Because b is a rotation of a, it also has length r. Because it is rotated an angle ϕ from a, b makes an angle (α + ϕ) with the x-axis. Using the trigonometric addition identities (Section 2.3.3):

Substituting x_a = r cos α and y_a = r sin α gives

In matrix form, the transformation that takes a to b is then

Because the norm of each row of a rotation matrix is one (sin² ϕ+cos² ϕ = 1), and the rows are orthogonal (cos ϕ(– sin ϕ)+sin ϕ cos ϕ = 0), we see that rotation matrices are orthogonal matrices (Section 6.2.4). By looking at the matrix, we can read off two pairs of orthonormal vectors: the two columns, which are the vectors to which the transformation sends the canonical basis vectors (1, 0) and (0, 1) ; and the rows, which are the vectors that the transformations sends to the canonical basis vectors.

7.1.4 Reflection

We can reflect a vector across either of the coordinate axes by using a scale with one negative scale factor (see Figures 7.8 and 7.9):

While one might expect that the matrix with –1 in both elements of the diagonal is also a reflection, in fact it is just a rotation by π radians.

7.1.5 Composition and Decomposition of Transformations

It is common for graphics programs to apply more than one transformation to an object. For example, we might want to first apply a scale S and then a rotation R. This would be done in two steps on a 2D vector v₁:

Another way to write this is

Because matrix multiplication is associative, we can also write

In other words, we can represent the effects of transforming a vector by two matrices in sequence using a single matrix of the same size, which we can compute by multiplying the two matrices: M = RS (Figure 7.10).

It is very important to remember that these transforms are applied from the right side first. So the matrix M = RS first applies S and then R.

Building up a transformation from rotation and scaling transformations actually works for any linear transformation, and this fact leads to a powerful way of thinking about these transformations, as explored in the next section.

7.1.6 Decomposition of Transformations

Sometimes, it’s necessary to “undo” a composition of transformations, taking a transformation apart into simpler pieces. For instance, it’s often useful to present a transformation to the user for manipulation in terms of separate rotations and scale factors, but a transformation might be represented internally simply as a matrix, with the rotations and scales already mixed together. This kind of manipulation can be achieved if the matrix can be computationally disassembled into the desired pieces, the pieces adjusted, and the matrix reassembled by multiplying the pieces together again.

It turns out that this decomposition, or factorization, is possible, regardless of the entries in the matrix—and this fact provides a fruitful way of thinking about transformations and what they do to geometry that is transformed by them.

Symmetric Eigenvalue Decomposition

Let’s start with symmetric matrices. Recall from Section 6.4 that a symmetric matrix can always be taken apart using the eigenvalue decomposition into a product of the form

where R is an orthogonal matrix and S is a diagonal matrix; we will call the columns of R (the eigenvectors) by the names v₁ and v₂, and we’ll call the diagonal entries of S (the eigenvalues) by the names λ₁ and λ₂.

In geometric terms, we can now recognize R as a rotation and S as a scale, so this is just a multi-step geometric transformation (Figure 7.12):

1. Rotate v₁ and v₂ to the x- and y-axes (the transform by R^T).

2. Scale in x and y by (λ₁,λ₂) (the transform by S).

3. Rotate the x- and y-axes back to v₁ and v₂ (the transform by R).

Looking at the effect of these three transforms together, we can see that they have the effect of a nonuniform scale along a pair of axes. As with an axis-aligned scale, the axes are perpendicular, but they aren’t the coordinate axes; instead, they are the eigenvectors of A. This tells us something about what it means to be a symmetric matrix: symmetric matrices are just scaling operations—albeit potentially nonuniform and non–axis-aligned ones.

Singular Value Decomposition

A very similar kind of decomposition can be done with nonsymmetric matrices as well: it’s the singular value decomposition (SVD), also discussed in Section 6.4.1. The difference is that the matrices on either side of the diagonal matrix are no longer the same:

The two orthogonal matrices that replace the single rotation R are called U and V, and their columns are called u_i (the left singular vectors) and v_i (the right singular vectors), respectively. In this context, the diagonal entries of S are called singular values rather than eigenvalues. The geometric interpretation is very similar to that of the symmetric eigenvalue decomposition (Figure 7.14):

1. Rotate v₁ and v₂ to the x- and y-axes (the transform by V^T).

2. Scale in x and y by (σ₁,σ₂) (the transform by S).

3. Rotate the x- and y-axes to u₁ and u₂ (the transform by U).

The principal difference is between a single rotation and two different orthogonal matrices. This difference causes another, less important, difference. Because the SVD has different singular vectors on the two sides, there is no need for negative singular values: we can always flip the sign of a singular value, reverse the direction of one of the associated singular vectors, and end up with the same transformation again. For this reason, the SVD always produces a diagonal matrix with all positive entries, but the matrices U and V are not guaranteed to be rotations—they could include reflection as well. In geometric applications like graphics, this is an inconvenience, but a minor one: it is easy to differentiate rotations from reflections by checking the determinant, which is +1 for rotations and –1 for reflections, and if rotations are desired, one of the singular values can be negated, resulting in a rotation–scale–rotation sequence where the reflection is rolled in with the scale, rather than with one of the rotations.

In summary, every matrix can be decomposed via SVD into a rotation times a scale times another rotation. Only symmetric matrices can be decomposed via eigenvalue diagonalization into a rotation times a scale times the inverse-rotation, and such matrices are a simple scale in an arbitrary direction. The SVD of a symmetric matrix will yield the same triple product as eigenvalue decomposition via a slightly more complex algebraic manipulation.

Paeth Decomposition of Rotations

Another decomposition uses shears to represent nonzero rotations (Paeth, 1990). The following identity allows this:

For example, a rotation by π/4 (45°) is (see Figure 7.16)

This particular transform is useful for raster rotation because shearing is a very efficient raster operation for images; it introduces some jagginess, but will leave no holes. The key observation is that if we take a raster position (i, j) and apply a horizontal shear to it, we get

If we round sj to the nearest integer, this amounts to taking each row in the image and moving it sideways by some amount—a different amount for each row. Because it is the same displacement within a row, this allows us to rotate with no gaps in the resulting image. A similar action works for a vertical shear. Thus, we can implement a simple raster rotation easily.

7.2.1 Arbitrary 3D Rotations

As in 2D, 3D rotations are orthogonal matrices. Geometrically, this means that the three rows of the matrix are the Cartesian coordinates of three mutually orthogonal unit vectors as discussed in Section 2.4.5. The columns are three, potentially different, mutually orthogonal unit vectors. There are an infinite number of such rotation matrices. Let’s write down such a matrix:

Here, u = x_ux + y_uy + z_uz and so on for v and w. Since the three vectors are orthonormal, we know that

We can infer some of the behavior of the rotation matrix by applying it to the vectors u, v and w. For example,

Note that those three rows of R_uvwu are all dot products:

Similarly, R_uvwv = y, and R_uvww = z. So R_uvw takes the basis uvw to the corresponding Cartesian axes via rotation.

If R_uvw is a rotation matrix with orthonormal rows, then R^T_uvw is also a rotation matrix with orthonormal columns and in fact is the inverse of R_uvw (the inverse of an orthogonal matrix is always its transpose). An important point is that for transformation matrices, the algebraic inverse is also the geometric inverse. So if R_uvw takes u to x, then R^T_uvw takes x to u. The same should be true of v and y as we can confirm:

So we can always create rotation matrices from orthonormal bases.

If we wish to rotate about an arbitrary vector a, we can form an orthonormal basis with w = a, rotate that basis to the canonical basis xyz, rotate about the z-axis, and then rotate the canonical basis back to the uvw basis. In matrix form, to rotate about the w-axis by an angle ϕ:

Here, we have w a unit vector in the direction of a (i.e., a divided by its own length). But what are u and v? A method to find reasonable u and v is given in Section 2.4.6.

If we have a rotation matrix and we wish to have the rotation in axis-angle form, we can compute the one real eigenvalue (which will be λ = 1), and the corresponding eigenvector is the axis of rotation. This is the one axis that is not changed by the rotation.

See Section 16.2.2 for a comparison of the few most-used ways to represent rotations, besides rotation matrices.

7.2.2 Transforming Normal Vectors

While most 3D vectors we use represent positions (offset vectors from the origin) or directions, such as where light comes from, some vectors represent surface normals. Surface normal vectors are perpendicular to the tangent plane of a surface. These normals do not transform the way we would like when the underlying surface is transformed. For example, if the points of a surface are transformed by a matrix M, a vector t that is tangent to the surface and is multiplied by M will be tangent to the transformed surface. However, a surface normal vector n that is transformed by M may not be normal to the transformed surface (Figure 7.17).

We can derive a transform matrix N which does take n to a vector perpendicular to the transformed surface. One way to attack this issue is to note that a surface normal vector and a tangent vector are perpendicular, so their dot product is zero, which is expressed in matrix form as

If we denote the desired transformed vectors as t_M = Mt and n_N = Nn, our goal is to find N such that . We can find N by some algebraic tricks.

First, we can sneak an identity matrix into the dot product and then take advantage of M^–1M = I:

Although the manipulations above don’t obviously get us anywhere, note that we can add parentheses that make the above expression more obviously a dot product:

This means that the row vector that is perpendicular to t_M is the left part of the expression above. This expression holds for any of the tangent vectors in the tangent plane. Since there is only one direction in 3D (and its opposite) that is perpendicular to all such tangent vectors, we know that the left part of the expression above must be the row vector expression for n_N ; i.e., it is n^T_N , so this allows us to infer N:

so we can take the transpose of that to get

Therefore, we can see that the matrix that correctly transforms normal vectors so they remain normal is N = (M^–1)^T, i.e., the transpose of the inverse matrix. Since this matrix may change the length of n, we can multiply it by an arbitrary scalar and it will still produce n_N with the right direction. Recall from Section 6.3 that the inverse of a matrix is the transpose of the cofactor matrix divided by the determinant. Because we don’t care about the length of a normal vector, we can skip the division and find that for a 3 × 3 matrix,

This assumes the element of M in row i and column j is m_ij. So the full expression for N is

8.1.1 The Viewport Transformation

We begin with a problem whose solution will be reused for any viewing condition. We assume that the geometry we want to view is in the canonical view volume, and we wish to view it with an orthographic camera looking in the – z direction. The canonical view volume is the cube containing all 3D points whose Cartesian coordinates are between –1 and +1—that is, (x, y, z) ∈ [–1, 1]³ (Figure 8.3). We project x = –1 to the left side of the screen, x = +1 to the right side of the screen, y = –1 to the bottom of the screen, and y = +1 to the top of the screen.

Recall the conventions for pixel coordinates from Chapter 3: each pixel “owns” a unit square centered at integer coordinates; the image boundaries have a half-unit overshoot from the pixel centers; and the smallest pixel center coordinates are (0, 0) . If we are drawing into an image (or window on the screen) that has n_x by n_y pixels, we need to map the square [–1, 1]² to the rectangle [–0.5,n_x – 0.5] × [–0.5,n_y – 0.5].

For now, we will assume that all line segments to be drawn are completely inside the canonical view volume. Later, we will relax that assumption when we discuss clipping.

Since the viewport transformation maps one axis-aligned rectangle to another, it is a case of the windowing transform given by Equation (7.6):

Note that this matrix ignores the z-coordinate of the points in the canonical view volume, because a point’s distance along the projection direction doesn’t affect where that point projects in the image. But before we officially call this the view-port matrix, we add a row and column to carry along the z-coordinate without changing it. We don’t need it in this chapter, but eventually, we will need the z values because they can be used to make closer surfaces hide more distant surfaces (see Section 9.2.3).

8.1.2 The Orthographic Projection Transformation

Of course, we usually want to render geometry in some region of space other than the canonical view volume. Our first step in generalizing the view will keep the view direction and orientation fixed looking along – z with +y up, but will allow arbitrary rectangles to be viewed. Rather than replacing the viewport matrix, we’ll augment it by multiplying it with another matrix on the right.

Under these constraints, the view volume is an axis-aligned box, and we’ll name the coordinates of its sides so that the view volume is [l, r] × [b, t] × [f, n] shown in Figure 8.4. We call this box the orthographic view volume and refer to the bounding planes as follows:

That vocabulary assumes a viewer who is looking along the minus z-axis with his head pointing in the y-direction.¹ This implies that n > f, which may be unintuitive, but if you assume the entire orthographic view volume has negative z values, then the z = n “near” plane is closer to the viewer if and only if n > f ; here, f is a smaller number than n, i.e., a negative number of larger absolute value than n.

This concept is shown in Figure 8.5. The transform from orthographic view volume to the canonical view volume is another windowing transform, so we can simply substitute the bounds of the orthographic and canonical view volumes into Equation (7.7) to obtain the matrix for this transformation:

To draw 3D line segments in the orthographic view volume, we project them into screen x-and y-coordinates and ignore z-coordinates. We do this by combining Equations (8.2) and (8.3). Note that in a program, we multiply the matrices together to form one matrix and then manipulate points as follows:

The z-coordinate will now be in [–1, 1]. We don’t take advantage of this now, but it will be useful when we examine z-buffer algorithms.

The code to draw many 3D lines with endpoints a_i and b_i thus becomes both simple and efficient:

construct Mvp
construct Morth
M = MvpMorth
for each line segment (ai, bi) do
  p = Mai
  q = Mbi
  drawline(xp, yp, xq, yq)

8.1.3 The Camera Transformation

We’d like to be able to change the viewpoint in 3D and look in any direction. There are a multitude of conventions for specifying viewer position and orientation. We will use the following one (see Figure 8.6):

• the eye position e,

• the gaze direction g,

• the view-up vector t.

The eye position is a location that the eye “sees from.” If you think of graphics as a photographic process, it is the center of the lens. The gaze direction is any vector in the direction that the viewer is looking. The view-up vector is any vector in the plane that both bisects the viewer’s head into right and left halves and points “to the sky” for a person standing on the ground. These vectors provide us with enough information to set up a coordinate system with origin e and a uvw basis, using the construction of Section 2.4.7:

Our job would be done if all points we wished to transform were stored in coordinates with origin e and basis vectors u, v,and w. But as shown in Figure 8.7, the coordinates of the model are stored in terms of the canonical (or world) origin o and the x-, y-, and z-axes. To use the machinery we have already developed, we just need to convert the coordinates of the line segment endpoints we wish to draw from xyz-coordinates into uvw-coordinates. This kind of transformation was discussed in Section 7.5, and the matrix that enacts this transformation is the canonical-to-basis matrix of the camera’s coordinate frame:

Alternatively, we can think of this same transformation as first moving e to the origin, then aligning u, v, w to x, y, z.

To make our previously z-axis-only viewing algorithm work for cameras with any location and orientation, we just need to add this camera transformation to the product of the viewport and projection transformations, so that it converts the incoming points from world to camera coordinates before they are projected:

construct Mvp
construct Morth
construct Mcam
M = MvpMorthMcam
for each line segment (ai, bi) do
  p = Mai
  q = Mbi
  drawline(xp, yp, xq, yq)

Again, almost no code is needed once the matrix infrastructure is in place.

9.1.1 Line Drawing

Most graphics packages contain a line drawing command that takes two endpoints in screen coordinates (see Figure 3.10) and draws a line between them. For example, the call for endpoints (1,1) and (3,2) would turn on pixels (1,1) and (3,2) and fill in one pixel between them. For general screen coordinate endpoints (x₀,y₀) and (x₁,y₁), the routine should draw some “reasonable” set of pixels that approximates a line between them. Drawing such lines is based on line equations, and we have two types of equations to choose from: implicit and parametric. This section describes the approach using implicit lines.

Line Drawing Using Implicit Line Equations

The most common way to draw lines using implicit equations is the midpoint algorithm ((Pitteway, 1967; van Aken & Novak, 1985)). The midpoint algorithm ends up drawing the same lines as the Bresenham algorithm (Bresenham, 1965), but it is somewhat more straightforward.

The first thing to do is find the implicit equation for the line as discussed in Section 2.7.2:

We assume that x₀ ≤ x₁. If that is not true, we swap the points so that it is true. The slope m of the line is given by

The following discussion assumes m ∈ (0, 1]. Analogous discussions can be derived for m ∈ (–∞, –1], m ∈ (–1, 0],and m ∈ (1, ∞). The four cases cover all possibilities.

For the case m ∈ (0, 1], there is more “run” than “rise” ; i.e., the line is moving faster in x than in y. If we have an API where the y-axis points downward, we might have a concern about whether this makes the process harder, but, in fact, we can ignore that detail. We can ignore the geometric notions of “up” and “down,” because the algebra is exactly the same for the two cases. Cautious readers can confirm that the resulting algorithm works for the y-axis downward case. The key assumption of the midpoint algorithm is that we draw the thinnest line possible that has no gaps. A diagonal connection between two pixels is not considered a gap.

As the line progresses from the left endpoint to the right, there are only two possibilities: draw a pixel at the same height as the pixel drawn to its left, or draw a pixel one higher. There will always be exactly one pixel in each column of pixels between the endpoints. Zero would imply a gap, and two would be too thick a line. There may be two pixels in the same row for the case we are considering; the line is more horizontal than vertical, so sometimes it will go right and sometimes up. This concept is shown in Figure 9.2, where three “reasonable” lines are shown, each advancing more in the horizontal direction than in the vertical direction.

The midpoint algorithm for m ∈ (0, 1] first establishes the leftmost pixel and the column number (x-value) of the rightmost pixel and then loops horizontally establishing the row (y-value) of each pixel. The basic form of the algorithm is

y = y0
for x = x0 to x1 do
  draw(x, y)
  if (some condition) then
     y = y +1

Note that x and y are integers. In words this says, “keep drawing pixels from left to right and sometimes move upward in the y-direction while doing so.” The key is to establish efficient ways to make the decision in the if statement.

An effective way to make the choice is to look at the midpoint of the line between the two potential pixel centers. More specifically, the pixel just drawn is pixel (x, y) whose center in real screen coordinates is at (x, y). The candidate pixels to be drawn to the right are pixels (x + 1,y) and (x + 1,y + 1). The midpoint between the centers of the two candidate pixels is (x +1,y +0.5). If the line passes below this midpoint, we draw the bottom pixel, and otherwise, we draw the top pixel (Figure 9.3).

To decide whether the line passes above or below (x+1,y +0.5),weevaluate f (x +1,y +0.5) in Equation (9.1). Recall from Section 2.7.1 that f (x, y) = 0 for points (x, y) on the line, f (x, y) > 0 for points on one side of the line, and f (x, y) < 0 for points on the other side of the line. Because –f (x, y) = 0 and f (x, y) = 0 are both perfectly good equations for the line, it is not immediately clear whether f (x, y) being positive indicates that (x, y) is above the line or whether it is below. However, we can figure it out; the key term in Equation (9.1) is the y term (x₁ – x₀)y . Note that (x₁ – x₀) is definitely positive because x₁ > x₀. This means that as y increases, the term (x₁ – x₀)y gets larger (i.e., more positive or less negative). Thus, the case f (x, +∞) is definitely positive, and definitely above the line, implying points above the line are all positive.

Another way to look at it is that the y component of the gradient vector is positive. So above the line, where y can increase arbitrarily, f (x, y) must be positive. This means we can make our code more specific by filling in the if statement:

if f (x + 1,y + 0.5) < 0 then
   y = y + 1

The above code will work nicely for lines of the appropriate slope (i.e., between zero and one). The reader can work out the other three cases which differ only in small details.

If greater efficiency is desired, using an incremental method can help. An incremental method tries to make a loop more efficient by reusing computation from the previous step. In the midpoint algorithm as presented, the main computation is the evaluation of f (x +1,y +0.5). Note that inside the loop, after the first iteration, either we already evaluated f (x – 1,y +0.5) or f (x – 1,y – 0.5) (Figure 9.4). Note also this relationship:

This allows us to write an incremental version of the code:

y = y0
d = f (x0 + 1,y0 + 0.5)
for x = x0 to x1 do
   draw(x, y)
   if d < 0 then
      y = y + 1
      d = d + (x1 – x0) + (y0 – y1)
   else
      d = d + (y0 – y1)

This code should run faster since it has little extra setup cost compared to the non-incremental version (that is not always true for incremental algorithms), but it may accumulate more numeric error because the evaluation of f (x, y + 0.5) may be composed of many adds for long lines. However, given that lines are rarely longer than a few thousand pixels, such an error is unlikely to be critical. Slightly longer setup cost, but faster loop execution, can be achieved by storing (x₁ – x₀)+(y₀ – y₁) and (y₀ – y₁) as variables. We might hope a good compiler would do that for us, but if the code is critical, it would be wise to examine the results of compilation to make sure.

9.1.2 Triangle Rasterization

We often want to draw a 2D triangle with 2D points p₀ = (x₀,y₀), p₁ = (x₁,y₁), and p₂ = (x₂,y₂) in screen coordinates. This is similar to the line drawing problem, but it has some of its own subtleties. As with line drawing, we may wish to interpolate color or other properties from values at the vertices. This is straightforward if we have the barycentric coordinates (Section 2.9). For example, if the vertices have colors c₀, c₁, and c₂, the color at a point in the triangle with barycentric coordinates (α, β, γ) is

This type of interpolation of color is known in graphics as Gouraud interpolation after its inventor (Gouraud, 1971).

Another subtlety of rasterizing triangles is that we are usually rasterizing triangles that share vertices and edges. This means we would like to rasterize adjacent triangles, so there are no holes. We could do this by using the midpoint algorithm to draw the outline of each triangle and then fill in the interior pixels. This would mean adjacent triangles both draw the same pixels along each edge. If the adjacent triangles have different colors, the image will depend on the order in which the two triangles are drawn. The most common way to rasterize triangles that avoids the order problem and eliminates holes is to use the convention that pixels are drawn if and only if their centers are inside the triangle; i.e., the barycentric coordinates of the pixel center are all in the interval (0, 1). Thisraises the issue of what to do if the center is exactly on the edge of the triangle. There are several ways to handle this as will be discussed later in this section. The key observation is that barycentric coordinates allow us to decide whether to draw a pixel and what color that pixel should be if we are interpolating colors from the vertices. So our problem of rasterizing the triangle boils down to efficiently finding the barycentric coordinates of pixel centers (Pineda, 1988). The brute force rasterization algorithm is

for all x do
  for all y do
    compute (α, β, γ) for (x, y)
    if (α ∈ [0, 1] and β ∈ [0, 1] and γ ∈ [0, 1]) then
       c = αc0 + βc1 + γc2
       drawpixel (x, y) with color c

The rest of the algorithm limits the outer loops to a smaller set of candidate pixels and makes the barycentric computation efficient.

We can add a simple efficiency by finding the bounding rectangle of the three vertices and only looping over this rectangle for candidate pixels to draw. We can compute barycentric coordinates using Equation (2.32). This yields the algorithm:

xmin = floor(xi)
xmax = ceiling(xi)
ymin = floor(yi)
ymax = ceiling(yi)
for y = ymin to ymax do
  for x = xmin to xmax do
    α = f12(x, y)/f12(x0,y0)
    β = f20(x, y)/f20(x1,y1)
    γ = f01(x, y)/f01(x2,y2)
    if (α > 0 and β > 0 and γ > 0) then
       c = αc0 + βc1 + γc2
       drawpixel (x, y) with color c

Here, f_ij is the line given by Equation (9.1) with the appropriate vertices:

Note that we have exchanged the test α ∈ (0, 1) with α > 0 etc., because if all of α, β, γ are positive, then we know they are all less than one because α + β + γ = 1. We could also compute only two of the three barycentric variables and get the third from that relation, but it is not clear that this saves computation once the algorithm is made incremental, which is possible as in the line drawing algorithms; each of the computations of α, β, and γ does an evaluation of the form f (x, y) = Ax + By + C. In the inner loop, only x changes, and it changes by one. Note that f (x +1,y) = f (x, y)+ A. This is the basis of the incremental algorithm. In the outer loop, the evaluation changes for f (x, y) to f (x, y + 1), so a similar efficiency can be achieved. Because α, β, and γ change by constant increments in the loop, so does the color c. So this can be made incremental as well. For example, the red value for pixel (x + 1,y) differs from the red value for pixel (x, y) by a constant amount that can be precomputed. An example of a triangle with color interpolation is shown in Figure 9.5.

Dealing with Pixels on Triangle Edges

We have still not discussed what to do for pixels whose centers are exactly on the edge of a triangle. If a pixel is exactly on the edge of a triangle, then it is also on the edge of the adjacent triangle if there is one. There is no obvious way to award the pixel to one triangle or the other. The worst decision would be to not draw the pixel because a hole would result between the two triangles. Better, but still not good, would be to have both triangles draw the pixel. If the triangles are transparent, this will result in a double-coloring. We would really like to award the pixel to exactly one of the triangles, and we would like this process to be simple; which triangle is chosen does not matter as long as the choice is well defined.

One approach is to note that any off-screen point is definitely on exactly one side of the shared edge and that is the edge we will draw. For two non-overlapping triangles, the vertices not on the edge are on opposite sides of the edge from each other. Exactly one of these vertices will be on the same side of the edge as the off-screen point (Figure 9.6). This is the basis of the test. The test if numbers p and q have the same sign can be implemented as the test pq > 0, which is very efficient in most environments.

Note that the test is not perfect because the line through the edge may also go through the off-screen point, but we have at least greatly reduced the number of problematic cases. Which off-screen point is used is arbitrary, and (x, y) = (–1, –1) is as good a choice as any. We will need to add a check for the case of a point exactly on an edge. We would like this check not to be reached for common cases, which are the completely inside or outside tests. This suggests

fα = f12(x0,y0)
fβ = f20(x1,y1)
fγ = f01(x2,y2)
for y = ymin to ymax do
  for x = xmin to xmax do
    α = f12(x, y)/fα
    β = f20(x, y)/fβ
    γ = f01(x, y)/fγ
    if (α ≥ 0 and β ≥ 0 and γ ≥ 0) then
       if (α > 0 or fαf12(–1, –1) > 0) and
          (β > 0 or fβf20(–1, –1) > 0) and
          (γ > 0 or fγf01(–1, –1) > 0) then
          c = αc0 + βc1 + γc2
          drawpixel (x, y) with color c

We might expect that the above code would work to eliminate holes and double-draws only if we use exactly the same line equation for both triangles. In fact, the line equation is the same only if the two shared vertices have the same order in the draw call for each triangle. Otherwise, the equation might flip in sign. This could be a problem depending on whether the compiler changes the order of operations. So if a robust implementation is needed, the details of the compiler and arithmetic unit may need to be examined. The first four lines in the pseudocode above must be coded carefully to handle cases where the edge exactly hits the pixel center.

In addition to being amenable to an incremental implementation, there are several potential early exit points. For example, if α is negative, there is no need to compute β or γ. While this may well result in a speed improvement, profiling is always a good idea; the extra branches could reduce pipelining or concurrency and might slow down the code. So as always, test any attractive-looking optimizations if the code is a critical section.

Another detail of the above code is that the divisions could be divisions by zero for degenerate triangles, i.e., if f_γ = 0. Either the floating point error conditions should be accounted for properly, or another test will be needed.

9.1.3 Perspective Correct Interpolation

There are some subtleties in achieving correct-looking perspective when interpolating quantities, such as texture coordinates or 3D positions, that need to vary linearly across the 3D triangles. We’ll use texture coordinates as an example of a quantity where perspective correction is important, but the same considerations apply to any attribute where linearity in 3D space is important.

The reason things are not straightforward is that just interpolating texture coordinates in screen space results in incorrect images, as shown for the grid texture in Figure 9.7. Because things in perspective get smaller as the distance to the viewer increases, the lines that are evenly spaced in 3D should compress in 2D image space. More careful interpolation of texture coordinates is needed to accomplish this.

We can implement texture mapping on triangles by interpolating the (u, v) coordinates, modifying the rasterization method of Section 9.1.2, but this results in the problem shown at the right of Figure 9.7. A similar problem occurs for triangles if screen space barycentric coordinates are used as in the following rasterization code:

for all x do
  for all y do
    compute (α, β, γ) for (x, y)
    if α ∈ (0, 1) and β ∈ (0, 1) and γ ∈ (0, 1) then
       t = αt0 + βt1 + γt2
       drawpixel (x, y) with color texture(t) for a solid texture
       or with texture(β, γ) for a 2D texture.

This code will generate images, but there is a problem. To unravel the basic problem, let’s consider the progression from world space q to homogeneous point r to homogenized point s:

The simplest form of the texture coordinate interpolation problem is when we have texture coordinates (u, v) associated with two points, q and Q, and we need to generate texture coordinates in the image along the line between s and S. Ifthe world-space point q′ that is on the line between q and Q projects to the screen-space point s′ on the line between s and S, then the two points should have the same texture coordinates.

The naïve screen-space approach, embodied by the algorithm above, says that at the point s′ = s + α(S – s), we should use texture coordinates u_s + α(u_S – u_s) and v_s + α(v_S – v_s) . This doesn’t work correctly because the world-space point q′ that transforms to s′ is not q + α(Q – q).

However, we know from Section 8.4 that the points on the line segment between q and Q do end up somewhere on the line segment between s and S; infact, in that section we showed that

The interpolation parameters t and α are not the same, but we can compute one from the other:¹

These equations provide one possible fix to the screen-space interpolation idea. To get texture coordinates for the screen-space point s′ = s + α(S – s), compute u′_s = u_s + t(α)(u_S – u_s) and v′_s = v_s + t(α)(v_S – v_s). These are the coordinates of the point q′ that maps to s′, so this will work. However, it is slow to evaluate t(α) for each fragment, and there is a simpler way.

The key observation is that because, as we know, the perspective transform preserves lines and planes, it is safe to linearly interpolate any attributes we want across triangles, but only as long as they go through the perspective transformation along with the points (Figure 9.8). To get a geometric intuition for this, reduce the dimension so that we have homogeneous points (x_r,y_r,w_r) and a single attribute u being interpolated. The attribute u is supposed to be a linear function of x_r and y_r,soifweplot u as a height field over (x_r,y_r), the result is a plane. Now, if we think of u as a third spatial coordinate (call it u_r to emphasize that it’s treated the same as the others) and send the whole 3D homogeneous point (x_r,y_r,u_r,w_r) through the perspective transformation, the result (x_s,y_s,u_s) still generates points that lie on a plane. There will be some warping within the plane, but the plane stays flat. This means that u_s is a linear function of (x_s,y_s)—which is to say, we can compute u_s anywhere by using linear interpolation based on the coordinates (x_s,y_s).

Returning to the full problem, we need to interpolate texture coordinates (u, v) that are linear functions of the world space coordinates (x_q,y_q,z_q). After transforming the points to screen space, and adding the texture coordinates as if they were additional coordinates, we have

The practical implication of the previous paragraph is that we can go ahead and interpolate all of these quantities based on the values of (x_s,y_s)—including the value z_s, used in the z-buffer. The problem with the naïve approach is simply that we are interpolating components selected inconsistently—as long as the quantities involved are from before or all from after the perspective divide, all will be well.

The one remaining problem is that (u/w_r,v/w_r) is not directly useful for looking up texture data; we need (u, v). This explains the purpose of the extra parameter we slipped into (9.3), whose value is always 1: once we have u/w_r, v/w_r,and 1/w_r, we can easily recover (u, v) by dividing.

To verify that this is all correct, let’s check that interpolating the quantity 1/w_r in screen space indeed produces the reciprocal of the interpolated w_r in world space. To see this is true, confirm (Exercise 2):

remembering that α(t) and t are related by Equation 9.2.

This ability to interpolate 1/w_r linearly with no error in the transformed space allows us to correctly texture triangles. We can use these facts to modify our scan-conversion code for three points t_i = (x_i,y_i,z_i,w_i) that have been passed through the viewing matrices, but have not been homogenized, complete with texture coordinates t_i = (u_i,v_i):

for all xs do
  for all ys do
   compute (α, β, γ) for (xs,ys)
   if (α ∈ [0, 1] and β ∈ [0, 1] and γ ∈ [0, 1]) then
       us = α(u0/w0) + β(u1/w1) + γ(u2/w2)
       vs = α(v0/w0) + β(v1/w1) + γ(v2/w2)
       1s = α(1/w0) + β(1/w1) + γ(1/w2)
       u = us/1s
       v = vs/1s
       drawpixel (xs,ys) with color texture(u, v)

Of course, many of the expressions appearing in this pseudocode would be precomputed outside the loop for speed.

In practice, modern systems interpolate all attributes in a perspective-correct way, unless some other method is specifically requested.

9.1.4 Clipping