Library/Graphics/Bump Mapping

Preface
This article describes one approach to bump mapping in GLSL - and then dives into a detailed explanation of the math behind the code.

It is important to note that this is only one approach. There are methods that are more efficient and produce better visual results. However, the approach described here should provide a good bit of detail behind bump mapping so that, from here, a reader can understand the math better and choose the method that works best for their needs.

Some advantages and disadvantages of the approach described below are:  The GLSL code is general, concise, and easily encapsulated into an independent function It does not require explicit tangent and binormals to be specified with the geometry It works with any height mapping function or texture (not necessarily only 2d textures) It is relatively computationally expensive on the GPU: much of which could be done higher in the pipeline There are other approaches that produce smoother visual results 

The GLSL Fragment shader
Let's start with some completely uncommented code and then provide the explanation.

Assume the input variable "fragVertex" is the fragment position and "fragNormal" is the fragment normal. The coordinate system in which these are expressed does not matter as long as they are consistent. The "value" parameter should be the height value as sampled for that fragment.

The above is the code to cut & paste if you want to simply want working code without bothering with why it works. Pass in the fragment position and normal in the same coordinate system and the height value at the fragment, then use the return value as the new normal. It should "just work."

Assuming however that you'd like to know why it works, let's explain...

Bump Mapping Perturbs Normals
The basic goal of bump mapping is to apply a height map to the surface of an object. This height map does not actually change the geometry of the object, but only perturbs the normal at that point of the surface as if the geometry were changed. This is computationally less expensive than actually changing the geometry but still yields a good visual effect.

Therefore, the goal is to compute the normal of the height map at the sample point, transform that into the coordinate system of the object, and use that normal as the surface normal in place of the "actual" geometric surface normal.

The first step is to find the normal at a given point on the height map...

Find the Normal via Discrete Samples
I'm going to rush through this section assuming the reader is okay with the math being here. Leave a comment if more detail is desired.

If you imagine the height map as texture, the normal at a point can be computed creating a triangle at that point connecting it to it's neighbors and computing the normal of that triangle.

In three-space, define three points where $$F(x,y)$$ is defined as the height map at $(x,y)$:


 * $$P_0 = < x,y,F(x,y) > $$
 * $$P_1 = < x+1,y,F(x+1,y) > $$
 * $$P_2 = < x,y+1,F(x,y+1)> $$

The normal of the surface is defined by the cross product of two perpendicular vectors lying along surface, thus:

$$N = normalize( (P_1 - P_0) \times (P_2 - P_0)) $$

(I'm assuming this is all pretty straightforward to the reader at this point.)

Another common approach would be to sample the neighbors in both directions, with the equation thus re-defined as:


 * $$U = < x + 1, y, F(x+1,y) > - < x -1, y, F(x-1,y) >$$
 * $$V = < x, y+1, F(x,y+1) > - < x, y-1, F(x,y-1) >$$
 * $$N = normalize(U \times V) $$

Find the Normal via a Derivative
Let's now think of the height map as a continuous function $$z = F(x,y)$$ rather than necessarily a discrete texture. We'll work through the steps to restate the above discrete sampling in terms of a derivatives of $F$.

Let's go back to this definition of finding the height map normal:


 * $$P_0 = < x,y,F(x,y) > $$
 * $$P_1 = < x+1,y,F(x+1,y) > $$
 * $$P_2 = < x,y+1,F(x,y+1)> $$
 * $$N = normalize( (P_1 - P_0) \times (P_2 - P_0)) $$

Step 1: Since we have a continuous function, let's not move "1" unit in x and y, but rather a small scalar delta $dx$ and $dy$.


 * $$P_0 = < x,y,F(x,y) > $$
 * $$P_1 = < x+dx,y,F(x+dx,y) > $$
 * $$P_2 = < x,y+dy,F(x,y+dy)> $$
 * $$N = normalize( (P_1 - P_0) \times (P_2 - P_0)) $$

Step 2: Substitute in for $P$ and simplify:

$N = normalize( (< x+dx,y,F(x+dx,y) > - < x,y,F(x,y) >) \times (< x,y+dy,F(x,y+dy)> - < x,y,F(x,y) >)) $

The equation is a bit long, but simple:

$$N = normalize( < dx,0,F(x+dx,y) - F(x,y)> \times < 0, dy, F(x,y+dy) - F(x,y) > )$$

Step 3: Reform the equation a bit by multiplying the last components by $dx/dx$ and $dy/dy$, respectively:


 * $$N = normalize( < dx,0, \frac{F(x+dx,y) - F(x,y)}{dx}dx> \times < 0, dy, \frac{F(x,y+dy) - F(x,y)}{dy}dy > )$$

We do this to use the fundamental theorem of calculus to state the equation in terms of derivatives. The fundamental theorem of calculus states:


 * $$\frac{dF}{dx}= \lim_{dx\rightarrow 0} \frac{F(x + dx) - F(x)))}{dx} $$

