From 669dd314fafdd1957303e06f5b6cb34d9daccaf8 Mon Sep 17 00:00:00 2001
From: hellerve <veit@veitheller.de>
Date: Mon, 5 Jun 2017 19:39:41 -0400
Subject: [PATCH] util: started implementing flong shaping functions

---
 utilities/flong.glsl | 94 ++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 94 insertions(+)
 create mode 100644 utilities/flong.glsl

diff --git a/utilities/flong.glsl b/utilities/flong.glsl
new file mode 100644
index 0000000..386af4b
--- /dev/null
+++ b/utilities/flong.glsl
@@ -0,0 +1,94 @@
+// these are shaping functions taken from various places of https://flong.com
+
+// Polynominal shaping
+
+float blinn_wyvill_cos_approx(float x) {
+  float x2 = x*x;
+  float x4 = x2*x2;
+  float x6 = x4*x2;
+
+  float fa = ( 4.0/9.0);
+  float fb = (17.0/9.0);
+  float fc = (22.0/9.0);
+
+  return fa*x6 - fb*x4 + fc*x2;
+}
+
+// TODO: factor out clamping
+float double_cubic_seat(float x, vec2 ab) {
+  float epsilon = 0.00001;
+  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;
+  float a = clamp(ab.x, min_param_a, max_param_a);
+  float b = clamp(ab.y, min_param_b, max_param_b);
+
+  if (x <= a){
+    return b - b*pow(1.-x/a, 3.0);
+  }
+  return b + (1.-b)*pow((x-a)/(1.-a), 3.0);
+}
+
+float double_cubic_seat_linear_blend(float x, vec2 ab) {
+  float epsilon = 0.00001;
+  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;
+  float a = clamp(ab.x, min_param_a, max_param_a);
+  float b = clamp(ab.y, min_param_b, max_param_b);
+  b = 1.0 - b;
+
+  if (x<=a){
+    return b*x + (1.-b)*a*(1.-pow(1.-x/a, 3.));
+  }
+  return b*x + (1.-b)*(a + (1.-a)*pow((x-a)/(1.-a), 3.));
+}
+
+float double_odd_polynomial_seat(float x, float a, float b, float n) {
+  float epsilon = 0.00001;
+  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 = clamp(a, min_param_a, max_param_a);
+  b = clamp(b, min_param_b, max_param_b);
+
+  float p = 2.*n + 1.;
+  if (x <= a){
+    return b - b*pow(1.-x/a, p);
+  }
+  return b + (1.-b)*pow((x-a)/(1.-a), p);
+}
+
+float double_poly_sigmoid(float x, float a, float b, float n) {
+  if(mod(n, 2.) == 0.) {
+    if(x <= 0.5) {
+      return pow(2.0*x, n)/2.0;
+    }
+    return 1.0 - pow(2.*(x-1.), n)/2.0;
+  } else {
+    if (x<=0.5){
+      return pow(2.0*x, n)/2.0;
+    }
+    return 1.0 + pow(2.0*(x-1.), n)/2.0;
+  }
+}
+
+float quad_through_point(float x, vec2 ab) {
+  float epsilon = 0.00001;
+  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;
+  float a = clamp(ab.x, min_param_a, max_param_a);
+  float b = clamp(ab.y, min_param_b, max_param_b);
+
+  float A = (1.-b)/(1.-a) - (b/a);
+  float B = (A*(a*a)-b)/a;
+  float y = A*(x*x) - B*(x);
+  y = clamp(0., 1., y);
+
+  return y;
+}