ShapingFunctions.glsl

// Collection of Open Source Shaping Functions in GLSL
// https://thebookofshaders.com/05/
// https://www.shadertoy.com/view/3lVyDD
//============================================================================================================

#define PI 3.14159265359
#define epsilon 0.00001

float sq(float x) {
    return x*x;
}


// Polynomial shaping functions
// http://www.flong.com/archive/texts/code/shapers_poly/index.html
//------------------------------------------------------------------------------------------------------------

// Blinn Wyvill raised inverted cosine approximation
float bwCos (float x) {
    float x2 = x*x;
    float x4 = x2*x2;
    float x6 = x4*x2;
    
    float fa = 4./9.;
    float fb = 17./9.;
    float fc = 22./9.;
    
    return fa*x6 - fb*x4 + fc*x2;
}

// Double cubic seat
// joining two cubic curves at a horizontal inflection point at coord (a,b) in the unit square
float dcSeat (float x, float a, float b) {
  float min_param_a = 0.0 + epsilon;
  float max_param_a = 1.0 - epsilon;
  float min_param_b = 0.0;
  float max_param_b = 1.0;
  a = min(max_param_a, max(min_param_a, a));  
  b = min(max_param_b, max(min_param_b, b)); 
  
  float y = 0.;
  if (x <= a){
    y = b - b*pow(1.-x/a, 3.0);
  } else {
    y = b + (1.-b)*pow((x-a)/(1.-a), 3.0);
  }
  return y;
}

// Double cubic seat with linear blend
// a controls inflection point across the diagonal, b controls blending with identity function y = x
float dcSeatBlend (float x, float a, float b) {
  float min_param_a = 0.0 + epsilon;
  float max_param_a = 1.0 - epsilon;
  float min_param_b = 0.0;
  float max_param_b = 1.0;
  a = min(max_param_a, max(min_param_a, a));  
  b = min(max_param_b, max(min_param_b, b)); 
  b = 1.0 - b; //reverse for intelligibility.
  
  float y = 0.;
  if (x<=a){
    y = b*x + (1.-b)*a*(1.-pow(1.-x/a, 3.0));
  } else {
    y = b*x + (1.-b)*(a + (1.-a)*pow((x-a)/(1.-a), 3.0));
  }
  return y;
}

// Double Odd Polynomial seat
// uses paramater n (between 1 to ~20) to determine the flatness of the inflection point
float dcSeatOdd(float x, float a, float b, float n) {
  float min_param_a = 0.0 + epsilon;
  float max_param_a = 1.0 - epsilon;
  float min_param_b = 0.0;
  float max_param_b = 1.0;
  a = min(max_param_a, max(min_param_a, a));  
  b = min(max_param_b, max(min_param_b, b)); 
  
  float p = 2.*n + 1.;
  float y = 0.;
  if (x <= a)
      y = b - b*pow(1.-x/a, p);
  else
      y = b + (1.-b)*pow((x-a)/(1.-a), p);
     
  return y;
}

// Symmetric Double Polynomial Sigmoids
// a,b inflection point, n (int from 0 to 10) the steepness of the wall in between.
// approximates raised inverted cos to within 2.8%
float doublePolySigmoid (float x, float a, float b, float n) {
  float y = 0.;
  if (mod(n,2.) == 0.){ 
    // even polynomial
    if (x <= 0.5)
      y = pow(2.0 * x, n) / 2.0;
    else
      y = 1.0 - pow(2.*(x-1.), n)/2.0;
  } 
  else { 
    // odd polynomial
    if (x <= 0.5)
      y = pow(2.0*x, n)/2.0;
    else
      y = 1.0 + pow(2.0*(x-1.), n)/2.0;
  }

  return y;
}

// Quadratic through a point
// not all combinations of (a,b) will return a correct quad
float quadAtPoint(float x, float a, float b) {
  float min_param_a = 0.0 + epsilon;
  float max_param_a = 1.0 - epsilon;
  float min_param_b = 0.0;
  float max_param_b = 1.0;
  a = min(max_param_a, max(min_param_a, a));  
  b = min(max_param_b, max(min_param_b, b)); 
  
  float A = (1.-b)/(1.-a) - (b/a);
  float B = (A*(a*a)-b)/a;
  float y = A*(x*x) - B*(x);
  y = min(1.,max(0.,y)); 
  
  return y;
}








