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:

Hrbf primitive

Once you computed an HRBF primitive (see source code HRBF Surface reconstruction code with Eigen lib) you need to be able to animate it and transform its global support to compact support. The code below is a simple wrapper allowing HRBF animation through a 4x4 matrix and conversion from global to local support.

/// @brief (-3/16)*(d/radius)^5 + (5/8)*(d/radius)^3 - (15/16)*(d/radius) + 1/2
/// junction is C2 at the boundary of the field [0, 1]
namespace  to_compact {

float f(float dist, float radius) {
    if(f < -radius)
        return 1.f;
    else if(f > radius)
        return 0.f;
    else {
        // (-3/16)*(f/radius)^5+(5/8)*(f/radius)^3-(15/16)*(f/radius) + 1/2
        const float fact   = (f/radius);
        const float fact2  = fact  * fact;
        const float fact4  = fact2 * fact2;
        return (-3.f/16.f)*fact4*fact + (5.f/8.f)*fact2*fact - (15.f/16.f)*fact + 1.f/2.f;
    }
}

void gf(float dist, float radius, Vec3& grad) {
    if(f < -r || f > r)
        grad = Vec3(0.f,0.f,0.f);
    else {
        // (-15/(16r))*(f/r)^4 + (15/(8*r))*(f/r)^2 - (15/(16r))
        const float fact  = f/r;
        const float fact2 = fact*fact;
        const float fact4 = fact2*fact2;
        const float tmp   = (-15.f/(16.f*r));
        const float scale = tmp*fact4 + (15.f/(8.f*r))*fact2 + tmp;
        grad = grad * scale;
    }
}

float fngf(float dist, float radius, Vec3& grad) {
    if(f < -radius) {
        grad = Vec3(0.f,0.f,0.f);
        return 1.f;
    }
    else if(f > radius) {
        grad = Vec3(0.f,0.f,0.f);
        return 0.f;
    } else {
        const float fact   = (f/radius);
        const float fact2  = fact  * fact;
        const float fact4  = fact2 * fact2;
        // Compute grad
        // (-15/(16r))*(f/r)^4 + (15/(8*r))*(f/r)^2 - (15/(16r))
        const float tmp   = (-15.f/(16.f*radius));
        const float scale = tmp*fact4 + (15.f/(8.f*radius))*fact2 + tmp;
        grad = grad * scale;
        // Compute potential
        // (-3/16)*(f/radius)^5+(5/8)*(f/radius)^3-(15/16)*(f/radius) + 1/2
        return (-3.f/16.f)*fact4*fact + (5.f/8.f)*fact2*fact - (15.f/16.f)*fact + 0.5f;
    }
}
}

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

/// Hermite radial basis function primitive with compact support
class Hrbf {
public:

    /// @param radius_support : distance from the 0.5f iso surface beyond which
    /// the field equals 0.f
    Hrbf(float radius_support = 10.f) : _support_radius(radius_support)
    {}

    void update(const std::vector& points,
                const std::vector& normals)
    {
        _hrbf_global.hermite_fit(points, normals);
    }

    float f(const Point3& pos) const
    {
        Point3 p = _global_frame.inverse() * pos;
        return to_compact::f(_hrbf_global.eval(p), _support_radius);
    }

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

    Vec3 gf (const Point3& pos) const
    {
        Point3 p =  _global_frame.inverse() * pos;
        float value =_hrbf_global.eval(p);
        Vec3 grad = _hrbf_global.grad(p);
        to_compact::gf(value, _support_radius, grad);
        return  _global_frame.inverse().transpose() * grad;
    }

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

    float fngf(const Vec3& pos, Vec3& grad) const
    {
        Point3 p =  _global_frame.inverse() * pos;
        float value = _hrbf_global.eval(p);
        grad = _hrbf_global.grad(p);
        float val = to_compact::fngf(value, _support_radius, grad);
        grad = _global_frame.inverse().transpose() * grad;
        return val;
    }

    /// Orientation and position of the Hrbf
    Matrix4x4 _global_frame;

private:
    /// Core of the HRBF math.
    /// Global support
    Hrbf_3d _hrbf_global;

    /// if( HRBF potential > _support_radius )
    ///     -> potential == 0
    float _support_radius;
};

Hrbf bounding box

Hrbf are very slow to evaluate, in order to accelerate implicit skinning it is necessary to cache the values of the hrbf into a 3d grid and evaluate its values trought tri-linear interpolation. To this end, one need to compute a bounding box around the hrbf to know which area should be discretized and sampled.

The code below is fairly ad hoc but should be robust enough for our usage. We take an HRBF with global support and compute the bounding box enclosing the iso surface specified by the user:

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_;
    }
};

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

struct Bbox3{

    Bbox3(): pmin(Point3( FLT_MAX,  FLT_MAX,  FLT_MAX)),
        pmax(Point3(-FLT_MAX, -FLT_MAX, -FLT_MAX))
    {     }

    void add_point(const Point3& p){
        pmin.x = fminf(p.x, pmin.x);
        pmin.y = fminf(p.y, pmin.y);
        pmin.z = fminf(p.z, pmin.z);
        pmax.x = fmaxf(p.x, pmax.x);
        pmax.y = fmaxf(p.y, pmax.y);
        pmax.z = fmaxf(p.z, pmax.z);
    }

    Vec3 lengths() const { return pmax-pmin; }

    void get_corners(std::vector& corners) const
    {
        corners.resize(8);
        for (int i = 0; i < 8; ++i)
            corners[i] = get_corner( i );
    }

    Point3 get_corner(int i) const
    {
        Vec3 diff = pmax - pmin;
        diff.x *= (i      & 0x1);
        diff.y *= (i >> 1 & 0x1);
        diff.z *= (i >> 2 & 0x1);
        return Point3( pmin.x + diff.x, pmin.y + diff.y, pmin.z + diff.z);
    }

    Point3 pmin;
    Point3 pmax;
};

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

struct Obbox{
    Obbox() : _transfo( Matrix4x4::identity() )
    {   }

    /// The local space bounding box (a normal axis bouding box)
    Bbox3 _local_bbox;
    /// The associated transformation to convert the local bbox
    /// to its global orientation
    Matrix4x4 _transfo;
};

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

const int GRID_RES = 8;

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

float dichotomic_search(const Ray& r,
                        float t0, float t1,
                        float target_iso,
                        const Hrbf_global& hrbf,
                        float eps = 0.01f)
{
    float t = t0;
    float f0 = hrbf.eval( r(t0).to_vec3() );
    float f1 = hrbf.eval( r(t1).to_vec3() );

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

    Point3 p;
    for(unsigned short i = 0 ; i < 15; ++i)
    {
        t = (t0 + t1) * 0.5f;
        p = r(t);

        f0 = hrbf.eval(p.to_vec3());

        if(f0 > target_iso){
            t1 = t;
            if((f0-target_iso) < eps) break;
        } else {
            t0 = t;
            if((target_iso-f0) < eps) break;
        }
    }
    return t;
}

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

/// Cast a ray and return the farthest point matching 'target_iso'.
/// We use newton iterations
/// @param start : origin of the ray must be inside the primitive
/// @param hrbf : field function we seek the iso-surface 'target-iso'
/// @param target_iso : iso_surface we want to reach while pushing the point start
/// @param hrbf_scale : overall scale of the hrbf.
/// Will be used for error detections
/// @param custom_dir : if true we don't follow the gradient but use 'dir'
/// defined by the user
/// @return the farthest point matching 'target_iso' along the ray.
Point3 push_point(const Point3& start,
                  const Hrbf_global& hrbf,
                  float target_iso,
                  float hrbf_scale,
                  bool custom_dir = false,
                  const Vec3 dir = Vec3())
{
    Vec3  grad;
    Vec3  n_dir = dir.normalized();
    Vec3  step;
    Point3 res  = start;
    Point3 prev = res;

    for(int i = 0; i < 20; i++)
    {
        grad = hrbf.grad(res);
        const float pot      = hrbf.eval(res);
        const float pot_diff = fabsf( pot - target_iso);
        const float norm     = grad.safe_normalize();

        // Check the norm are large enough
        if( custom_dir && n_dir.norm() > 0.00001f)
            step = n_dir;
        else if( norm > 0.00001f )
            step = grad;
        else{
            // Sometime moving a bit the evaluted point will help falling back
            // in an area where the gradient is not null
            res += Vec3(0.00001f * hrbf_scale, 0.f, 0.f);
            continue;
        }

        float scale = (pot_diff / norm) * 0.4f;

        // check our step is not large compared to the overall scale of the object
        if( scale >  hrbf_scale * 4.f * target_iso ) {
            // Move a bit the point see if we get somewhere better...
            res += Vec3(0.00001f * hrbf_scale, 0.f, 0.f);
            continue;
        }

        if( pot > target_iso)
        {
            Ray r( prev, step);
            float t = dichotomic_search(r, 0.f, (res-prev).norm(), target_iso, hrbf);
            res = r( t );
            break;
        }

        prev = res;
        res = res + step * scale;

        if( pot_diff <= 0.01f || scale < 0.00001f)
            break;
    }

    return res;
}

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

