Implicit skinning routines

Here I will progressively gather bits and pieces of code related to Implicit skinning.

Other pages

Some pieces are described more in depth in separate posts:

Bulge Free DQS

There is a small trick that can be used to correct the bulge of Dual Quaternion Skinning prior to the Implicit Skinning. It may help increase performance and avoid some miss-projection during implicit skinning. There is also the Disney paper Real-time Skeletal Skinning with Optimized Centers of Rotation which also reduce the bulge artefact in a more clean and robust way, haven't tested out yet but it looks like the nicest input for implicit skinning 2013 so far.

Projection algorithm

The projection algorithm corrects the shape of the mesh according to an implicit surface. The following excerpt of code can be used to reproduce the 2013 paper. Note that in this version we project using a fixed step's length, newton iterations (steps' length base on the gradient's norm) should be more efficient.

Inputs:

struct Ray {
    Point3 _pos;
    Vec3  _dir;

    Ray() : _pos( Point3(0.f, 0.f, 0.f) ), _dir( Vec3::zero() ) {  }

    Ray(const Point3& pos_, const Vec3& dir_) :
        _pos( pos_ ), _dir( dir_ )
    { }

    void set_pos(const Point3& pos_) { _pos = pos_; }
    void set_dir(const Vec3& dir_) { _dir = dir_; }
    Point3 operator ()(const float t_) const {
        return _pos + _dir * t_;
    }
};

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

class Implicit_base  {
public:
    Implicit_base() { }
    virtual ~Implicit_base(){ }

    // -------------------------------------------------------------------------
    /// @name evaluate potential field and gradient:
    
    /// @brief evaluate the scalar 'field' -> (f)
    virtual float f(const Vec3& global_pos) const = 0;

    /// @brief evaluate 'gradient field' -> (gf)
    virtual Vec3 gf(const Vec3& global_pos) const = 0;

    /// @brief evaluate 'field and gradient field' -> (fngf)
    virtual float fngf(const Vec3& global_pos, Vec3& grad) const = 0;
};

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

class Implicit_skinning_2013 : public Implicit_projection_base {
public:

    /// Flags defining the states a vertex go through when projected
    /// This helps to see what is happening and debug.
    enum Vertex_state_e {
        /// The vertex is going away from the target potential instead of getting closer
        ePOTENTIAL_PIT = 0,     
        /// A collision is detected (skin contact).
        eGRADIENT_DIVERGENCE,   
        /// stopped before reaching base potential. 
        /// (Not enough iteration to march until we reach the target iso)
        eNB_ITER_MAX,           
        /// The vertex is already on the correct iso-surface
        eNOT_DISPLACED,         
        /// The vertex has been succesfully projected on its base potential
        eFITTED,                
        /// Gradient norm is null -> we can't set a proper marching direction
        eNORM_GRAD_NULL,        
        eUNDEF,        
        eNB_CASES
    };

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

    Implicit_skinning_2013()
        : _max_nb_iter(1000)
        , _gradient_threshold(0.8f)
        , _step_length(0.01f)
        , _potential_pit(true)
    {

    }

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

    void project_vertex(
            const Implicit_base& csg_tree,
            const Vec3& in_vertex,
            float base_potential,
            Vec3& out_vertex,
            Vec3& out_gradient,
            Vertex_state_e& vertex_state,
            float eps = 0.001f);

    // -------------------------------------------------------------------------
    /// @name Projection parameters
    // -------------------------------------------------------------------------
    /// max number of steps to march towards the base potential
    unsigned _max_nb_iter;         
    /// maximum gradient angle authorized between two projection steps
    float    _gradient_threshold;  
    /// size of a single step
    float    _step_length;         
    /// wether the projection stops when the vertex marches backwards to its goal.
    bool     _potential_pit;       
};

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