// Exponential shaping functions
// http://www.flong.com/archive/texts/code/shapers_exp/index.html
//------------------------------------------------------------------------------------------------------------

// Exponential Ease-in and Ease-out
// a varies the function from an ease out form to an ease in form
float exponentialEase(float x, float a) {
    float min_param_a = 0. + epsilon;
    float max_param_a = 1. - epsilon;
    a = max(min_param_a, min(max_param_a, a));
    
    if (a < 0.5)
        return pow(x, 2.*a);
    else
        return pow(x, 1./(1.-(2.*(a-.5))));
}

// Double Exponential Seat
// better derivatives than cubic, a from (0, 1)
float doubleExpSeat(float x, float a) {
  float min_param_a = 0.0 + epsilon;
  float max_param_a = 1.0 - epsilon;
  a = min(max_param_a, max(min_param_a, a)); 

  if (x <= 0.5)
    return (pow(2.0*x, 1.-a))/2.0;
  else
    return 1.0 - (pow(2.0*(1.0-x), 1.-a))/2.0;
}

// Double Exponential Sigmoid
// useful for adjustable-contrast functions, approximates raised inverted cosine to within 1% with a = 0.426
float doubleExpSigmoid(float x, float a) {
  float min_param_a = 0.0 + epsilon;
  float max_param_a = 1.0 - epsilon;
  a = min(max_param_a, max(min_param_a, a));
  a = 1.0-a; // for sensible results
  
  if (x <= 0.5)
    return (pow(2.0*x, 1.0/a))/2.0;
  else
    return 1.0 - (pow(2.0*(1.0-x), 1.0/a))/2.0;
}

// Logistic Sigmoid
// can represent the growth of organic populations/ other natural phenomena. Often used for weighting 
// signal-response functions in neural networks. param a regulates the slope/ growth rate of the sigmoid
// during its rise. a = 0 leads to y = x. Expensive to calculate. 
float logisticSigmoid(float x, float a) {
  float min_param_a = 0.0 + epsilon;
  float max_param_a = 1.0 - epsilon;
  a = max(min_param_a, min(max_param_a, a));
  a = (1./(1.-a) - 1.);

  float A = 1.0 / (1.0 + exp(0. -((x-0.5)*a*2.0)));
  float B = 1.0 / (1.0 + exp(a));
  float C = 1.0 / (1.0 + exp(0.-a)); 
  float y = (A-B)/(C-B);
  return y;
}








// Circular and Elliptical Shaping Functions
// http://www.flong.com/archive/texts/code/shapers_circ/index.html
//------------------------------------------------------------------------------------------------------------

float circularEaseIn (float x){
  float y = 1. - sqrt(1. - x*x);
  return y;
}

float circularEaseOut (float x){
  float y = sqrt(1. - sq(1. - x));
  return y;
}

// Double Circle Seat
// param a [0, 1] governs the location of the horizontal tangent
float doubleCircleSeat(float x, float a) {
    float min_param_a = 0.;
    float max_param_a = 1.;
    a = max(min_param_a, min(max_param_a, a));
    
    if (x <= a)
        return sqrt(sq(a) - sq(x - a));
    else
        return 1. - sqrt(sq(1.-a) - sq(x-a));
}

// Double Circle Sigmoid
// two circular arts meet at a vertical tangent 
float doubleCircleSigmoid(float x, float a) {
    float min_param_a = 0.0;
    float max_param_a = 1.0;
    a = max(min_param_a, min(max_param_a, a));
    
    if (x <= a) 
        return a - sqrt(a*a - x*x);
    else
        return a + sqrt(sq(1.-a) - sq(x-1.));
}


// Double Elliptic Seat
// double circle seat meeting at coord (a,b)
float doubleEllipticSeat(float x, float a, float b) {
  float min_param_a = 0.0 + epsilon;
  float max_param_a = 1.0 - epsilon;
  float min_param_b = 0.0;
  float max_param_b = 1.0;
  a = max(min_param_a, min(max_param_a, a)); 
  b = max(min_param_b, min(max_param_b, b)); 

  float y = 0.;
  if (x<=a){
    y = (b/a) * sqrt(sq(a) - sq(x-a));
  } else {
    y = 1. - ((1.-b)/(1.-a))*sqrt(sq(1.-a) - sq(x-a));
  }
  return y;
}

// Double Elliptic Sigmoid
// double circle sigmoid meeting at coord (a,b)
float doubleEllipticSigmoid(float x, float a, float b) {
  float min_param_a = 0.0 + epsilon;
  float max_param_a = 1.0 - epsilon;
  float min_param_b = 0.0;
  float max_param_b = 1.0;
  a = max(min_param_a, min(max_param_a, a)); 
  b = max(min_param_b, min(max_param_b, b));
 
  float y = 0.;
  if (x<=a){
    y = b * (1. - (sqrt(sq(a) - sq(x))/a));
  } else {
    y = b + ((1.-b)/(1.-a))*sqrt(sq(1.-a) - sq(x-1.));
  }
  return y;
}

// Double-Linear with Circular Fillet
// joins two straight lines with a circular arc whose radius is adjustable. param R is the radius, 
// coord of intersection is (a, b). Adapted from Robert D. Miller in Graphics Gems III
float arcStartAngle;
float arcEndAngle;
float arcStartX,  arcStartY;
float arcEndX,    arcEndY;
float arcCenterX, arcCenterY;
float arcRadius;
void computeFilletParameters (
  float p1x, float p1y, 
  float p2x, float p2y, 
  float p3x, float p3y, 
  float p4x, float p4y,
  float r);
float linetopoint (float a, float b, float c, float ptx, float pty);

//--------------------------------------------------------
float circularFillet (float x, float a, float b, float R){
  
  float min_param_a = 0.0 + epsilon;
  float max_param_a = 1.0 - epsilon;
  float min_param_b = 0.0 + epsilon;
  float max_param_b = 1.0 - epsilon;
  a = max(min_param_a, min(max_param_a, a)); 
  b = max(min_param_b, min(max_param_b, b)); 

  computeFilletParameters (0., 0., a, b, a, b, 1., 1.,  R);
  float t = 0.;
  float y = 0.;
  x = max(0., min(1., x));
  
  if (x <= arcStartX){
    t = x / arcStartX;
    y = t * arcStartY;
  } else if (x >= arcEndX){
    t = (x - arcEndX)/(1. - arcEndX);
    y = arcEndY + t*(1. - arcEndY);
  } else {
    if (x >= arcCenterX){
      y = arcCenterY - sqrt(sq(arcRadius) - sq(x-arcCenterX)); 
    } else{
      y = arcCenterY + sqrt(sq(arcRadius) - sq(x-arcCenterX)); 
    }
  }
  return y;
}

//------------------------------------------
// Return signed distance from line Ax + By + C = 0 to point P.
float linetopoint (float a, float b, float c, float ptx, float pty){
  float lp = 0.0;
  float d = sqrt((a*a)+(b*b));
  if (d != 0.0){
    lp = (a*ptx + b*pty + c)/d;
  }
  return lp;
}

