Fast Normal Calculations for Planets

Up till now I've been using a brute force method of normal calculation for the patch normals. I would use the cross product to find the normal for each triangle and then smooth the triangle normals when calculation the per-vertex normals.

I was aware of the simplified normal calculations for flat heightfields but wasn't sure how to translate those to the sphere - until now. Aurelio Reis and I had been corresponding and I emailed him about how I didn't think his method for mapping normals to the sphere was correct (turns out I completely misunderstood his method) and I recommended translating the normal on the cube face to a normal on the sphere by getting the rotation matrix from the normal of the cube face to the vertex position on the unit sphere.

The fast method for deriving a rotation matrix given a start and an end vector can be found here: http://www.cs.brown.edu/~jfh/papers/Moller-EBA-1999/main.htm - Tomas Moller & John F. Hughes "Efficiently Building a Matrix to Rotate One Vector to Another and their sample code here: http://jgt.akpeters.com/papers/MollerHughes99/code.html

The funny thing is my implementation failed and I had given up, when Aurelio wrote me back with news of his success, and some helpful advice about calculating the heightfield normals. With that encouragement I attacked it again, fixed my bugs and got it working.

Thanks to Aurelio's suggestion, I'm using the Game Programming Gems 3 method for calculating heightfield normals. The tricky thing with this method of calculating normals is that it assumes the distance between the vertices is 1 unit. When the camera approaches the terrain and it splits the quadtree, the new patches are twice the resolution of the old patch, so the vertices are now 0.5 units apart. Thankfully the distance of the normal vector doesn't matter in the calculations so it simplifies things and you can see below that for level N of the quadtree to calculate a normal for a vertex at that level the equation is:

Normal = ( h3 - h1, 2.0 / (1 << level), h4 - h2 )

Where h1-h4 are the heights of the neighboring vertices, and level is the depth in the quadtree (starting at depth 0).

Here's what that equations looks like for levels 0, 1, and 2:

Level 0:
Nv = ( h3 - h1, 2, h4 - h2 )

Level 1:
Nv = ( h3 - h1, 1, h4 - h2 )

Level 2:
Nv = ( h3 - h1, 0.5 , h4 - h2 )

I normalise the normals before applying the rotation matrix that rotates them from cube space to sphere space.

Here's the basic code for calculating the normals (hMap is a RGBA image where the RGB values are the XYZ for the vertex position on the unit sphere and the A value is the height in the heightmap)

    // Get the Y factor for the normal calculations
    Real f = 2.0  / ((Real)(1 << getDepth()));
    Matrix3 m;
    
    float* hMap = (float*)mHeightMap->getData();

    // skip the outer edge and first vertex
    hMap += (patchSize+2 + 1) << 2;

    // loop through the vertices and save the normals to the normal map
    float* nMap = (float*)mNormalMap->getData(); 
    for(int y = 0; y < patchSize; ++y) {
        for(int x = 0; x < patchSize; ++x) {
            h1 = *(hMap + 7);
            h2 = *(hMap + ((patchSize + 2) * 4) + 3);
            h3 = *(hMap - 4 + 3);
            h4 = *(hMap - ((patchSize + 2) * 4) + 3);
                    
            spherePos.x = hMap[0];
            spherePos.y = hMap[1];
            spherePos.z = hMap[2];
        
            normal = Vector3( (h3 - h1), f, (h4 - h2)).normalisedCopy();
            
            buildRotationMatrix(Vector3::UNIT_Y, spherePos, m);
            normal = m * normal;
            
            *nMap++ = normal.x;
            *nMap++ = normal.y;
            *nMap++ = normal.z;
            
            hMap += 4;           
        }

        // skip edge vertices
        hMap += 8;
    }

Of course I plan to clean up the code when/if I make the loops OpenMP compatible.