cubic_ease.cpp

cubic_ease.cpp

float cubic_ease(float p1x, float p1y, float t, float p2x, float p2y) {
    
    if (t>=1) return 1;
    if (t<=0) return 0;
    
    
    static const double epsilon = 1e-7;
    auto cx_ = 3.0 * p1x;
    auto bx_ = 3.0 * (p2x - p1x) - cx_;
    auto ax_ = 1.0 - cx_ - bx_;
    
    auto cy_ = 3.0 * p1y;
    auto by_ = 3.0 * (p2y - p1y) - cy_;
    auto ay_ = 1.0 - cy_ - by_;
    
#define SampleCurveX(m_t) (((ax_ * m_t + bx_) * m_t + cx_) * m_t)
    // `ax t^3 + bx t^2 + cx t' expanded using Horner's rule.
#define SampleCurveY(m_t) (((ay_ * m_t + by_) * m_t + cy_) * m_t)
#define SampleCurveDerivativeX(m_t) ((3.0 * ax_ * m_t + 2.0 * bx_) * m_t + cx_)
    float t0, t1, t2, x2, d2;
    int i;
    
    auto x = t;
    
    // First try a few iterations of Newton's method -- normally very fast.
    for (t2 = x, i = 0; i < 8; i++) {
        
        x2 = SampleCurveX(t2) - x;
        if (fabsf(x2) < epsilon)
            return SampleCurveY(t2);
        d2 = SampleCurveDerivativeX(t2);
        if (fabsf(d2) < 1e-6)
            break;
        t2 = t2 - x2 / d2;
    }
    
    // Fall back to the bisection method for reliability.
    t0 = 0.0;
    t1 = 1.0;
    t2 = x;
    
    while (t0 < t1) {
        
        x2 = SampleCurveX(t2);
        if (fabsf(x2 - x) < epsilon)
            return SampleCurveY(t2);
        if (x > x2)
            t0 = t2;
        else
            t1 = t2;
        t2 = (t1 - t0) * .5 + t0;
    }
    
    // Failure.
    return SampleCurveY(t2);
}