Laplacian smoothing (C++ code to smooth a mesh)

Dropping a procedure to compute the Laplacian smoothing of a 3D mesh (with cotangent weights).

Note that the cotangent weights can be replaced by w=1.0. In this case the smoothing will be said to be a "uniform Laplacian smoothing" and the initial distribution of the triangles will not be preserved. Of course you can test other barycentric weights. The code below use the explicit scheme to compute the Laplacian smoothing, explicit scheme is known to be unstable (especially with cotangent weights) for t > 1. On the other hand Implicit scheme of the Laplacian smoothing is unconditionally stable for any t! But you need to solve a linear system of equations

C++ code for the explicit and implicit scheme ].

/// @brief Smoothing while minimizing the texture distortions.
class Weighted_laplacian {
    // Some internal buffer for the computation
    std::vector<Vec3> _buffer_vertices;
    // Cache the cotangent weights
    std::vector< std::vector<float> > _cotan_weights;
    // Store the list of first ring neighbor for each vertex.
    // the first ring neighbor is the list of vertices directly adjacent to another vertex
    const Mesh_topology* _topo;
    /// @note avoid odd numbers of iterations they imply a full copy of
    /// the vertex buffer at the end of the smoothing process.
    int _nb_iter;

    /// If the **rest** position of the mesh changes you need to update the
    /// pre-computation step.
    void update(const Mesh_topology& rest_pose_mesh)
        _topo = &topo;
        mesh_utils::laplacian_cotan_weights(topo, _cotan_weights);

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

    /// @param vertices : in place smoothing of this vertex buffer
    void process(Vec3 vertices[],
                 float smooth_factors[],
                 int nb_vertices)
        _buffer_vertices.resize( nb_vertices );
        if(_nb_iter == 0){
        // first_ring[a_vertex][ith_neighbor of a_vertex] = 
        // "index of a neighbor directly connected to a_vertex"
        const std::vector< std::vector< mesh_t::Vert_idx > >& first_ring = _topo->first_ring();
        Vec3* src_vertices = vertices;
        Vec3* dst_vertices =;
        for(int k = 0; k < _nb_iter; k++)
            for( int i = 0; i < nb_vertices; i++)
                Vec3 cog(0.f);
                float sum = 0.f;
                size_t nb_neighs = first_ring[i].size();
                for(size_t n = 0; n < nb_neighs; n++)
                    mesh_t::Vert_idx neigh = first_ring[i][n];
                    // Note: if we set w = 1.0f then the procedures operates
                    // an uniform smoothing, the initial distribution of 
                    // the triangles will not be preserved
                    float w = _cotan_weights[i][n]; 
                    cog += src_vertices[neigh] * w;
                    sum += w;
                // When using cotan weights smoothing may be unstable 
                // in this case we need to set t < 1
                // sometimes you even need to get as low as t < 0.5
                float t = smooth_factors[i];                
                dst_vertices[i] = (cog / sum) * t + src_vertices[i] * (1.f - t);
            std::swap(dst_vertices, src_vertices);
        if(_nb_iter%2 == 1)            
            utils::copy( vertices/*destination*/,*source*/, nb_vertices );

(draft notes)

Smooth a vertex:

cog_i = \sum\limits_{{j \in \mathcal N(i)}}  w_{ij} p_j

Damping (linear interpolation with previous position) \( cog_i \lambda+ p_i (1-\lambda) \)

re-write the expression so that we add a displacement vector to the original vertex:
cog_i \lambda+ p_i (1-\lambda) \\
cog_i \lambda+ p_i - p_i \lambda \\
p_i + cog_i \lambda - p_i \lambda \\
p_i + \lambda (cog_i - p_i) \\

Let \( \mathbf{L} \) be the Laplacian matrix and define \( \tilde { \mathbf{L} }= \mathbf D_l^{-1} \mathbf L \) with \(  \mathbf D_l \) the diagonal elements of \( \mathbf L \) and \( \mathbf V_k \) be the vertices after the k iteration of smoothing.

One iteration of Laplacian expressed with matrices:

(\mathbf I +\lambda \tilde{ \mathbf L })  \mathbf V_k =  \mathbf V_{k+1}

This explicit scheme is unstable

Note: some authors chose to use \( -\mathbf L \) (positive diagonal, negative adjacent elements) then the expression becomes \( (\mathbf I - \lambda \tilde{ \mathbf L })  \mathbf V \)

You can precompute p iterations by multiplying the matrix to itself \( (\mathbf I +\lambda \tilde{ \mathbf L })^p  \)

Now if you solve the linear system of equations:
(\mathbf I - \lambda \tilde{ \mathbf L })  \mathbf V_k = \mathbf V_{k-1} \\
(\mathbf I - \lambda \tilde{ \mathbf L })  \mathbf X = \mathbf V

this is unconditionally stable for any \( \lambda \) :) (i.e. no ugly mesh explosion)

Related smoothing method

Taubin smoothing


Implicit Fairing of Irregular Meshes using Diffusion and Curvature Flow (Mathieu Desbrun & al)


Mean curvature flow and pitfalls.
Curvature flow C++ implementation
Conformal Willmore flow.

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.
Anti-spam question: