GameDev.netNormal Computations for Heightfield Lighting

Normal Computations for Heightfield Lighting
## Introduction
The goal of this post is to explain the most common algorithms used for computing surface and vertex normals. After a survey of the most common algorithms we'll take an in-depth look at the two most common algorithms. These in-depth looks will include a logical analysis of why the method is effective, how the math can be simplified based on what we know about heightfields, and the relative performance of the particular algorithm. By the end of the article the reader should be comfortable with computing normals as they pertain to heightfields and be familiar with which method would be best suited for the desired look and performance.
## Quick History of Normal Calculations of Polyhedra
Gouraud's algorithm would remain the defacto algorithm for computing per-vertex lighting for more than 20 years. In 1997 however, two German Computer Scientists Thόrmer & Wόthrich, proposed a new method for computing normals which was mathematically more accurate. Their algorithm, which is referred to as "Mean Weighted by Angle" (MWA), suggested that facets which are "attached" to a vertex normal should have their contribution to the vertex normal Not long after Thόrmer & Wόthrich proposed their new algorithm a man named Nelson Max published his findings in the Journal of Graphics Tools describing four additional algorithms. I will list these here for brevity, however experiments have been done which show these methods are neither as accurate as MWA or as fast as MWE. For this reason, we will not be exploring these methods further. "Mean Weighted by Sine and Edge Length Reciprocal," "Mean Weighted by Areas of Adjacent Triangles," "Mean Weighted by Edge Length Reciprocals," and "Mean Weighted by Square Root of Edge Length Reciprocals." ## Heightfields
index ═ z * numVertsWide + x The value in memory can be stored as either a floating point value or an integer. In the case of an integer, the values are often stored as a 1 byte (8 bits) or 4 byte (32 bits) block of memory. When the height values are stored as 8-bit integers the heightfield can easily be saved and loaded to disk in the form of a grayscale bitmap. This makes it possible to use a standard image editing application to modify the terrain offline. The same can be done with a 32 bit value, but the alpha value and colors become less intuitive to work with in an image editing program. When the value in memory is stored as an 8 bit integer it is also common to multiply the height values by a vertical scale component. This enables the heightfield to cover elevations of greater and lower than 255. Depending on the resolution of your terrain, however, 8 bit integers can often seem blocky and unrealistic. For the purpose of this article and demo, we will use 32 bit, single precision Floating Point values to represent height. In addition to having a vertical scale factor in the case of 8 bit heightfield, it is also common in all cases to have a horizontal scale factor. The reason for this is quite simple. When using indices for your x,z pairs it would force all of your vertices to reside along 1 meter (world units) intervals (i.e. X = 0, 1, 2, 3 .Z = 0, 1, 2, 3). As with using 8 bit values for height, forcing your vertices to reside 1 meter apart could cause the terrain to seem unrealistic in many cases. The horizontal scale factor, sometimes called "Scale" or "Units Per Vertex" allows your vertices to be spaced in distances greater or smaller than 1.0. For example, if your scale factor is 0.5f, then all of your vertices would be 0.5 meters apart, leading to more vertices in the same amount of space, and causing your terrain to be a higher "Resolution." It should be noted here that "Scale" and "Resolution" are inversely proportional. As your scale increases, causing a greater distance between your vertices, your resolution, or the complexity of the terrain, decreases. Finally, although we represent a heightfield as a single array of floating point values in this article, it is common for culling and rendering purposes to break a large heightfield up into smaller "tiles." If this is what you're trying to do, you can think of the heightfield in this article as being a single tile of the larger system. Ultimately, tiling your heightfield means creating a mapping between tiles and vertices, and also extending your normal calculations to take into account the heights of vertices in tiles adjacent to the current tile. You will see later that our algorithm for computing surface normals already takes tiling into account, and with little modification you can get these algorithms to work within a tiled terrain system. ## Normal Calculations for Heightfields## Method 1: Mean Weighted Equally (MWE)Mean Weighted Equally is a two stage algorithm for computing the normal at a vertex in a heightfield (vertex normal). The first stage in the process is to compute the surface normals of each of the triangles in the heightfield using cross-product calculations. The computed surface normals can either be normalized or not, depending on the quality of vertex normals you are looking for. The second stage of the process is to add up all of the surface normals attached to a given vertex and then either average the resultant vertex normal or normalize it. The following two equations are given as the mathematical definitions for the surface and vertex normals. Surface Normal: B
Vertex Normal: NΣ (_{nwe} = N_{si}) / 6.0f | i = 0
5The equation for the surface normal above can be read as "The Surface Normals is equal to the vector The equation for the vertex normal above can be read as "The vertex normal is the sum of all surface normals which this vertex is a part of, divided by 6, where the surface normal index is in the range of 0 to 5, inclusive. I should take a moment to note that in Gouraud's algorithm he does not strictly state that there are up to 6 surfaces connected to a single vertex, however this is a mathematical certainty for a heightfield where all of the triangles are oriented the same direction. As well, based on the algorithm we use to implement the above equation each vertex will have
In the diagram to the right each quad is made up of two triangles, one in the lower left of the quad and one in the top right. For the sake of this article you can assume that we're using DirectX's default left-handed coordinate system and that triangles are oriented positively in a clock-wise ordering. You can also assume that our heightfield lies in the xz-plane, and that y is the coordinate of height. From the diagram you can see that quad 0 is made up of triangles 0-3-1 and 1-3-4. Mathematically speaking a plane or surface can be computed from three points using the cross-product of the vectors formed by two of the edges of the triangle. Equation for a Normal Vector from a Left-Handed Triangle:
As an example, given triangles 031, 134, and height values as a function of (x,z) pairs you can compute the normals to the triangles with the following equations: // Triangle 0 Normal
_{3}, 1 ) ( 0, h_{0}, 0 ) ] Χ [( 1, h_{1}, 0 ) ( 0, h_{0}, 0 ) ]
Optionally normalize N
_{0}// Triangle 1 Normal A = 3 1; A = ( 0, h_{3}, 1 ) ( 1, h_{1}, 0 )
B = 4 1; B = ( 1, h_{4}, 1 ) ( 1, h_{1}, 0 )
N[ ( 0, h_{1} = _{3}, 1 ) ( 1, h_{1}, 0 ) ] Χ [ ( 1, h_{4}, 1 ) ( 1, h_{1}, 0 ) ]
Optionally normalize N_{1}The above equations illustrate how surface normals can be computed for a single quad of the heightfield, however, this is a special case and we need a general form of the equation if we're going to apply it to the terrain as a whole. In general, we'd like an equation that uses i, and j notation, instead of specific x,z pairs to allow us to iterate over the entire terrain and maintain the same equation. So given that i is an index in the range of 0 to for( j = 0; j < numVertsDeep; j++ ) for( i = 0; i < numVertsWide; i++ ) { ... } For the triangle in the bottom left, hereafter referred to as Triangle 0, let P be the point in the top left and _{1}P be the point in the bottom right. If _{2}P, _{0}P, and _{1}P were a function of the current [i, j] pair, the points could be represented as follows:_{2}
_{0}, j ]
P = [ i, h_{1}_{1}, j+1 ]
P = [ i+1, h_{2}_{2}, j ]// Triangle 0
P
_{0}B = P _{2}P_{0}Now that we can represent any point on the heightfield as a function of i' and j', let's re-examine the equation for the normal. Even though we are rephrasing the equations with (i, j) coordinates, the pairs will drop out leaving us with simpler calculations at the end. // Compute the two vectors used in the cross product
// Compute the normal to Triangle 0, using the generic form i,j and the full definition of cross-product
A Χ B
N = [ A_{0}_{y}B_{z} − A_{z}B_{y}, A_{z}B_{x} A_{x}B_{z}, A_{x}B_{y} A_{y}B_{x} ]
_{1 }- h_{0}, 1 ) ] Χ ( 1, h_{2} - h_{0}, 0 ) ]
N[ 0( h_{0} = _{1 }- h_{0 }) 1( h_{2} - h_{0 }), 1(1) 0(0), 0( h_{1 }- h_{0 }) 1( h_{1 }- h_{0 }) ]
N[ h_{0} = _{0 }- h_{2}, 1.0f, h_{0} - h_{1} ]In the equation above the cross product of two orthogonal vectors can be simplified dramatically. The calculation of a cross product normally requires 6 multiplications and 3 subtractions, however once we've computed it to the end we see that many of the terms drop out and that only 2 subtractions are required. Before you reach for the D3DXVec3Cross() function remember that DirectX knows nothing about the nature of your vectors. In this case, based on the orthogonal nature of our vectors it is far faster to compute the cross product ourselves than to let DirectX do it for us. Now that we've computed the cross product of Triangle 0 in general form, let's quickly compute the cross product of Triangle 1 using generic form. Once we're done we'll have the two equations necessary to iterate over the entire heightfield and compute surface normals. // The 3 points used to compute the top right triangle
P = [ i, j+1 ]
_{1}P = [ i+1, j+1 ]_{2}// Compute the two vectors used in the cross product
P; _{0}A = ( i, h_{1}, j+1 ) ( i+1, h_{0}, j ); A = ( -1, h_{1 }- h_{0}, 1 )
B = P _{2}P; _{0}B = ( i+1, h_{2}, j+1 ) ( i+1, h_{0}, j ); B = ( 0, h_{2} - h_{0}, 1 )
_{y}B_{z} − A_{z}B_{y}, A_{z}B_{x} A_{x}B_{z}, A_{x}B_{y} A_{y}B_{x} ]
N[( -1, h_{0} = _{1 }- h_{0}, 1 ) ] Χ ( 0, h_{2} - h_{0}, 1 ) ]
N[ 1( h_{0} = _{1 }- h_{0 }) 1( h_{2} - h_{0 }), 1(0) (-1)(1), (-1)( h_{2} - h_{0 }) (0)( h_{1 }- h_{0}) ]
N[ h_{0} = _{1 }- h_{2}, 1.0f, h_{0 }- h_{2} ]As with the equation for Triangle 0, Triangle 1 also simplifies down to 2 subtractions as one might expect. You are now armed with everything you need to compute the surface normals for every triangle on the heightfield. Simply iterate over the quads, and for each quad compute the two surface normals using the equations above. Before we look at the code for computing the surface normals, lets quickly evaluate the range over which we will iterate. A heightfield will have voidHeightfield::ComputeSurfaceNormals( void ) { INT normalIndex = 0; // Iterate through the z components for( int z = -1; z <= m_NumVertsDeep - 1; z++ ) { // Iterate through the x components for( int x = -1; x <= m_NumVertsWide - 1; x++ ) { // Does the current x,z lie on the terrain, or just outside? if( !IsValidCoord( x, z ) || !IsValidCoord( x+1, z ) || !IsValidCoord( x, z+1 ) || !IsValidCoord( x+1, z+1 ) ) { normalIndex += 2; continue; } // Set/Compute the normal for triangle 0 float height0 = GetHeight( x, z ); float height1 = GetHeight( x, z+1 ); float height2 = GetHeight( x+1, z ); D3DXVECTOR3 normal1( height0 - height2, 1.0f, height0 - height1 ); // Set/Compute the normal for triangle 1 height0 = height2; height2 = GetHeight( x+1, z+1 ); D3DXVECTOR3 normal2( height1 - height2, 1.0f, height0 - height2 ); // Set the surface normals in the array m_pSurfaceNormals[normalIndex++] = normal1; m_pSurfaceNormals[normalIndex++] = normal2; } } }
Vertex Normal: N_{si}) / 6.0f | i = 0
5In order to determine which of the surfaces are adjacent to the current vertex let's once again take a look at the image to the right. In the drawing to the right each of the vertices has a different number of surrounding surfaces. Vertex 0 for example has only one surface adjacent to it. Vertex 1 has three surfaces adjacent, two in the top left, and one in the top right. Vertex 2 has two adjacent surfaces to its upper left. And vertex 4 has six adjacent surfaces three above it, and three below it. In order to create an algorithm in this example that accurately computes the vertex normal of every vertex we would have to take into account each of the special cases. Is the vertex on an edge top, left, bottom, right? Is it on a corner top left, top right, bottom left, or bottom right? Such an algorithm would first have to compute conditional statements to determine which of the above statements is true. Once the condition of the vertex was known it would likely progress through a series of if-statements where certain adjacent surfaces would be summed. There's got to be a cleaner way! The primary reason we have this problem is because the heightfield has edges. If the heightfield had no edges then all vertices would have 6 surface normals such as the middle vertex to the image to the right. In which case, it might be possible to create an algorithm that reliably determined the six surrounding surfaces for any given vertex. So let's remove the edges from the heightfield! If you recall from the section above on computing the surface normals I illustrated how you should add an extra row and column of quads to each edge. The quads would not be updated but would merely have the UP vector as its normal. The reason we did that is so we could change the previous image to the one just below. In the image to the right we no longer have edges to our heightfield because of the additional sets of quads. Every vertex on our heightfield now has 6 surrounding surfaces. Based on the direction of our triangles, each vertex has two surfaces to the top left, one to the top right, one to the bottom left, and two to the bottom right. Now imagine there was a small square around any given vertex. The square would be a part of 6 surfaces. If you slide the square to the left, right, up, or down, it is always located in six triangles. Therefore, we will call this algorithm The idea behind the algorithm is to put the rectangle in the bottom left corner, determine which surfaces it's part of, and then add those surfaces to a list to be summed for the current vertex index. Let's try it with a few vertices to help you get the hang of it. Vertex 0
Vertex 1
Vertex 2
I wanted to get three vertices above because I wanted you to see the pattern evolve. Each row of triangles is sequential from one vertex to the next. Vertex 0 had surfaces 1, 2, and 3, where vertex 1 had surfaces 3, 4, 5 and so on. This would have continued all the way across the heightfield. The same was true for the top row. So if you could get a starting index for each vertex, we'll call it the triangleIndex, then we'd just need to start with that index and move along two to the right, and the corresponding three to the top and we'd be done. The only question would be how to get the surface which is immediately above the triangle index. That can be obtained from the following equation: rowOffset = 2( numVertsWide + 1 ) - 1; Here's the algorithm as implemented in the attached demo: void Heightfield::ComputeVertexNormals( void ) { D3DXVECTOR3 vec[6]; D3DXVECTOR3 vertexNormal; INT triIndex = 1; INT indexPlusOne = 2; INT indexPlusTwo = 3; INT vertIndex = 0; // Computes the offset based on the number of triangles across INT rowOffset = ( ( m_NumVertsWide + 1 ) * 2 ) - 1; for( INT z = 0; z < m_NumVertsDeep; z++ ) { for( INT x = 0; x < m_NumVertsWide; x++ ) { indexPlusOne = triIndex + 1; indexPlusTwo = triIndex + 2; // Get the three triangles below the vertex vec[0] = m_pSurfaceNormals[ triIndex ]; vec[1] = m_pSurfaceNormals[ indexPlusOne ]; vec[2] = m_pSurfaceNormals[ indexPlusTwo ]; // Get the three triangles above the vertex vec[3] = m_pSurfaceNormals[ rowOffset + triIndex ]; vec[4] = m_pSurfaceNormals[ rowOffset + indexPlusOne ]; vec[5] = m_pSurfaceNormals[ rowOffset + indexPlusTwo ]; // Sum the vectors and then divide by 6 to average them vertexNormal = vec[0] + vec[1] + vec[2] + vec[3] + vec[4] + vec[5]; vertexNormal /= 6.0f; m_pVertexNormals[vertIndex] = vertexNormal; triIndex += 2; vertIndex++; } triIndex += 2; } } ## Method 2: Mean Weighted by AngleAs with the algorithm above, Mean Weighted by Angle (MWA) is a two stage algorithm. The first stage is to compute the surface normals, much as in the algorithm above, except that each surface normal is post-processed before it is used in the computation of vertex normals. Unlike in MWE, the surface normals MUST be normalized or artifacts may appear in the lighting. The second stage of the process is mostly the same as it was in MWE, sum up the associated surface normals and then divide them by 6. Here is the modified formula used in computing the weighted surface normal: Surface Normal: (A ΧB)Since computation of the vertex normal is
Definition of Cross-Product:
Θ ═ sin The above equation, which is generally used as the definition of a cross-product can be re-written as the line just below it. The line below it tells us that the angle between vectors │ V_{x })^{2} + ( ∆V_{Y })^{2} + ( ∆V_{Z })^{2} )
│ V│ = √( ( i i )_{0}^{2} + ( h_{1}-h_{0 })^{2} + ( j+1-j_{ })^{2} )
│ V│ = √( ( 0_{0}^{2} + ( h_{1}-h_{0 })^{2} + ( 1_{ })^{2} )
│ V│ = √( ( h_{0}_{1}-h_{0 })^{2} + 1 )Through simplification we have reduced the magnitude of the vector down to a single multiply, addition, square, and square root. This is NOT an inexpensive equation. The square kills the performance of this algorithm as you will see later. Although we've computed the magnitude for the left edge of the triangle, there are still two more edges, each which have their own magnitudes. Below are the simplified equations for the magnitude of the two other edges of Triangle 0. │ _{2}-h_{1 })^{2} + 2 ) // Hypotenuse of Triangle 0
│ V│ = √( ( h_{2}_{2}-h_{0 })^{2} + 1 ) // Bottom edge of Triangle 0As with the equations for normal calculations, the orientation of Triangle 1 is slightly different and as such, has slightly different magnitude equations as well. Here are the equations for the magnitude of the three edges of Triangle 1. │ _{1}-h_{0 })^{2} + 2 ) // Hypotenuse of Triangle 1
│ V│ = √( ( h_{1}_{2}-h_{1 })^{2} + 1 ) // Top edge of Triangle 1
│ V│ = √( ( h_{2}_{2}-h_{0 })^{2} + 1 ) // Right edge of Triangle 1With the six equations above you have the ability to compute the magnitude of each of the edges of both the bottom left and top right triangles of every quad on the heightfield. The next thing to do is to plug them into the equation for Theta(angle) above, and compute the angle of each corner of the triangles. We will do Triangle 0 here for demonstration, Triangle 1 will be left up to the reader to prove, however the source code will be provided. Θ ═ sin // Compute the angle in the bottom left corner of Triangle 0
_{1}-h_{0 })^{2} + 1 ) // Left edge of Triangle 0
│ V│ = √( ( h_{2}_{2}-h_{0 })^{2} + 1 ) // Bottom edge of Triangle 0Θ ═ sin V││_{0} V│ ) where Ns is the already computed normal to the surface of the triangle_{2}To compute the weighted surface normal for the vertex in the bottom left corner of Triangle 0 do the following:
Listed here is the modified code for computing surface normals which iterates over every quad on the heightfield and computes the weighted surface normals for Triangle 0 and Triangle 1. void Heightfield::ComputeWeightedSurfaceNormals( void ) { INT normalIndex = 0; // Iterate through the z components for( int z = -1; z <= m_NumVertsDeep - 1; z++ ) { // Iterate through the x components for( int x = -1; x <= m_NumVertsWide - 1; x++ ) { if( !IsValidCoord( x, z ) || !IsValidCoord( x+1, z ) || !IsValidCoord( x, z+1 ) || !IsValidCoord( x+1, z+1 ) ) { normalIndex += 6; continue; } // Set/Compute the normal for triangle 1 float height0 = GetHeight( x, z ); float height1 = GetHeight( x, z+1 ); float height2 = GetHeight( x+1, z ); D3DXVECTOR3 normal1( height0 - height2, 1.0f, height0 - height1 ); D3DXVec3Normalize( &normal1, &normal1 ); // Get the Magnitude of the normal to the surface and the weights float normal1Mag = D3DXVec3Length( &normal1 ); float abMag = sqrt( (height1 - height0) * (height1 - height0) + 1 ); float acMag = sqrt( (height2 - height0) * (height2 - height0) + 1 ); float bcMag = sqrt( (height2 - height1) * (height2 - height1) + 2 ); float theta0 = asin( normal1Mag / ( abMag * acMag ) ); float theta1 = asin( normal1Mag / ( bcMag * abMag ) ); float theta2 = asin( normal1Mag / ( bcMag * acMag ) ); m_pWeightedSurfaceNormals[normalIndex++] = normal1 * theta0; m_pWeightedSurfaceNormals[normalIndex++] = normal1 * theta1; m_pWeightedSurfaceNormals[normalIndex++] = normal1 * theta2; // Set/Compute the normal for triangle 2 height0 = height2; height2 = GetHeight( x+1, z+1 ); D3DXVECTOR3 normal2( height1 - height2, 1.0f, height0 - height2 ); D3DXVec3Normalize( &normal2, &normal2 ); // Get the Magnitude of the normal to the surface and the weights float normal2Mag = D3DXVec3Length( &normal2 ); abMag = sqrt( (height1 - height0) * (height1 - height0) + 2 ); acMag = sqrt( (height2 - height0) * (height2 - height0) + 1 ); bcMag = sqrt( (height2 - height1) * (height2 - height1) + 1 ); theta0 = asin( normal2Mag / ( abMag * acMag ) ); theta1 = asin( normal2Mag / ( bcMag * abMag ) ); theta2 = asin( normal2Mag / ( bcMag * acMag ) ); m_pWeightedSurfaceNormals[normalIndex++] = normal2 * theta0; m_pWeightedSurfaceNormals[normalIndex++] = normal2 * theta1; m_pWeightedSurfaceNormals[normalIndex++] = normal2 * theta2; } } // end surface normal calculations } ## Analysis of the performance of each approachDuring the course of writing this article I decided to create a demo which illustrated the different topics in action. As I was researching the best way to implement the various heightfield algorithms I came upon some contradictions. Some people said you should normalize the surface normals before using them in any vertex normal calculation, some said it didn't matter. Also, many said the best way to compute the vertex normal was to average the surface normals, and others said just add up the surface normals and then normalize the resulting vector. Finally, there's the question of which algorithm looks best and performs best between Mean Weighed Equally and Mean Weighted by Angle. This section will attempt to compare and contrast the differences between the various methods. It should be understood that the performance analysis I give below is based upon only have a single heightfield in the scene, with no other objects, which is animating on a Sine curve every frame. Actual performance in real-time situations could be both better and worse than those given below based upon whether your terrain is static or dynamic, and what else you have going on in your scenes. The information provided below is strictly for the purpose of comparing computational complexity of the different algorithms.
## Future research / other areas of interestIn addition to the topics I have covered in this article, I was asked to cover two additional areas of interest. Unfortunately, I did not have time to fully explore these topics. I list them here so the reader is aware of them and can do their own research if interested.
In the case of heightfields, smoothing groups could be useful at the peaks of tall mountains or at the edges of cliffs and drop offs. Either of the algorithms above would cause a smooth curve down the steep drop. In reality, there are often sudden drops in which there is not a smooth curve. Smoothing groups can then be used to give your terrain more realistic edges.
The application of this on a heightfield is less obvious. I can envision this only being useful when your terrain is procedurally generated from a well known function. In the demo I provided the positions of all the vertices are computed from a Sine function. In this example, it should be possible, knowing that the height of my vertices range from 0 to 1 as a function of 0 to 360 degrees, to compute the normal at each of the vertices based on the angle between 0 and 360 of that particular vertex. This will be left up to the reader to implement. ## ConclusionDuring the past two weeks I have had the opportunity to write a demo which explores the implementation of lighting normals of a heightfield while at the same time being able to make empirical observations about the usefulness and performance of each of the different popular algorithms. This article contains all that I have discovered about computing normals for lighting of a heightfield and I hope that it has been both informative and helpful. Based on my findings I will implement my heightfields using normalized, weighted surfaces, and averaged - not normalized, vertex normals whenever possible. For those that are part of the GameDev.net community I hope to see you on the forums and I look forward to your feedback regarding this article. Feel free to send your code corrections and/or comments directly to me via PM. ## Appendix A: Demo Description, Links, and Key bindingsThe screen shots, performance comparisons, code snippets, etc. of this article were all based upon information gathered from the attached demo. The demo was an opportunity for me to test various types of functionality and to determine the visual and performance impacts of those features. Because of that, the demo has a limited amount of functionality.
## References- "A Comparison of Algorithms for Vertex Normal Computation", Shuangshuang Jin, Robert R. Lewis, David West, 2003, Washington State University
© 1999-2011 Gamedev.net. All rights reserved. |