C++ code for cotangent weights over a triangular mesh

Here is some code to compute the cotangent weights described in the article Discrete Differential-Geometry Operators for Triangulated 2-Manifolds. Cotangent weights have many application such as mesh smoothing, computing harmonic weights over the mesh, or diffusing values on its surface. I'm leaving this here for other tutorial to reference.

A cotangent weight can be seen as barycentric weights, for each vertex we associate a weight to each of its direct neighbor.

/// compute the cotangent weight of the neighbor edge_curr of curr_vert
float laplacian_cotan_weight(const Mesh_topology& m,
                             mesh_t::Vert_idx curr_vert,
                             int edge_curr )
{
    const Mesh_geometry& geom = m.geom();
    const int nb_neighs = m.first_ring()[curr_vert].size();
    const Vec3 c_pos = geom.vertex( curr_vert );

    const int edge_next = (edge_curr+1) % nb_neighs;

    const int id_curr = m.first_ring()[curr_vert][ edge_curr ];
    const int id_next = m.first_ring()[curr_vert][ edge_next ];
    const int id_prev = m.first_ring()[curr_vert][ (edge_curr-1) >= 0 ? (edge_curr-1) : nb_neighs-1 ];

    const Vec3 v1 = c_pos                - geom.vertex(id_prev);
    const Vec3 v2 = geom.vertex(id_curr) - geom.vertex(id_prev);
    const Vec3 v3 = c_pos                - geom.vertex(id_next);
    const Vec3 v4 = geom.vertex(id_curr) - geom.vertex(id_next);

    // wij = (cot(alpha) + cot(beta)),
    // for boundary edge, there is only one such edge
    float cotan1 = 0.0f;
    float cotan2 = 0.0f;
    // If the mesh is not a water-tight closed volume 
    // we must check for edges lying on the sides of wholes
    if( !m.is_side_edge( m.first_ring_edges(curr_vert)[edge_curr] ) )
    {
        // general case: not a boundary
        cotan1 = (v1.dot(v2)) / (v1.cross(v2)).norm();
        cotan2 = (v3.dot(v4)) / (v3.cross(v4)).norm();
    }
    else // boundary edge, only have one such angle
    {
        if( id_next == id_prev )
        {
            // two angles are the same, e.g. corner of a square
            cotan1 = (v1.dot(v2)) / (v1.cross(v2)).norm();
        }
        else
        {
            // find the angle not on the boundary
            if( !m.is_side_edge( m.first_ring_edges(curr_vert)[edge_next] ) )
                cotan2 = (v3.dot(v4)) / (v3.cross(v4)).norm();
            else
                cotan1 = (v1.dot(v2)) / (v1.cross(v2)).norm();
        }
    }

    return (cotan1 + cotan2);
}

// -----------------------------------------------------------------------------

void laplacian_cotan_weights(const Mesh_topology& m,
                             std::vector< std::vector >& per_edge_cotan_weights)
{
    int nb_verts = m.geom().nb_vertices();
    per_edge_cotan_weights.resize( nb_verts );
    for(int i = 0; i < nb_verts; ++i)
    {
        const int nb_neighs = (int)m.first_ring()[i].size();
        per_edge_cotan_weights[i].resize( nb_neighs );
        for(int n = 0; n < nb_neighs; ++n)
        {
            float w = laplacian_cotan_weight(m, i, n);
            //std::cout << "vert: " << i << " cotan wij " << w << std::endl;
            per_edge_cotan_weights[i][n] = w;
        }
    }
}

// -----------------------------------------------------------------------------

struct Mesh_topology {

    /// Per vertex list of the 1st ring of vertex neighbors
    /// first_ring[id_vertex][ith_vert_neighbor] == vert_neighbor_id
    /// @note each list of 1st ring neighbors respect their original geometric order.
    std::vector< std::vector<mesh_t::Vert_idx> > first_ring();


    /// Get a list of outward edges from the vertex 'i'. 
    /// @note same order as with first ring of vertices (first_ring())
    std::vector& first_ring_edges(mesh_t::Vert_idx i);
};


// -----------------------------------------------------------------------------

class Mesh_geometry {
public:
    inline const Vec3& vertex(mesh_t::Vert_idx idx) const { return _vertices[idx]; }
    inline int nb_vertices() const { return (int)_vertices.size(); }
public:
    /// List of vertices
    std::vector _vertices;
    /// List of triangles
    std::vector<mesh_t::Tri_face> _triangles;
};

This piece of code does not compile but comes from a perfectly working project. It should help you rapidely adapt this code to your own project.

Some other papers than can make use of those weights:

No comments

All html tags except <b> and <i> will be removed from your comment. You can make links by just typing the url or mail-address.
Spam bot question: