Method 2: Mean Weighted by Angle

As 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: N_s ═ Θ(A ×B)

Since computation of the vertex normal is mostly the same in this algorithm as it is in MWE, we will not cover it again. For specific implementation details and to see how it was slightly modified to access the weighted surface normals see the code in the attached demo.

Computation of surface normals
The first part of the algorithm for MWA is to compute the surface normals, exactly as we did with MWE. After we have computed the surface normal for a given triangle, however, we must post-process it by multiplying it with all three angles of the surface triangle. In doing so we will have created three new "weighted" surface normals that will replace the single surface normal used before. We already know how to get the surface normal so let's explore how to get the angle. The following equation is a commonly used method for getting the angle between two vectors:

Definition of Cross-Product:
│A × B│ = │A││B│sin(Θ)

Θ ═ sin^-1( │A × B│ / │A││B│ )

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 A and B can be determined by taking the Inverse Sine of the Magnitude of the Cross-Product between A and B, and dividing all that by the magnitude of vector A multiplied by the magnitude of vector B. As mentioned above, we already have the cross-product of a specific A and B for the given triangle, however we also need the magnitude of every edge of the triangle. The next step is to determine the most efficient way to compute the magnitude of a vector. The following equations derive the general form of computing the magnitude for the vector which lies along the left edge of Triangle 0.

│V₀│ = √( ( ∆V_x )² + ( ∆V_Y )² + ( ∆V_Z )² )
│V₀│ = √( ( i – i )² + ( h₁-h₀ )² + ( j+1-j )² )
│V₀│ = √( ( 0² + ( h₁-h₀ )² + ( 1 )² )
│V₀│ = √( ( h₁-h₀ )² + 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.

│V₁│ = √( ( h₂-h₁ )² + 2 ) // Hypotenuse of Triangle 0
│V₂│ = √( ( h₂-h₀ )² + 1 ) // Bottom edge of Triangle 0

As 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.

│V₀│ = √( ( h₁-h₀ )² + 2 ) // Hypotenuse of Triangle 1
│V₁│ = √( ( h₂-h₁ )² + 1 ) // Top edge of Triangle 1
│V₂│ = √( ( h₂-h₀ )² + 1 ) // Right edge of Triangle 1

With 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^-1( │A × B│ / │A││B│ )

// Compute the angle in the bottom left corner of Triangle 0
│V₀│ = √( ( h₁-h₀ )² + 1 ) // Left edge of Triangle 0
│V₂│ = √( ( h₂-h₀ )² + 1 ) // Bottom edge of Triangle 0

Θ ═ sin^-1( │N_s│ / │ V₀││ V₂│ ) where Ns is the already computed normal to the surface of the triangle

To compute the weighted surface normal for the vertex in the bottom left corner of Triangle 0 do the following:

N_w0 = Θ * N_s

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 approach