Thus our previous equation can be restated as:


 * $$N = normalize( < dx,0, \frac{dF}{dx}dx> \times < 0, dy, \frac{dF}{dy}dy > )$$

Thus we have the normal of a arbitrary height map function $F$ stated in terms of the partial derivatives of $F$.

Step 4: Simplify the math a bit further:

Multiplying the cross product out yields further simplifications:


 * $N = normalize( < dx,0, \frac{dF}{dx}dx> \times < 0, dy, \frac{dF}{dy}dy > )$
 * $N = normalize( < -dy \frac{dF}{dx}dx, -dx \frac{dF}{dy}dy, dx dy >  )$

Since the vector is being normalized, we can safely scale each component by the same factor $\frac{1}{dx dy}$


 * $N = normalize( \frac{1}{dx dy} < -dy \frac{dF}{dx}dx, -dx \frac{dF}{dy}dy, dx dy >  )$
 * $N = normalize( < -\frac{dF}{dx}, -\frac{dF}{dy}, 1 > )$

Voila! We have a simple equation stating the height map normal in terms of partial derivatives:

$$N = normalize( < -\frac{dF}{dx}, -\frac{dF}{dy}, 1 > )$$

Transforming the Normal to the Correct Coordinate System
The equation:

$$N = normalize( < -\frac{dF}{dx}, -\frac{dF}{dy}, 1 > )$$

yields the normal in terms of the x and y of the flat planar surface of the 2D function $F(x,y)$. Stated differently, z=1 means "up" in the coordinate system of $F$, but we need to transform this to the surface of the object so that "up" is in the direction of the surface's normal prior to any bump mapping.

We'll do this transformation via derivatives as well.

Transformation via Derivatives
Let's define an arbitrary space `A`, whose axes we will label $x,y,z$. It could represent eye-coordinates, object coordinates, or any other orthonormal coordinate system. Let's define a second coordinate system $B$ with axes $u,v,w$.

Let's now see how we can express the axes of u,v,w in terms of $A$. Let's take an identify function $G(x,y,z)$ which takes coordinates from space $A$ and returns them exactly:

$G(< x,y,z >) = < x,y,z > $

Now, let's take the derivative with respect to $u$:

$\frac{dG}{du} = \frac{d < x, y, z >}{du} $

$\frac{dG}{du} = < \frac{dx}{du}, \frac{dy}{du}, \frac{dz}{du} > $

In other words, the derivative of the identity function in $A$ by an axis in $B$ is equal to that axis of $B$ expressed in terms of the coordinate system of $A$.

So! That sounds a bit complicated, but a concrete example should clarify what this implies:

If we have a function $P$ that is equal to a position in, say, eye-space and a set of function that computes derivatives in another space, say window coordinates (which we'll call $\frac{d}{du}$, $\frac{d}{dv}$, and $\frac{d}{dw}$), then the axes of that second coordinate system can be expressed in terms of the first via the following:


 * $U = \frac{dP}{du}$
 * $V = \frac{dP}{dv}$
 * $W = \frac{dP}{dw}$

(The mapping is quite straightforward when you think of the above as, "How does the position in eye-coordinates change as we along the window coordinate axis u?" Then ask the same for v and w.)

...and thus transformations can be expressed as:

$N' = N_u U + N_v V + N_w W$

$N' = N_u \frac{dP}{du} + N_v \frac{dP}{dv} + N_w \frac{dP}{dw}$

Transformation of N from u,v,w to x,y,z space thus can be defined as:

$$N_{x,y,z} = N_u \frac{dP}{du} + N_v \frac{dP}{dv} + N_w \frac{dP}{dw}$$

Back to the GLSL Fragment Shader
That was a lot of math. Let's get back to the simple half dozen lines of GLSL code to implement bump mapping.

Step 1: Compute the height map normal in window coordinates

The first step of the shader is to compute the partial derivatives of the height value. These values of partial derivatives are with respect to window coordinates. GLSL computes the partial derivative approximations by looking at the values in neighboring fragments and computing the deltas; thus the coordinate system of these derivatives is window coordinates.

This gives us enough information to compute the height map normal in window coordinates:

$$N = normalize( < -\frac{dF}{dx}, -\frac{dF}{dy}, 1 > )$$

Step 2: Compute the window coordinate axes in the desired coordinate system.

We now need the normal in the same coordinate system as the fragment shader is primarily operating in.

The fragment position in that coordinate system acts as an identify function for that coordinate system. Thus, the partial derivatives of that function in window coordinates yields the window coordinate axes expressed in the coordinate system of the fragment position !

But, of course, there's no  GLSL function. However, we're trying to create a orthogonal basis, therefore dPdz must be perpendicular to the other two axes and thus can be found by a cross product:

Step 3: Transform the normal

From the equations above:

$$N_{u,v,w} = normalize( < -\frac{dF}{du}, -\frac{dF}{dv}, 1 > )$$

and

$$N_{x,y,z} = N_u \frac{dP}{du} + N_v \frac{dP}{dv} + N_w \frac{dP}{dw}$$

therefore:

$$N_{x,y,z} = -\frac{dF}{du} \frac{dP}{du} -\frac{dF}{dv} \frac{dP}{dv} + \frac{dP}{dw}$$

and written as GLSL using the variables defined in the prior steps:

Compute the height map normal and then transform it using partial derivatives.

We're done! We have the height map normal from an arbitrary function $F$ re-expressed in the coordinate system of the vertex position!

But wait, there's one more thing...

Account for Smooth Shading
The above math is completely correct given the surface normal as determined by the vertex position. However, the surface normal usually isn't determined solely by the positions of the vertices.

Remember smooth versus flat shading? Smooth shading is already layering a normal perturbation "approximation" of more complex geometry than is really there. By using $dPdx$ and other associated variables, the height map normal is being applied relative to the flat surface of the triangle. Therefore, where there is no indentation or bump (i.e. the height map function is constant), the result will be flat shading. This is not what is desired.

Therefore, we actually have to transform the height map normal onto the space as defined by the "pseudo-surface" defined by the smooth normal. Given all the math we've already been through, we're going to throw a solution out there without much of an explanation. There's still this one chunk of code from the middle of the function that hasn't been mentioned yet:

Essentially all this is doing is "recreating" the window coordinate to fragment coordinate basis such that it's rotated such that Z=1 maps to fragNormal rather than the implied "flat" normal of the triangle. Since the smooth and flat normally should be relatively similar this technique works correctly. Maybe someday I'll post on how to improve that last little hack to be more mathematically accurate, but for now, this post is already way too long!

There you have it: bump mapping!

How Does Tangent Space Fit In?
Tangent space is defined as the coordinate system in which the height map's normal is defined. It is not the same as the coordinate system as the object.

In the GLSL code provided above, tangent space is implicitly window coordinates. By computing the height map normal via the  and   functions, the resulting normal is in window coordinates. Therefore, the "tangent space" being used here is window coordinates.

More "traditionally", bump maps (and normal maps) are stored in a 2D texture. As such, the computation for the normals in the map are usually done - not with  and   - but rather by sampling (as discussed at the beginning of this article). The sampling is occuring in the uv-space of the texture. Therefore, the transformation of the resulting normal when applying to the object must apply the texture UV space to the object's surface in object or eye coordinates (whichever is being used in the shader).

In other words, the transformation needs to know how the U and V texture coordinates change relative to the desired coordinate system. Most discussions about tangent space revolve around the computation of the texture space U and V axes in terms of the surface coordinate system. These U and V axes are given the special names of the tangent vector and binormal (or bi-tangent) vectors.

Using the equations above, if the normal is computed in the u,v,w space of the texture (rather than window coordinates), $dP/du$ would need to be computed where $u$ represents a unit in the texture UV space - not window coordinates. GLSL doesn't provide such a function!

In any case, the computation of tangent space is covered in many, many other places in books and on the web; therefore, I leave it to those sources to explain.