//------------------------------------------
// Compute the parameters of a circular arc 
// fillet between lines L1 (p1 to p2) and
// L2 (p3 to p4) with radius R.  
// 
void computeFilletParameters (
  float p1x, float p1y, 
  float p2x, float p2y, 
  float p3x, float p3y, 
  float p4x, float p4y,
  float r) {

  float c1   = p2x*p1y - p1x*p2y;
  float a1   = p2y-p1y;
  float b1   = p1x-p2x;
  float c2   = p4x*p3y - p3x*p4y;
  float a2   = p4y-p3y;
  float b2   = p3x-p4x;
  if ((a1*b2) == (a2*b1)) {  /* Parallel or coincident lines */
    return;
  }

  float d1, d2;
  float mPx, mPy;
  mPx = (p3x + p4x)/2.0;
  mPy = (p3y + p4y)/2.0;
  d1 = linetopoint(a1,b1,c1,mPx,mPy);  /* Find distance p1p2 to p3 */
  if (d1 == 0.0) {
    return; 
  }
  mPx = (p1x + p2x)/2.0;
  mPy = (p1y + p2y)/2.0;
  d2 = linetopoint(a2,b2,c2,mPx,mPy);  /* Find distance p3p4 to p2 */
  if (d2 == 0.0) {
    return; 
  }

  float c1p, c2p, d;
  float rr = r;
  if (d1 <= 0.0) {
    rr= -rr;
  }
  c1p = c1 - rr*sqrt((a1*a1)+(b1*b1));  /* Line parallel l1 at d */
  rr = r;
  if (d2 <= 0.0) {
    rr = -rr;
  }
  c2p = c2 - rr*sqrt((a2*a2)+(b2*b2));  /* Line parallel l2 at d */
  d = (a1*b2)-(a2*b1);

  float pCx = (c2p*b1-c1p*b2)/d; /* Intersect constructed lines */
  float pCy = (c1p*a2-c2p*a1)/d; /* to find center of arc */
  float pAx = 0.;
  float pAy = 0.;
  float pBx = 0.;
  float pBy = 0.;
  float dP,cP;

  dP = (a1*a1) + (b1*b1);        /* Clip or extend lines as required */
  if (dP != 0.0) {
    cP = a1*pCy - b1*pCx;
    pAx = (-a1*c1 - b1*cP)/dP;
    pAy = ( a1*cP - b1*c1)/dP;
  }
  dP = (a2*a2) + (b2*b2);
  if (dP != 0.0) {
    cP = a2*pCy - b2*pCx;
    pBx = (-a2*c2 - b2*cP)/dP;
    pBy = ( a2*cP - b2*c2)/dP;
  }

  float gv1x = pAx-pCx; 
  float gv1y = pAy-pCy;
  float gv2x = pBx-pCx; 
  float gv2y = pBy-pCy;

  float arcStart = atan(gv1y, gv1x); // may need to swap these
  float arcAngle = 0.0;
  float dd = sqrt(((gv1x*gv1x)+(gv1y*gv1y)) * ((gv2x*gv2x)+(gv2y*gv2y)));
  if (dd != 0.0) {
    arcAngle = (acos((gv1x*gv2x + gv1y*gv2y)/dd));
  } 
  float crossProduct = (gv1x*gv2y - gv2x*gv1y);
  if (crossProduct < 0.0) { 
    arcStart -= arcAngle;
  }

  float arc1 = arcStart;
  float arc2 = arcStart + arcAngle;
  if (crossProduct < 0.0) {
    arc1 = arcStart + arcAngle;
    arc2 = arcStart;
  }

  arcCenterX    = pCx;
  arcCenterY    = pCy;
  arcStartAngle = arc1;
  arcEndAngle   = arc2;
  arcRadius     = r;
  arcStartX     = arcCenterX + arcRadius*cos(arcStartAngle);
  arcStartY     = arcCenterY + arcRadius*sin(arcStartAngle);
  arcEndX       = arcCenterX + arcRadius*cos(arcEndAngle);
  arcEndY       = arcCenterY + arcRadius*sin(arcEndAngle);
}

// Circular Arc through a point
// not every point in a unit square works, works best around the diaganol. Adapted from Paul Bourke's
// Equation of a Circle From 3 Points
float m_Centerx;
float m_Centery;
float m_dRadius;
void calcCircleFrom3Points(float pt1x, float pt1y, float pt2x, float pt2y, float pt3x, float pt3y);
bool IsPerpendicular(float pt1x, float pt1y, float pt2x, float pt2y, float pt3x, float pt3y);
void calcCircleFrom3Points(float pt1x, float pt1y, float pt2x, float pt2y, float pt3x, float pt3y);


float circularArcAtPoint (float x, float a, float b){  
  float min_param_a = 0.0 + epsilon;
  float max_param_a = 1.0 - epsilon;
  float min_param_b = 0.0 + epsilon;
  float max_param_b = 1.0 - epsilon;
  a = min(max_param_a, max(min_param_a, a));
  b = min(max_param_b, max(min_param_b, b));
  x = min(1.0-epsilon, max(0.0+epsilon, x));
  
  float pt1x = 0.;
  float pt1y = 0.;
  float pt2x = a;
  float pt2y = b;
  float pt3x = 1.;
  float pt3y = 1.;

  if (!IsPerpendicular(pt1x,pt1y, pt2x,pt2y, pt3x,pt3y))		
     calcCircleFrom3Points (pt1x,pt1y, pt2x,pt2y, pt3x,pt3y);	
  else if (!IsPerpendicular(pt1x,pt1y, pt3x,pt3y, pt2x,pt2y))		
     calcCircleFrom3Points (pt1x,pt1y, pt3x,pt3y, pt2x,pt2y);	
  else if (!IsPerpendicular(pt2x,pt2y, pt1x,pt1y, pt3x,pt3y))		
     calcCircleFrom3Points (pt2x,pt2y, pt1x,pt1y, pt3x,pt3y);	
  else if (!IsPerpendicular(pt2x,pt2y, pt3x,pt3y, pt1x,pt1y))		
     calcCircleFrom3Points (pt2x,pt2y, pt3x,pt3y, pt1x,pt1y);	
  else if (!IsPerpendicular(pt3x,pt3y, pt2x,pt2y, pt1x,pt1y))		
     calcCircleFrom3Points (pt3x,pt3y, pt2x,pt2y, pt1x,pt1y);	
  else if (!IsPerpendicular(pt3x,pt3y, pt1x,pt1y, pt2x,pt2y))		
     calcCircleFrom3Points (pt3x,pt3y, pt1x,pt1y, pt2x,pt2y);	
  else { 
    return 0.;
  }

  // constrain
  if ((m_Centerx > 0.) && (m_Centerx < 1.)) {
     if (a < m_Centerx) {
       m_Centerx = 1.;
       m_Centery = 0.;
       m_dRadius = 1.;
     } else {
       m_Centerx = 0.;
       m_Centery = 1.;
       m_dRadius = 1.;
     }
  }
  
  float y = 0.;
  if (x >= m_Centerx) {
    y = m_Centery - sqrt(sq(m_dRadius) - sq(x-m_Centerx)); 
  } else {
    y = m_Centery + sqrt(sq(m_dRadius) - sq(x-m_Centerx)); 
  }
  return y;
}

//----------------------
bool IsPerpendicular(float pt1x, float pt1y, float pt2x, float pt2y, float pt3x, float pt3y)
{
  // Check the given point are perpendicular to x or y axis 
  float yDelta_a = pt2y - pt1y;
  float xDelta_a = pt2x - pt1x;
  float yDelta_b = pt3y - pt2y;
  float xDelta_b = pt3x - pt2x;
  float epsilon2 = 0.000001;

  // checking whether the line of the two pts are vertical
  if (abs(xDelta_a) <= epsilon2 && abs(yDelta_b) <= epsilon2){
    return false;
  }
  if (abs(yDelta_a) <= epsilon2) {
    return true;
  }
  else if (abs(yDelta_b) <= epsilon2) {
    return true;
  }
  else if (abs(xDelta_a)<= epsilon2) {
    return true;
  }
  else if (abs(xDelta_b)<= epsilon2) {
    return true;
  }
  else return false;
}

//--------------------------
void calcCircleFrom3Points(float pt1x, float pt1y, float pt2x, float pt2y, float pt3x, float pt3y)
{
    float yDelta_a = pt2y - pt1y;
    float xDelta_a = pt2x - pt1x;
    float yDelta_b = pt3y - pt2y;
    float xDelta_b = pt3x - pt2x;
    float epsilon2 = 0.000001;

    if (abs(xDelta_a) <= epsilon2 && abs(yDelta_b) <= epsilon2){
        m_Centerx = 0.5*(pt2x + pt3x);
        m_Centery = 0.5*(pt1y + pt2y);
        m_dRadius = sqrt(sq(m_Centerx-pt1x) + sq(m_Centery-pt1y));
        return;
    }

    // IsPerpendicular() assure that xDelta(s) are not zero
    float aSlope = yDelta_a / xDelta_a; 
    float bSlope = yDelta_b / xDelta_b;
    if (abs(aSlope-bSlope) <= epsilon2) {	
        // checking whether the given points are colinear. 	
        return;
    }

    // calc center
    m_Centerx = (
        aSlope*bSlope*(pt1y - pt3y) + 
        bSlope*(pt1x + pt2x) - 
        aSlope*(pt2x+pt3x) )
        /(2.* (bSlope-aSlope));
    m_Centery = -1.*(m_Centerx - (pt1x+pt2x)/2.)/aSlope +  (pt1y+pt2y)/2.;
    m_dRadius = sqrt(sq(m_Centerx-pt1x) + sq(m_Centery-pt1y));
}








// Bezier and Other Parametric Shaping Functions
// http://www.flong.com/archive/texts/code/shapers_bez/index.html
//------------------------------------------------------------------------------------------------------------

// Quadratic Bezier
// 2nd order bezier with a single spline control point at (a,b). Same rates of change/ slope at its endpoints
float quadraticBezier(float x, float a, float b) {
    a = max(0., min(1., a)); 
    b = max(0., min(1., b)); 
    if (a == 0.5){
        a += epsilon;
    }
    // solve t from x (an inverse operation)
    float om2a = 1. - 2.*a;
    float t = (sqrt(a*a + om2a*x) - a)/om2a;
    float y = (1.-2.*b)*(t*t) + (2.*b)*t;
    return y;   
}

