Calculating normals in a triangle mesh

c++ opengl computational-geometry normals

65,782

Solution 1

Does each vertex have its own normals or does each triangle have its own normals?

Like so often the answer is: "It depends". Since a normal is defined as being the vector perpendicular to all vectors within a given plane (in N dimensions), you need a plane to calculate a normal. A vertex position is just a point and thus singular, so you actually need a face to calculate the normal. Thus, naively, one could assume that normals are per face as the first step in normal calculation is determining the face normals, by evaluating the cross product of the faces edges.

Say you have a triangle with points A, B, C, then these points have position vectors ↑A, ↑B, ↑C and the edges have vectors ↑B - ↑A and ↑C - ↑A so the face normal vector is ↑N_f = (↑B - ↑A) × (↑C - ↑A)

Note that the magnitude of ↑N_f as it's stated above is directly proportional to the face's area.

In smooth surfaces vertices are shared between faces (or you could say those faces share a vertex). In that case the normal at the vertex is not one of the face normals of the faces it is part of, but a linear combination of them:

↑N_v = ∑ p ↑N_f ; where p is a weighting for each face.

One could either assume a equal weighting between the participating face normals. But it makes more sense to assume that the larger a face is, the more it contributes to the normal.

Now recall that you normalize by a vector ↑v by scaling it with it's recipocal length: ↑v_i = ↑v/|↑v|. But as already told the length of the face normals already depends on the face's area. So the weighting factor p given above is already contained in the vector itself: Its length, aka magnitude. So we can get the vertex normal vector by simply summing up all the face normals.

In lighting calculations the normal vector must be unit length, i.e. normalized to be useable. So after summing up, we normalize the newly found vertex normal and use that.

The carefull reader may have noticed I specifically said smooth surfaces share vertices. And in fact, if you have some creases / hard edges in your geometry, then the faces on either side don't share vertices. In OpenGL a vertex is the whole combination of

position
normal
(colour)
N texture coordinates
M further attributes

You change one of these and you got a completely different vertex. Now some 3D modelers see a vertex only as a point's position and store the rest of those attributes per face (Blender is such a modeler). This saves some memory (or considerable memory, depending on the number of attributes). But OpenGL needs the whole thing, so if working with such a mixed paradigm file you will have to decompose it into OpenGL compatible data first. Have a look at one of Blender's export scripts, like the PLY exporter to see how it's done.

Now to cover some other thing. In your code you have this:

 glIndexPointer( GL_UNSIGNED_BYTE, 0, indices );

The index pointer has nothing to do with vertex array indices! This is an anachronsim from the days, when graphics still used palettes instead of true color. A pixels colour wasn't set by giving it's RGB values, but by a single number offsetting into a limited palette of colours. Palette colours can still be found in several graphics file formats, but no decent piece of hardware uses them anymore.

Please erase glIndexPointer (and glIndex) from your memory and your code, they don't do what you think they do The whole indexed color mode is arcane to used, and frankly I don't know of any hardware built after 1998 that still supported it.

Solution 2

Thumbs up for datenwolf! I completely agree with his approach. Adding the normal vectors of the adjacent triangles for each vertex and then normalising is the way to go. I just want to push the answer a little bit and have a closer look at the particular but quite common case of a rectangular, smooth mesh that has a constant x/y step. In other words, a rectangular x/y grid with a variable height at each point.

Such a mesh is created by looping over x and y and setting a value for z and can represent things like the surface of a hill. So each point of the mesh is represented by a vector

P = (x, y, f(x,y))

where f(x,y) is a function giving the z of each point on the grid.

Usually to draw such a mesh we use a TriangleStrip or a TriangleFan but any technique should give a similar topography for the resulting triangles.

     |/   |/   |/   |/
...--+----U----UR---+--...
    /|   /| 2 /|   /|           Y
   / |  / |  / |  / |           ^
     | /  | /  | /  | /         |
     |/ 1 |/ 3 |/   |/          |
...--L----P----R----+--...      +-----> X
    /| 6 /| 4 /|   /|          
   / |  / |  / |  / |         
     | /5 | /  | /  | /      
     |/   |/   |/   |/
...--DL---D----+----+--...
    /|   /|   /|   /|

For a triangleStrip each vertex P=(x0, y0, z0) has 6 adjacent vertices denoted

up       = (x0     , y0 + ay, Zup)
upright  = (x0 + ax, y0 + ay, Zupright) 
right    = (x0 + ax, y0     , Zright) 
down     = (x0     , y0 - ay, Zdown)
downleft = (x0 - ax, y0 - ay, Zdownleft) 
left     = (x0 - ax, y0     , Zleft)

where ax/ay is the constant grid step on the x/y axis respectively. On a square grid ax = ay.

ax = width / (nColumns - 1)
ay = height / (nRows - 1)

Thus each vertex has 6 adjacent triangles each one with its own normal vector (denoted N1 to N6). These can be calculated using the cross product of the two vectors defining the side of the triangle and being careful on the order in which we do the cross product. If the normal vector points in the Z direction towards you :

N1 = up x left =
   = (Yup*Zleft - Yleft*Zup, Xleft*Zup - Xup*ZLeft, Xleft*Yup - Yleft*Xup) 

   =( (y0 + ay)*Zleft - y0*Zup, 
      (x0 - ax)*Zup   - x0*Zleft, 
      x0*y0 - (y0 + ay)*(x0 - ax) ) 

N2 = upright  x up
N3 = right    x upright
N4 = down     x right
N5 = downleft x down
N6 = left     x downleft

And the resulting normal vector for each point P is the sum of N1 to N6. We normalise after summing. It's very easy to create a loop, calculate the values of each normal vector, add them and then normalise. However, as pointed out by Mr. Shickadance, this can take quite a while, especially for large meshes and/or on embedded devices.

If we have a closer look and perform the calculations by hand, we will find out that most of the terms cancel out each other, leaving us with a very elegant and easy to calculate final solution for the resulting vector N. The point here is to speed up calculations by avoiding calculating the coordinates of N1 to N6, doing 6 cross-products and 6 additions for each point. Algebra helps us to jump straight to the solution, use less memory and less CPU time.

I will not show the details of the calculations as it is long but straight-forward and will jump to the final expression of the Normal vector for any point on the grid. Only N1 is decomposed for the sake of clarity, the other vectors look alike. After summing we obtain N which is not yet normalized :

N = N1 + N2 + ... + N6

  = .... (long but easy algebra) ...

  = ( (2*(Zleft - Zright) - Zupright + Zdownleft + Zup - Zdown) / ax,
      (2*(Zdown - Zup)    + Zupright + Zdownleft - Zup - Zleft) / ay,
       6 )

There you go! Just normalise this vector and you have the normal vector for any point on the grid, provided you know the Z values of its surrounding points and the horizontal/vertical step of your grid.

Note that this is the weighed average of the surrounding triangles' normal vectors. The weight is the area of the triangles and is already included in the cross product.

You can even simplify it more by only taking into account the Z values of four surrounding points (up,down,left and right). In that case you get :

                                             |   \|/   |
N = N1 + N2 + N3 + N4                    ..--+----U----+--..
  = ( (Zleft - Zright) / ax,                 |   /|\   |
      (Zdown -  Zup  ) / ay,                 |  / | \  |
       2 )                                 \ | / 1|2 \ | /
                                            \|/   |   \|/
                                         ..--L----P----R--...
                                            /|\   |   /|\
                                           / | \ 4|3 / | \
                                             |  \ | /  |
                                             |   \|/   |
                                         ..--+----D----+--..
                                             |   /|\   |

which is even more elegant and even faster to calculate.

Hope this will make some meshes faster. Cheers

Solution 3

Per-vertex.

Use cross-products to calculate the face normals for the triangles surrounding a given vertex, add them together, and normalize.

Solution 4

As simple as it may seem, calculating the normal of a triangle is only part of the problem. The cross product of 2 sides of the polygon is sufficient in triangular cases, unless the triangle is collapsed onto itself and degenerate; in that case there is no one valid normal, so you can select one to your liking.

So why is the normalized cross product only part of the problem? The winding order of the vertices in that polygon defines the direction of the normal, i.e. if one pair of vertices is swapped in place, the normal will point in the opposite direction. So in fact this can be problematic if the mesh itself contains inconsistencies in that regard, i.e. parts of it assume one ordering, while other parts assume different orderings. One famous example is the original Stanford Bunny model, where some parts of the surface will point inwards, while others point outwards. The reason for that is because the model was constructed using a scanner, and no care was taken to produce triangles with regular winding patterns. (obviously, clean versions of the bunny also exist)

The winding problem is even more prominent if polygons can have multiple vertices, because in that case you would be averaging partial normals of the semi-triangulation of that polygon. Consider the case where partial normals are pointing in opposite directions, resulting in normal vectors of length 0 when taking the mean!

In the same sense, disconnected polygon soups and point clouds present challenges for accurate reconstruction due to the ill-defined winding number.

One potential strategy that is often used to solve this problem is to shoot random rays from outward to the center of each semi-triangulation (i.e. ray-stabbing). But one cannot assume that the triangulation is valid if polygons can contain multiple vertices, so rays may miss that particular sub-triangle. If a ray hits, then normal opposite to the ray direction, i.e. with dot(ray, n) < .5 satisfied, can be used as the normal for the entire polygon. Obviously this is rather expensive and scales with the number of vertices per polygon.

Thankfully, there's great new work that describes an alternative method that is not only faster (for large and complex meshes) but also generalizes the 'winding order' concept for constructions beyond polygon meshes, such as point clouds and polygon soups, iso-surfaces, and point-set surfaces, where connectivity may not even be defined!

As outlined in the paper, the method constructs a hierarchical splitting tree representation that is refined progressively, taking the parent 'dipole' orientation into account at every split operation. A polygon normal would then simply be an integration (mean) over all di-poles (i.e. point+normal pairs) of the polygon.

For people who are dealing with unclean mesh/pcl data from Lidar scanners or other sources, this could def. be a game-changer.

Solution 5

For those like me who came across this question, your answer might be this :

// Compute Vertex Normals
std::vector<sf::Glsl::Vec3> verticesNormal;
verticesNormal.resize(verticesCount);

for (i = 0; i < indices.size(); i += 3)
{
    // Get the face normal
    auto vector1 = verticesPos[indices[(size_t)i + 1]] - verticesPos[indices[i]];
    auto vector2 = verticesPos[indices[(size_t)i + 2]] - verticesPos[indices[i]];
    auto faceNormal = sf::VectorCross(vector1, vector2);
    sf::Normalize(faceNormal);

    // Add the face normal to the 3 vertices normal touching this face
    verticesNormal[indices[i]] += faceNormal;
    verticesNormal[indices[(size_t)i + 1]] += faceNormal;
    verticesNormal[indices[(size_t)i + 2]] += faceNormal;
}

// Normalize vertices normal
for (i = 0; i < verticesNormal.size(); i++)
    sf::Normalize(verticesNormal[i]);

View more solutions

65,782

Author by

Vince

Updated on September 04, 2020

Comments

Vince over 3 years

I have drawn a triangle mesh with 10000 vertices(100x100) and it will be a grass ground. I used gldrawelements() for it. I have looked all day and still can't understand how to calculate the normals for this. Does each vertex have its own normals or does each triangle have its own normals? Can someone point me in the right direction on how to edit my code to incorporate normals?

struct vertices {
    GLfloat x;
    GLfloat y;
    GLfloat z;
}vertices[10000];

GLuint indices[60000];

/*
99..9999
98..9998
........
01..9901
00..9900
*/

void CreateEnvironment() {
    int count=0;
    for (float x=0;x<10.0;x+=.1) {
        for (float z=0;z<10.0;z+=.1) {
            vertices[count].x=x;
            vertices[count].y=0;
            vertices[count].z=z;
            count++;
        }
    }
    count=0;
    for (GLuint a=0;a<99;a++){
        for (GLuint b=0;b<99;b++){
            GLuint v1=(a*100)+b;indices[count]=v1;count++;
            GLuint v2=(a*100)+b+1;indices[count]=v2;count++;
            GLuint v3=(a*100)+b+100;indices[count]=v3;count++;
        }
    }
    count=30000;
    for (GLuint a=0;a<99;a++){
        for (GLuint b=0;b<99;b++){
            indices[count]=(a*100)+b+100;count++;//9998
            indices[count]=(a*100)+b+1;count++;//9899
            indices[count]=(a*100)+b+101;count++;//9999
        }
    }
}

void ShowEnvironment(){
    //ground
    glPushMatrix();
    GLfloat GroundAmbient[]={0.0,0.5,0.0,1.0};
    glMaterialfv(GL_FRONT,GL_AMBIENT,GroundAmbient);
    glEnableClientState(GL_VERTEX_ARRAY);
    glIndexPointer( GL_UNSIGNED_BYTE, 0, indices );
    glVertexPointer(3,GL_FLOAT,0,vertices);
    glDrawElements(GL_TRIANGLES,60000,GL_UNSIGNED_INT,indices);
    glDisableClientState(GL_VERTEX_ARRAY);
    glPopMatrix();
}

EDIT 1 Here is the code I have written out. I just used arrays instead of vectors and I stored all of the normals in the struct called normals. It still doesn't work however. I get an unhandled exception at *indices.

struct Normals {
    GLfloat x;
    GLfloat y;
    GLfloat z;
}normals[20000];
Normals* normal = normals;
//***************************************ENVIRONMENT*************************************************************************
struct vertices {
    GLfloat x;
    GLfloat y;
    GLfloat z;
}vertices[10000];

GLuint indices[59403];

/*
99..9999
98..9998
........
01..9901
00..9900
*/

void CreateEnvironment() {
    int count=0;
    for (float x=0;x<10.0;x+=.1) {
        for (float z=0;z<10.0;z+=.1) {
            vertices[count].x=x;
            vertices[count].y=rand()%2-2;;
            vertices[count].z=z;
            count++;
        }
    }
    //calculate normals 
    GLfloat vector1[3];//XYZ
    GLfloat vector2[3];//XYZ
    count=0;
    for (int x=0;x<9900;x+=100){
        for (int z=0;z<99;z++){
            vector1[0]= vertices[x+z].x-vertices[x+z+1].x;//vector1x
            vector1[1]= vertices[x+z].y-vertices[x+z+1].y;//vector1y
            vector1[2]= vertices[x+z].z-vertices[x+z+1].z;//vector1z
            vector2[0]= vertices[x+z+1].x-vertices[x+z+100].x;//vector2x
            vector2[1]= vertices[x+z+1].y-vertices[x+z+100].y;//vector2y
            vector2[2]= vertices[x+z+1].z-vertices[x+z+100].z;//vector2z
            normals[count].x= vector1[1] * vector2[2]-vector1[2]*vector2[1];
            normals[count].y= vector1[2] * vector2[0] - vector1[0] * vector2[2];
            normals[count].z= vector1[0] * vector2[1] - vector1[1] * vector2[0];count++;
        }
    }
    count=10000;
    for (int x=100;x<10000;x+=100){
        for (int z=0;z<99;z++){
            vector1[0]= vertices[x+z].x-vertices[x+z+1].x;//vector1x -- JUST ARRAYS
            vector1[1]= vertices[x+z].y-vertices[x+z+1].y;//vector1y
            vector1[2]= vertices[x+z].z-vertices[x+z+1].z;//vector1z
            vector2[0]= vertices[x+z+1].x-vertices[x+z-100].x;//vector2x
            vector2[1]= vertices[x+z+1].y-vertices[x+z-100].y;//vector2y
            vector2[2]= vertices[x+z+1].z-vertices[x+z-100].z;//vector2z
            normals[count].x= vector1[1] * vector2[2]-vector1[2]*vector2[1];
            normals[count].y= vector1[2] * vector2[0] - vector1[0] * vector2[2];
            normals[count].z= vector1[0] * vector2[1] - vector1[1] * vector2[0];count++;
        }
    }

    count=0;
    for (GLuint a=0;a<99;a++){
        for (GLuint b=0;b<99;b++){
            GLuint v1=(a*100)+b;indices[count]=v1;count++;
            GLuint v2=(a*100)+b+1;indices[count]=v2;count++;
            GLuint v3=(a*100)+b+100;indices[count]=v3;count++;
        }
    }
    count=30000;
    for (GLuint a=0;a<99;a++){
        for (GLuint b=0;b<99;b++){
            indices[count]=(a*100)+b+100;count++;//9998
            indices[count]=(a*100)+b+1;count++;//9899
            indices[count]=(a*100)+b+101;count++;//9999
        }
    }
}

void ShowEnvironment(){
    //ground
    glPushMatrix();
    GLfloat GroundAmbient[]={0.0,0.5,0.0,1.0};
    GLfloat GroundDiffuse[]={1.0,0.0,0.0,1.0};
    glMaterialfv(GL_FRONT,GL_AMBIENT,GroundAmbient);
    glMaterialfv(GL_FRONT,GL_DIFFUSE,GroundDiffuse);
    glEnableClientState(GL_VERTEX_ARRAY);
    glEnableClientState(GL_NORMAL_ARRAY);
    glNormalPointer( GL_FLOAT, 0, normal);
    glVertexPointer(3,GL_FLOAT,0,vertices);
    glDrawElements(GL_TRIANGLES,60000,GL_UNSIGNED_INT,indices);
    glDisableClientState(GL_VERTEX_ARRAY);
    glDisableClientState(GL_NORMAL_ARRAY);
    glPopMatrix();
}
//***************************************************************************************************************************

Wroclai almost 13 years

Very correct. :-) Instead of adding my own answer I'll just add a comment about this answer in which I have given a detailed description on calculating the normals.
Mr. Shickadance about 12 years

Is there a more efficient way of doing this? It works OK for simple models, but once I have a large number of vertices it takes a long time to calculate all the per-vertex normals.
jozxyqk over 10 years

Then there's all the edge cases, for example triangles with zero area or thin triangles that skews the weight of the normal - weighting the normals by triangle area (summing the cross product un-normalized and only normalizing at the end) is possible but doesn't always give perfect results.
unlut over 2 years

I don't get why, for example, normal vector for first triangle is equal to "N1 = up x left". Shouldn't the calculation performed by cross-producting two edges of triangle, these edges could be E1 = (U-P) and E2 = (L-P), then normal should be equal to "N1 = (U-P) x (L-P)".