Off by one error when creating face indices

Question

I seem to have an off by one error when creating the indices for the faces of a sphere and I can't understand where.

The exact problem is that when I try to calculate the normals of each face a dark region of nans appears in the region where the vertices wrap around.

I generate the vertices and face information like this:

pair<vector<vec3>, vector<uint>> GenerateSphere(
    float radius,
    uint u_resolution,
    uint v_resolution)
{
    vector<vec3> points;
    vector<uint> indices;
    for(uint i=0; i< u_resolution; i++)
    {
        float u = (float(i) / float(u_resolution - 1)) * M_PI - M_PI / 2.f;
        for(uint j=0; j < v_resolution; j++)
        {
            float v = float(j) * 2.f * M_PI / float(v_resolution - 1);

            points.push_back(vec3(cos(u) * cos(v), cos(u) * sin(v), sin(u)) * radius);
        }
    }

    for(uint i=0; i < u_resolution - 1; i++)
    {
        for(uint j=0; j < v_resolution; j++)
        {
            // Grab the 4 points in a quad face in the sphere
            uint p00 = (i * (v_resolution) + j);
            uint p01 = (i * (v_resolution) + ((j + 1) % v_resolution));
            uint p10 = ((i + 1) * v_resolution + j);
            uint p11 = ((i + 1) * v_resolution + ((j + 1) % v_resolution));

            assert((i + 1) < u_resolution);
            assert(p00 < points.size());
            assert(p01 < points.size());
            assert(p10 < points.size());
            assert(p11 < points.size());
            assert(
                p00 != p01 && p00 != p11 && p00 != p10 &&
                p01 != p10 && p01 != p11 &&
                p10 != p11
            );

            // Add the first triangle
            indices.push_back(p00);
            indices.push_back(p01);
            indices.push_back(p11);

            // Add the second triangle
            indices.push_back(p00);
            indices.push_back(p11);
            indices.push_back(p10);
        }
    }

    //for(auto& v : points) cout << to_string(v) << endl;
    //cout << endl;
    return {points, indices};
}

And immediately after calculate the normals like this:

std::vector<glm::vec3> GenerateNormals(const std::vector<glm::vec3>& vertices,
    const std::vector<uint>& indices)
{
    assert(indices.size() % 3 ==0);

    map<uint, vector<vec3>> normal_map;
    for(uint i=0; i < indices.size(); i += 3)
    {
        vec3 v1 = vertices[indices[i]];
        vec3 v2 = vertices[indices[i + 1]];
        vec3 v3 = vertices[indices[i + 2]];

        vec3 d1 = v2 - v1;
        vec3 d2 = v3 - v1;

        assert(indices[i] < vertices.size());
        assert(indices[i + 1] < vertices.size());
        assert(indices[i + 2] < vertices.size());
        assert(indices[i] != indices[i+1] && indices[i] != indices[i+2] &&
            indices[i+1] != indices[i+2]);

        if(abs(dot(d1, d2)) < 0.0001)
        {
            d1 = v3 - v2;
            d2 = v1 - v2;
        }

        if(abs(dot(d1, d2)) < 0.0001)
        {
            d1 = v1 - v3;
            d2 = v2 - v3;
        }

        vec3 normal = normalize(cross(d1, d2));

        normal_map[indices[i]].push_back(normal);
        normal_map[indices[i + 1]].push_back(normal);
        normal_map[indices[i + 2]].push_back(normal);
    }

    assert(normal_map.size() == vertices.size());
    std::vector<glm::vec3> normals(vertices.size());
    for(auto& pair : normal_map)
    {
        uint index = pair.first;
        vector<vec3> v_normals = pair.second;

        vec3 normal(0);
        for(auto& n : v_normals)
            normal += n;

        normal /= float(v_normals.size());
        normals[index] = normal;
    }
    assert(normals.size() == vertices.size());
    return normals;
}

Note: I know that I could just calculate the normals by taking the difference between the current vertex and the center of the sphere. These 2 functions I plan to use as templates to code other paramtetric surfaces, so that's not a solution for me I am trying to understand what I did wrong more than trying to fix it.

I have added lots of assertions through the code to try to catch the issue and I have also stepped through gdb but I can't seem to find the issue.

Answer 1

The problem with your code is, that it produces a lot of degenerated triangles. The not so obvious ones are at the top and bottom. The obvious ones are at the "wrapping line".

So what exactly is the problem? Each of your faces is a rectangle that consists of 2 triangles and 4 points. The 2 triangles share 2 points and each one has 1 individual point. But at the bottom/top 2 points are identical, so your rectangle degenerates to a triangle. As a result, one triangle degenerates to a line because of the duplicate points.

Now, why is this problematic for you? It is this line, that causes the trouble:

vec3 normal = normalize(cross(d1, d2));

The cross product of a vector with itself is the zero vector. If you normalize a vector, you divide by its length, which is 0 in this case. So this is where your NaNs come from. Until now, I just spoke about the top and bottom, but why is it happening at the "wrapping line"? The problem is a result of your index calculation:

    for(uint i=0; i < u_resolution - 1; i++)
    {
        for(uint j=0; j < v_resolution; j++)
        {
            // Grab the 4 points in a quad face in the sphere
            uint p00 = (i * (v_resolution) + j);
            uint p01 = (i * (v_resolution) + ((j + 1) % v_resolution));
            uint p10 = ((i + 1) * v_resolution + j);
            uint p11 = ((i + 1) * v_resolution + ((j + 1) % v_resolution));
   ...

The j loop is running until v_resolution. So j+1 gets problematic. I know your intention here is, that you need the triangles from the "other side", to calculate your normals correctly. That is why you performed the % v_resolution to get the corresponding points. But the problem is, you are not getting the points you want. With your code, you are including the last and first point of a "horizontal line", but they are identical in your case. The result is that you get a face that degenerates to a line, which means, that you will divide by 0 during normalization of both triangles surface normals.

A quick way to fix the NaNs at the "wrapping line" would be to make a temporary variable and store j+1 in it. Then check if it is equal or higher than v_resolution. If true, increase it by one. Use this variable everywhere you used j+1 during the index generation. However, this will still cause an issue at the wrapping line, since your start and end vertices are not the same and you only considered "both sides" for the last vertices of a horizontal line, but not for the first. This will lead to incorrect normals for the start vertices and you should see a jump in the shading. You might fix it, by just copying the normals of the last vertices to the first ones.

Also, you still need to address the problems at the top and bottom. This can be fixed by making a special loop only for top and bottom triangles. However, these are suboptimal fixes. A better approach would be to write code that does not produce duplicate points since it lowers the risk of creating degenerated triangles and is more memory efficient.