// Cubic Bezier
float slopeFromT (float t, float A, float B, float C);
float xFromT (float t, float A, float B, float C, float D);
float yFromT (float t, float E, float F, float G, float H);

float cubicBezier (float x, float a, float b, float c, float d){
    float y0a = 0.00; // initial y
    float x0a = 0.00; // initial x 
    float y1a = b;    // 1st influence y   
    float x1a = a;    // 1st influence x 
    float y2a = d;    // 2nd influence y
    float x2a = c;    // 2nd influence x
    float y3a = 1.00; // final y 
    float x3a = 1.00; // final x 

    float A =   x3a - 3.*x2a + 3.*x1a - x0a;
    float B = 3.*x2a - 6.*x1a + 3.*x0a;
    float C = 3.*x1a - 3.*x0a;   
    float D =   x0a;

    float E =   y3a - 3.*y2a + 3.*y1a - y0a;    
    float F = 3.*y2a - 6.*y1a + 3.*y0a;             
    float G = 3.*y1a - 3.*y0a;             
    float H =   y0a;

    // Solve for t given x (using Newton-Raphelson), then solve for y given t.
    // Assume for the first guess that t = x.
    float currentt = x;
    int nRefinementIterations = 5;
    for (int i=0; i < nRefinementIterations; i++){
        float currentx = xFromT (currentt, A,B,C,D); 
        float currentslope = slopeFromT (currentt, A,B,C);
        currentt -= (currentx - x)*(currentslope);
        currentt = clamp(currentt, 0., 1.);
    } 

    float y = yFromT (currentt,  E,F,G,H);
    return y;
}

// Helper functions:
float slopeFromT (float t, float A, float B, float C) {
    float dtdx = 1.0/(3.0*A*t*t + 2.0*B*t + C); 
    return dtdx;
}

float xFromT (float t, float A, float B, float C, float D) {
    float x = A*(t*t*t) + B*(t*t) + C*t + D;
    return x;
}

float yFromT (float t, float E, float F, float G, float H) {
    float y = E*(t*t*t) + F*(t*t) + G*t + H;
    return y;
}

// Cubic Bezier (nearly) at two points

// Helper functions. 
float B0 (float t){
    return (1.-t)*(1.-t)*(1.-t);
}
float B1 (float t){
    return  3.*t* (1.-t)*(1.-t);
}
float B2 (float t){
    return 3.*t*t* (1.-t);
}
float B3 (float t){
    return t*t*t;
}
float  findx (float t, float x0, float x1, float x2, float x3){
    return x0*B0(t) + x1*B1(t) + x2*B2(t) + x3*B3(t);
}
float  findy (float t, float y0, float y1, float y2, float y3){
    return y0*B0(t) + y1*B1(t) + y2*B2(t) + y3*B3(t);
}

float cubicBezierAt(float x, float a, float b, float c, float d) {
    float y = 0.;
    float min_param_a = 0.0 + epsilon;
    float max_param_a = 1.0 - epsilon;
    float min_param_b = 0.0 + epsilon;
    float max_param_b = 1.0 - epsilon;
    a = max(min_param_a, min(max_param_a, a));
    b = max(min_param_b, min(max_param_b, b));

    float x0 = 0.;  
    float y0 = 0.;
    float x4 = a;  
    float y4 = b;
    float x5 = c;  
    float y5 = d;
    float x3 = 1.;  
    float y3 = 1.;
    float x1,y1,x2,y2; // to be solved.

    // arbitrary but reasonable 
    // t-values for interior control points
    float t1 = 0.3;
    float t2 = 0.7;

    float B0t1 = B0(t1);
    float B1t1 = B1(t1);
    float B2t1 = B2(t1);
    float B3t1 = B3(t1);
    float B0t2 = B0(t2);
    float B1t2 = B1(t2);
    float B2t2 = B2(t2);
    float B3t2 = B3(t2);

    float ccx = x4 - x0*B0t1 - x3*B3t1;
    float ccy = y4 - y0*B0t1 - y3*B3t1;
    float ffx = x5 - x0*B0t2 - x3*B3t2;
    float ffy = y5 - y0*B0t2 - y3*B3t2;

    x2 = (ccx - (ffx*B1t1)/B1t2) / (B2t1 - (B1t1*B2t2)/B1t2);
    y2 = (ccy - (ffy*B1t1)/B1t2) / (B2t1 - (B1t1*B2t2)/B1t2);
    x1 = (ccx - x2*B2t1) / B1t1;
    y1 = (ccy - y2*B2t1) / B1t1;

    x1 = max(0.+epsilon, min(1.-epsilon, x1));
    x2 = max(0.+epsilon, min(1.-epsilon, x2));

    // Note that this function also requires cubicBezier()!
    y = cubicBezier (x, x1,y1, x2,y2);
    y = max(0., min(1., y));
    return y;
}





