curve: de casteljau implementation for solving cubic bezier

CREDITS

@@ -55,3 +55,9 @@ https://github.com/erich666/GraphicsGems/blob/master/gems/BoxSphere.c
 
 10. Horizontal add
 https://stackoverflow.com/questions/6996764/fastest-way-to-do-horizontal-float-vector-sum-on-x86
+
+11. de casteljau implementation and comments
+https://forums.khronos.org/showthread.php/10264-Animations-in-1-4-1-release-notes-revision-A/page2?highlight=bezier
+https://forums.khronos.org/showthread.php/10644-Animation-Bezier-interpolation
+https://forums.khronos.org/showthread.php/10387-2D-Tangents-in-Bezier-Splines?p=34164&viewfull=1#post34164
+https://forums.khronos.org/showthread.php/10651-Animation-TCB-Spline-Interpolation-in-COLLADA?highlight=bezier

include/cglm/bezier.h

@@ -53,4 +53,77 @@ glm_bezier(float s, float p0, float c0, float c1, float p1) {
   return a + s * (c1 * xs3 + p1 * ss - a);
 }
 
+/*!
+ * @brief iterative way to solve cubic equation
+ *
+ * @param[in]  s    parameter between 0 and 1
+ * @param[in]  p0   begin point
+ * @param[in]  c0   control point 1
+ * @param[in]  c1   control point 2
+ * @param[in]  p1   end point
+ *
+ * @return parameter to use in cubic equation
+ */
+CGLM_INLINE
+float
+glm_decasteljau(float prm, float p0, float c0, float c1, float p1) {
+  float u, v, a, b, c, d, e, f;
+  int   i;
+
+  if (prm - p0 < CGLM_DECASTEL_SMALL)
+    return 0.0f;
+
+  if (p1 - prm < CGLM_DECASTEL_SMALL)
+    return 1.0f;
+
+  u  = 0.0f;
+  v  = 1.0f;
+
+  for (i = 0; i < CGLM_DECASTEL_MAX; i++) {
+    /* de Casteljau Subdivision */
+    a  = (p0 + c0) * 0.5f;
+    b  = (c0 + c1) * 0.5f;
+    c  = (c1 + p1) * 0.5f;
+    d  = (a  + b)  * 0.5f;
+    e  = (b  + c)  * 0.5f;
+    f  = (d  + e)  * 0.5f; /* this one is on the curve! */
+
+    /* The curve point is close enough to our wanted t */
+    if (fabsf(f - prm) < CGLM_DECASTEL_EPS)
+      return glm_clamp_zo((u  + v) * 0.5f);
+
+    /* dichotomy */
+    if (f < prm) {
+      p0 = f;
+      c0 = e;
+      c1 = c;
+      u  = (u  + v) * 0.5f;
+    } else {
+      c0 = a;
+      c1 = d;
+      p1 = f;
+      v  = (u  + v) * 0.5f;
+    }
+  }
+
+  return glm_clamp_zo((u  + v) * 0.5f);
+}
+
+/*!
+ * @brief solve cubic bezier equation
+ *
+ * @param[in]  s    parameter between 0 and 1
+ * @param[in]  p0   begin point
+ * @param[in]  c0   control point 1
+ * @param[in]  c1   control point 2
+ * @param[in]  p1   end point
+ *
+ * @return parameter to use in cubic equation
+ */
+CGLM_INLINE
+float
+glm_bezier_solve(float prm, float p0, float c0, float c1, float p1) {
+  return glm_decasteljau(prm, p0, c0, c1, p1);
+}
+
 #endif /* cglm_bezier_h */

include/cglm/call.h

@@ -28,6 +28,7 @@ extern "C" {
 #include "call/sphere.h"
 #include "call/ease.h"
 #include "call/curve.h"
+#include "call/bezier.h"
 
 #ifdef __cplusplus
 }

include/cglm/call/bezier.h

@@ -5,8 +5,8 @@
  * Full license can be found in the LICENSE file
  */
 
-#ifndef cglmc_curve_h
+#ifndef cglmc_bezier_h
-#define cglmc_curve_h
+#define cglmc_bezier_h
 #ifdef __cplusplus
 extern "C" {
 #endif
@@ -17,7 +17,15 @@ CGLM_EXPORT
 float
 glmc_bezier(float s, float p0, float c0, float c1, float p1);
 
+CGLM_EXPORT
+float
+glmc_decasteljau(float prm, float p0, float c0, float c1, float p1);
+
+CGLM_EXPORT
+float
+glmc_bezier_solve(float prm, float p0, float c0, float c1, float p1);
+
 #ifdef __cplusplus
 }
 #endif
-#endif /* cglmc_curve_h */
+#endif /* cglmc_bezier_h */

src/bezier.c

@@ -13,3 +13,15 @@ float
 glmc_bezier(float s, float p0, float c0, float c1, float p1) {
   return glm_bezier(s, p0, c0, c1, p1);
 }
+
+CGLM_EXPORT
+float
+glmc_decasteljau(float prm, float p0, float c0, float c1, float p1) {
+  return glm_decasteljau(prm, p0, c0, c1, p1);
+}
+
+CGLM_EXPORT
+float
+glmc_bezier_solve(float prm, float p0, float c0, float c1, float p1) {
+  return glm_bezier_solve(prm, p0, c0, c1, p1);
+}