How to Calculate Surface Normals for Generated Geometry

Question

I have a class that generates a 3D shape based on inputs from the calling code. The inputs are things like length, depth, arc, etc. My code generates the geometry perfectly, however I am running into trouble when calculating the surface normals. When lit, my shape has very bizarre coloring/texture from the incorrect surface normals that are being calculated. From all my research I believe my math is correct, it seems that something is wrong with my technique or method.

At a high level how does one go about programmatically calculating the surface normals for a generated shape? I am using Swift/SceneKit on iOS for my code but a generic answer is fine.

I have two arrays that represent my shape. One is an array of 3D points that represents the vertices that make up the shape. The other array is a list of indexes of the first array that map the vertices into triangles. I need to take that data and generate a 3rd array that is a set of surface normals that aid in the lighting of the shape. (see SCNGeometrySourceSemanticNormal in SceneKit`)

The list of vertices and indexes is always different depending on the inputs to the class so I cannot pre-calculate or hard code the surface normals.

Answer 1

You simply dont want fully smooth results. While the commented method by Nathan Reed: "Calculate each vertex to face normal, sum them, normalize sum", generally works it sometimes fails spectacularly. But that is of no importance here, we can use that method by adding a rejection clause to it.

In this case you simply want certain parts not to be smoothed against certain other parts. You want selective hard edges. So for example the flat top and bottom is separate form the triangle strip on the side, as is each flat area.

Image 1: The result you want.

In effect you only want to average the vertices of the curved area all others can use the normal they get form their triangle alone. So you are better of thinking of the mesh as 9 separate regions that are handled without the others.

Image 2: Image showing the mesh structure and the normals.

You can certainly automatically deduce this by not including normals that are outside certain angle from the primary vertexes normal. Pseudocode:

For vertex in faceVertex:
    normal = vertex.normal
    For adjVertex in adjacentVertices:
        if anglebetween(vertex.normal, adjVertex.normal )  < treshold:
            normal += adjVertex.normal
    normal = normalize(normal)

That works, but you can simply avoid all of this at creation time because you understand that separate planes are working differently. So only the curved sides need normal direction merging. And in fact you can just directly calulate them from the underlying mathematical shape.

Answer 2

I see mainly three ways of computing normals for a generated shape.

Analytic normals

In some cases you have enough information about the surface to generate the normals. For example, the normal of any point on a sphere is trivial to compute. Put simply, when you know the derivative of the function, you also know the normal.

If your case is narrow enough to allow you to use analytic normals, they will probably give the best result in terms of precision. The technique doesn't scale too well though: if you also need to handle cases where you cannot use analytic normals, it may be easier to keep the technique that handles the general case and drop the analytic one altogether.

Vertex normals

The cross product of two vectors gives a vector perpendicular to the plane they belong to. So getting the normal of a triangle is straightforward:

vec3 computeNormal(vec3 a, vec3 b, vec3 c)
{
    return normalize(crossProduct(b - a, c - a));
}

Moreover, in the above example, the length of the cross product is proportional to the area inside abc. So the smoothed normal at a vertex shared by several triangles can be computed by summing up the cross products and normalizing as a last step, thus weighting each triangle by its area.

vec3 computeNormal(vertex a)
{
    vec3 sum = vec3(0, 0, 0);
    list<vertex> adjacentVertices = getAdjacentVertices(a);
    for (int i = 1; i < adjacentVertices; ++i)
    {
        vec3 b = adjacentVertices[i - 1];
        vec3 c = adjacentVertices[i];
        sum += crossProduct(b - a, c - a);
    }
    if (norm(sum) == 0)
    {
        // Degenerate case
        return sum;
    }
    return normalize(sum);
}

If you are working with quads, there is a nice trick you can use: for a quad abcd, use crossProduct(c - a, d - b) and it will handle nicely cases where the quad is in fact a triangle.

Iñigo quilez wrote a few short articles on the topic: clever normalization of a mesh, and normal and area of n sided polygons.

Normals from partial derivatives

Normals can be computed in the fragment shader from the partial derivatives. The math behind is the same, except this time it is done in screen space. This article by Angelo Pesce describes the technique: Normals without normals.