// Inigo Quilez Useful Functions
// https://www.iquilezles.org/www/articles/functions/functions.htm

// Almost Identity
// smoothly blend a value against a lowerbound, anything above m remains unchanged,
// n is the value to be taken when input is zero
float almostIdentity(float x, float m, float n) {
    if (x > m) return x;
    float a = 2.*n - m;
    float b = 2.*m - 3.*n;
    float t = x/m;
    return (a*t + b)*t*t + n;
}

// Almost Unit Identity
// maps the unit interval, 0 deriv at origin and 1 deriv at 1. Ideal for transitioning from stationary to motion,
// equivalent to almostIdentity(x, 1, 0)
float almostIdentity(float x) {
    return x*x*(2.-x);
}

// Almost Identity II
// not as fast as cubics above, can also be used as a smooth-abs(), has a zero derivative, but has a non zero second
// derivative
float almostIdentity(float x, float n) {
    return sqrt(x*x+n);
}

// Exponential Impulse
// great for triggering behaviours or envelopes for music or animation or anything that grows fast then slowly decays.
// use k to control the stretching of the function. its maximum of 1 is at x = 1/k
float expImpulse(float x, float k) {
    float h = k*x;
    return h*exp(1.-h);
}

// Sustained Impulse
float expSustainImpulse(float x, float f, float k) {
    float s = max(x - f, .0);
    return min(x*x/ (f*f), 1. + (2./f) * s * exp(-k*s));
}

// Polynomial Impulse
// peaks at x = sqrt(1/k)
float quaImpulse(float k, float x) {
    return 2.0*sqrt(k)*x/(1.0+k*x*x);
}
// n is the degree of the polynomial 
float polyImpulse(float k, float n, float x) {
    return (n/(n-1.0))*pow((n-1.0)*k,1.0/n) * x/(1.0+k*pow(x,n));
}

// Cubic Pulse
// smoothstep(c-w,c,x)-smoothstep(c,c+w,x) replaces this since it isolates some features in a signal, 
// replaces a gaussian as well
float cubicPulse(float c, float w, float x) {
    x = abs(x - c);
    if(x>w) return 0.0;
    x /= w;
    return 1.0 - x*x*(3.0-2.0*x);
}

// Exponential Step
// the higher the power n, the sharper s shaped curve you get until you get step()
float expStep(float x, float k, float n) {
    return exp(-k*pow(x,n));
}

// Gain
// remaps the unit interval by expanding the sides and compressing the center (0.5 maps to 0.5). k=1 is the identity curve, k<1
// produces classic gain shape, k>1 produces s curves. symmetric and inverse for k=a and k=1/a
float gain(float x, float k) {
    float a = 0.5*pow(2.0*((x<0.5)?x:1.0-x), k);
    return (x<0.5)?a:1.0-a;
}

// Parabola
// maps x = 0, x = 1 to 0 and x = 0.5 to 1
float parabola(float x, float k) {
    return pow(4.0*x*(1.0-x), k);
}

// Power Curve
// generalization of parabola, except a controls shape of left side of curve and b controls shape of the right side.
// slow calculation of k ensures the curve's max is at 1
float pcurve_uniform(float x, float a, float b) {
    float k = pow(a+b,a+b) / (pow(a,a)*pow(b,b));
    return k * pow(x,a) * pow(1.0-x, b);
}

// Power Curve
float pcurve(float x, float a, float b) {
    return pow(x,a) * pow(1.0-x, b);
}

// Sinc Curve
// phase shifted sinc curve, can be useful for bouncing behaviours, k increases num of bounces
float sinc(float x, float k) {
    float a = PI*(k*x-1.0);
    return sin(a)/a;
}