void Implicit_skinning_2013::project_vertex(
        const Implicit_base& csg_tree,
        const Vec3& in_vertex,
        float base_potential,
        Vec3& out_vertex,
        Vec3& out_gradient,
        Implicit_skinning_2013::Vertex_state_e& vertex_state,
        float eps)
{
    const float ptl = base_potential;

    Vec3 v0 = in_vertex;
    Vec3 gf0;
    float f0;
    f0 = csg_tree.fngf(v0, gf0) - ptl;

    out_gradient = gf0;
    // STOP CASE : Gradient is null we can't know where to march
    if(gf0.norm() <= 0.00001f){
        vertex_state = eNORM_GRAD_NULL;
        return;
    }

    // STOP CASE : Point already near enough the isosurface
    if( std::abs( f0 ) < eps ){
        vertex_state = eNOT_DISPLACED;
        return;
    }

    vertex_state = eNB_ITER_MAX;

    // Inside we march along the inverted gradient
    // outside along the gradient :
    const float dl = (f0 > 0.f) ? -_step_length : _step_length;

    Ray r;
    r.set_pos( v0.to_point3() );
    float t = 0.f;

    Vec3  gfi = gf0;
    float fi  = f0;
    Vec3  vi  = v0;
    for(unsigned i = 0; i < _max_nb_iter; ++i)
    {

        r.set_pos( v0.to_point3() );
        r.set_dir( gf0.normalized() );
        t = dl;

        vi = r(t);
        fi = csg_tree.fngf(vi, gfi) - ptl;

        // STOP CASE 1 : Initial iso-surface reached
        if( fi * f0 <= 0.f)
        {
            /// todo update gradient to match the new position of dichotomic search
            t = dichotomic_search( r, 0.f, t, ptl, csg_tree, eps, 20);
            vertex_state = eFITTED;
            break;
        }

        // STOP CASE 2 : Gradient divergence
        if( (gf0.normalized()).dot(gfi.normalized()) < _gradient_threshold)
        {
            vertex_state = eGRADIENT_DIVERGENCE;
            break;
        }

        // STOP CASE 3 : Potential pit
        if( ((fi - f0)*dl < 0.f) && _potential_pit )
        {
            vertex_state = ePOTENTIAL_PIT;
            break;
        }

        v0  = vi;
        f0  = fi;
        gf0 = gfi;

        if(gf0.norm_squared() < (0.001f*0.001f))
            break;
    }

    const Point3 res = r(t);
    out_gradient = gfi;
    out_vertex   = res.to_vec3();
}

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

template< class Object >
float dichotomic_search(const Ray& ray,
                        float t0, float t1,
                        float target_iso,
                        const Object& obj,
                        float eps,
                        unsigned max_steps)
{
    float f_mid;
    float t = t0;
    float f0 = obj.f( ray(t0).to_vec3() );
    float f1 = obj.f( ray(t1).to_vec3() );
    assert( std::max( f0, f1 ) > target_iso &&
            std::min( f0, f1 ) < target_iso &&
            "target_iso must be within [f0, f1]" );

    if(f0 > f1){
        t0 = t1;
        t1 = t;
    }

    Vec3 p;
    for(unsigned i = 0 ; i < max_steps; ++i)
    {
        t = (t0 + t1) * 0.5f;
        p = ray(t).to_vec3();
        f_mid = obj.f( p );
        if(f_mid > target_iso){
            t1 = t;
            if((f_mid-target_iso) < eps) break;
        } else {
            t0 = t;
            if((target_iso-f_mid) < eps) break;
        }
    }
    return t;
}

Weitghed Laplacian smoothing

We use a simple smoothing scheme by iteratively replacing vertices to the barycenter of their 1st-ring neighbors. Computing the barycenter of the 1st ring neighbor can be computing with different barycentric weights. I recommand cotangent weights, they gave a good tradeoff between quality and simplicity of implementation.

You can find cotangent weights codes here.
Plug the cotangent weights in this simple smoothing procedure here.

Notes

A few things that are missing that I migth post later (mainly here as a reminder for me):

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: