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.

[![enter image description here][1]][1]

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


  [1]: https://i.stack.imgur.com/yh3pn.png