/// Cast several rays along a square grid and return the farthest point which
/// potential is zero.
/// @param obbox : oriented bounding box which we add casted points to
/// @param org : origin point of the grid (top left corner)
/// @param x : extrimity point of the grid in x direction (top right corner)
/// @param y : extrimity point of the grid in y direction (bottom left corner)
/// @param hrbf_scale : overall scale of the hrbf.
/// @warning org, x and y parameters must lie in the same plane. Also this
/// (org-x) dot (org-y) == 0 must be true at all time. Otherwise the behavior
/// is undefined.
void push_face(Obbox& obbox,
               const Point3& org,
               const Point3& x,
               const Point3& y,
               const Hrbf_global& hrbf,
               float target_iso,
               float hrbf_scale)
{
    int res = GRID_RES;

    Vec3 axis_x = x-org;
    Vec3 axis_y = y-org;
    Vec3 udir = (axis_x.cross(axis_y)).normalized();

    float len_x = axis_x.normalize();
    float len_y = axis_y.normalize();

    float step_x = len_x / (float)res;
    float step_y = len_y / (float)res;

    Matrix4x4 tr = obbox._transfo.fast_inverse();
    for(int i = 0; i < res; i++)
    {
        for(int j = 0; j < res; j++)
        {
            Point3 p = org +
                       axis_x * (step_x * (float)i + step_x * 0.5f) +
                       axis_y * (step_y * (float)j + step_y * 0.5f);

            Point3 res = push_point(p, hrbf, target_iso, hrbf_scale, true, udir);

            obbox._local_bbox.add_point( tr * res);
        }
    }
}

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

Obbox compute_obbox( const Hrbf_global& hrbf, float target_iso, const Matrix4x4& orientation ) const
{
    Obbox obbox;
    obbox._transfo = orientation;
    Matrix4x4 bbox_tr_inv = obbox._transfo.fast_inverse();

    // No samples, no obbox then we choose default values
    if(_points.size() == 0)
    {
        obbox._local_bbox.pmin = Point3(-1.f);
        obbox._local_bbox.pmax = Point3( 1.f);
        return obbox;
    }

    // Seek target_iso along samples normals of the HRBF
    for(unsigned i = 0; i < hrbf._points.size(); i++)
        obbox._local_bbox.add_point(bbox_tr_inv * Point3(hrbf._points[i]));

    float bbox_scale = obbox._local_bbox.lengths().get_max();
    for(unsigned i = 0; i < _points.size(); i++){
        Point3 pt = push_point(Point3(_points[i]), hrbf, target_iso, bbox_scale);
        obbox._local_bbox.add_point(bbox_tr_inv * pt);
    }

    // Push obbox faces by discretizing them into regulars 2d grid.
    // Each point of the grid is pushed to the target iso
    std::vector corners;
    /**
        @code

            6 +----+ 7
             /|   /|
          2 +----+3|
            |4+--|-+5
            |/   |/
            +----+
           0      1
        // Vertex 0 is pmin and vertex 7 pmax
        @endcode
    */
    for (int i = 0; i < 2; ++i)
    {
        // Get obox in world coordinates
        obbox._local_bbox.get_corners(corners);
        for(int i = 0; i < 8; ++i)
            corners[i] = obbox._transfo * corners[i];


        Point3 list[6][3] =   {{corners[2], corners[3], corners[0]}, // FRONT
                               {corners[0], corners[1], corners[4]}, // BOTTOM
                               {corners[3], corners[7], corners[1]}, // RIGHT
                               {corners[6], corners[2], corners[4]}, // LEFT
                               {corners[7], corners[6], corners[5]}, // REAR
                               {corners[6], corners[7], corners[2]}, // TOP
                              };

        // Pushing according a grid from the box's faces
        for(int i = 0; i < 6; ++i) {
            push_face(obbox, list[i][0], list[i][1], list[i][2],
                    hrbf, target_iso, bbox_scale);
        }
    }

    return obbox;
}

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;
}

Weighted 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: