recp · recp · Feb 3, 2019 · Jan 21, 2019 · Jan 22, 2019 · Jan 22, 2019
diff --git a/.gitignore b/.gitignore
@@ -69,3 +69,4 @@ win/cglm_test_*
 win/x64
 win/x85
 win/Debug
+cglm-test-ios*
diff --git a/CREDITS b/CREDITS
@@ -52,3 +52,12 @@ https://gamedev.stackexchange.com/questions/28395/rotating-vector3-by-a-quaterni
 
 9. Sphere AABB intersect
 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
diff --git a/README.md b/README.md
@@ -82,7 +82,11 @@ Currently *cglm* uses default clip space configuration (-1, 1) for camera functi
 - inline or pre-compiled function call
 - frustum (extract view frustum planes, corners...)
 - bounding box  (AABB in Frustum (culling), crop, merge...)
+- bounding sphere
 - project, unproject
+- easing functions
+- curves
+- curve interpolation helpers (S*M*C, deCasteljau...)
 - and other...
 
 <hr />

diff --git a/configure.ac b/configure.ac
@@ -29,6 +29,7 @@ LT_INIT
 # Checks for libraries.
 AC_CHECK_LIB([m], [floor])
 
+m4_ifdef([AM_SILENT_RULES], [AM_SILENT_RULES([yes])])
 AC_SYS_LARGEFILE
 
 # Checks for header files.

diff --git a/docs/source/api.rst b/docs/source/api.rst
@@ -46,3 +46,5 @@ Follow the :doc:`build` documentation for this
    io
    call
    sphere
+   curve
+   bezier
diff --git a/docs/source/bezier.rst b/docs/source/bezier.rst
@@ -0,0 +1,89 @@
+.. default-domain:: C
+
+Bezier
+================================================================================
+
+Header: cglm/bezier.h
+
+Common helpers for cubic bezier and similar curves.
+
+Table of contents (click to go):
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Functions:
+
+1. :c:func:`glm_bezier`
+2. :c:func:`glm_hermite`
+3. :c:func:`glm_decasteljau`
+
+Functions documentation
+~~~~~~~~~~~~~~~~~~~~~~~
+
+.. c:function:: float glm_bezier(float s, float p0, float c0, float c1, float p1)
+
+    | cubic bezier interpolation
+    | formula:
+
+    .. code-block:: text
+
+      B(s) = P0*(1-s)^3 + 3*C0*s*(1-s)^2 + 3*C1*s^2*(1-s) + P1*s^3
+
+    | similar result using matrix:
+
+    .. code-block:: text
+
+      B(s) = glm_smc(t, GLM_BEZIER_MAT, (vec4){p0, c0, c1, p1})
+
+    | glm_eq(glm_smc(...), glm_bezier(...)) should return TRUE
+
+    Parameters:
+      | *[in]*  **s**   parameter between 0 and 1
+      | *[in]*  **p0**  begin point
+      | *[in]*  **c0**  control point 1
+      | *[in]*  **c1**  control point 2
+      | *[in]*  **p1**  end point
+
+    Returns:
+        B(s)
+
+.. c:function:: float glm_hermite(float s, float p0, float t0, float t1, float p1)
+
+    | cubic hermite interpolation
+    | formula:
+
+    .. code-block:: text
+
+      H(s) = P0*(2*s^3 - 3*s^2 + 1) + T0*(s^3 - 2*s^2 + s) + P1*(-2*s^3 + 3*s^2) + T1*(s^3 - s^2)
+
+    | similar result using matrix:
+
+    .. code-block:: text
+
+      H(s) = glm_smc(t, GLM_HERMITE_MAT, (vec4){p0, p1, c0, c1})
+
+    | glm_eq(glm_smc(...), glm_hermite(...)) should return TRUE
+
+
+    Parameters:
+      | *[in]*  **s**   parameter between 0 and 1
+      | *[in]*  **p0**  begin point
+      | *[in]*  **t0**  tangent 1
+      | *[in]*  **t1**  tangent 2
+      | *[in]*  **p1**  end point
+
+    Returns:
+        B(s)
+
+.. c:function:: float glm_decasteljau(float prm, float p0, float c0, float c1, float p1)
+
+    | iterative way to solve cubic equation
+
+    Parameters:
+      | *[in]*  **prm** parameter between 0 and 1
+      | *[in]*  **p0**  begin point
+      | *[in]*  **c0**  control point 1
+      | *[in]*  **c1**  control point 2
+      | *[in]*  **p1**  end point
+
+    Returns:
+        parameter to use in cubic equation
diff --git a/docs/source/curve.rst b/docs/source/curve.rst
@@ -0,0 +1,41 @@
+.. default-domain:: C
+
+Curve
+================================================================================
+
+Header: cglm/curve.h
+
+Common helpers for common curves. For specific curve see its header/doc
+e.g bezier
+
+Table of contents (click to go):
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Functions:
+
+1. :c:func:`glm_smc`
+
+Functions documentation
+~~~~~~~~~~~~~~~~~~~~~~~
+
+.. c:function:: float  glm_smc(float s, mat4 m, vec4 c)
+
+    | helper function to calculate **S** * **M** * **C** multiplication for curves
+
+    | this function does not encourage you to use SMC, instead it is a helper if you use SMC.
+
+    | if you want to specify S as vector then use more generic glm_mat4_rmc() func.
+
+    | Example usage:
+
+    .. code-block:: c
+
+       Bs = glm_smc(s, GLM_BEZIER_MAT, (vec4){p0, c0, c1, p1})
+
+    Parameters:
+      | *[in]*  **s**  parameter between 0 and 1 (this will be [s3, s2, s, 1])
+      | *[in]*  **m**  basis matrix
+      | *[out]* **c**  position/control vector
+
+    Returns:
+        scalar value e.g. Bs
diff --git a/docs/source/mat4.rst b/docs/source/mat4.rst
@@ -45,6 +45,7 @@ Functions:
 #. :c:func:`glm_mat4_inv_fast`
 #. :c:func:`glm_mat4_swap_col`
 #. :c:func:`glm_mat4_swap_row`
+#. :c:func:`glm_mat4_rmc`
 
 Functions documentation
 ~~~~~~~~~~~~~~~~~~~~~~~
@@ -270,3 +271,20 @@ Functions documentation
       | *[in, out]*  **mat**   matrix
       | *[in]*       **row1**  row1
       | *[in]*       **row2**  row2
+
+.. c:function:: float  glm_mat4_rmc(vec4 r, mat4 m, vec4 c)
+
+    | **rmc** stands for **Row** * **Matrix** * **Column**
+
+    | helper for  R (row vector) * M (matrix) * C (column vector)
+
+    | the result is scalar because S * M = Matrix1x4 (row vector),
+    | then Matrix1x4 * Vec4 (column vector) = Matrix1x1 (Scalar)
+
+    Parameters:
+      | *[in]*  **r**  row vector or matrix1x4
+      | *[in]*  **m**  matrix4x4
+      | *[in]*  **c**  column vector or matrix4x1
+
+    Returns:
+        scalar value e.g. Matrix1x1
diff --git a/docs/source/vec4.rst b/docs/source/vec4.rst
@@ -58,11 +58,7 @@ Functions:
 #. :c:func:`glm_vec4_minv`
 #. :c:func:`glm_vec4_clamp`
 #. :c:func:`glm_vec4_lerp`
-#. :c:func:`glm_vec4_isnan`
-#. :c:func:`glm_vec4_isinf`
-#. :c:func:`glm_vec4_isvalid`
-#. :c:func:`glm_vec4_sign`
-#. :c:func:`glm_vec4_sqrt`
+#. :c:func:`glm_vec4_cubic`
 
 Functions documentation
 ~~~~~~~~~~~~~~~~~~~~~~~
@@ -401,3 +397,11 @@ Functions documentation
       | *[in]*  **to**     to value
       | *[in]*  **t**      interpolant (amount) clamped between 0 and 1
       | *[out]* **dest**   destination
+
+.. c:function:: void  glm_vec4_cubic(float s, vec4 dest)
+
+    helper to fill vec4 as [S^3, S^2, S, 1]
+
+    Parameters:
+      | *[in]*  **s**      parameter
+      | *[out]* **dest**   destination
diff --git a/include/cglm/bezier.h b/include/cglm/bezier.h
@@ -0,0 +1,152 @@
+/*
+ * Copyright (c), Recep Aslantas.
+ *
+ * MIT License (MIT), http://opensource.org/licenses/MIT
+ * Full license can be found in the LICENSE file
+ */
+
+#ifndef cglm_bezier_h
+#define cglm_bezier_h
+
+#define GLM_BEZIER_MAT_INIT  {{-1.0f,  3.0f, -3.0f,  1.0f},                   \
+                              { 3.0f, -6.0f,  3.0f,  0.0f},                   \
+                              {-3.0f,  3.0f,  0.0f,  0.0f},                   \
+                              { 1.0f,  0.0f,  0.0f,  0.0f}}
+#define GLM_HERMITE_MAT_INIT {{ 2.0f, -3.0f,  0.0f,  1.0f},                   \
+                              {-2.0f,  3.0f,  0.0f,  0.0f},                   \
+                              { 1.0f, -2.0f,  1.0f,  0.0f},                   \
+                              { 1.0f, -1.0f,  0.0f,  0.0f}}
+/* for C only */
+#define GLM_BEZIER_MAT  ((mat4)GLM_BEZIER_MAT_INIT)
+#define GLM_HERMITE_MAT ((mat4)GLM_HERMITE_MAT_INIT)
+
+#define CGLM_DECASTEL_EPS   1e-9
+#define CGLM_DECASTEL_MAX   1000
+#define CGLM_DECASTEL_SMALL 1e-20
+
+/*!
+ * @brief cubic bezier interpolation
+ *
+ * Formula:
+ *  B(s) = P0*(1-s)^3 + 3*C0*s*(1-s)^2 + 3*C1*s^2*(1-s) + P1*s^3
+ *
+ * similar result using matrix:
+ *  B(s) = glm_smc(t, GLM_BEZIER_MAT, (vec4){p0, c0, c1, p1})
+ *
+ * glm_eq(glm_smc(...), glm_bezier(...)) should return TRUE
+ *
+ * @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 B(s)
+ */
+CGLM_INLINE
+float
+glm_bezier(float s, float p0, float c0, float c1, float p1) {
+  float x, xx, ss, xs3, a;
+
+  x   = 1.0f - s;
+  xx  = x * x;
+  ss  = s * s;
+  xs3 = (s - ss) * 3.0f;
+  a   = p0 * xx + c0 * xs3;
+
+  return a + s * (c1 * xs3 + p1 * ss - a);
+}
+
+/*!
+ * @brief cubic hermite interpolation
+ *
+ * Formula:
+ *  H(s) = P0*(2*s^3 - 3*s^2 + 1) + T0*(s^3 - 2*s^2 + s)
+ *            + P1*(-2*s^3 + 3*s^2) + T1*(s^3 - s^2)
+ *
+ * similar result using matrix:
+ *  H(s) = glm_smc(t, GLM_HERMITE_MAT, (vec4){p0, p1, c0, c1})
+ *
+ * glm_eq(glm_smc(...), glm_hermite(...)) should return TRUE
+ *
+ * @param[in]  s    parameter between 0 and 1
+ * @param[in]  p0   begin point
+ * @param[in]  t0   tangent 1
+ * @param[in]  t1   tangent 2
+ * @param[in]  p1   end point
+ *
+ * @return H(s)
+ */
+CGLM_INLINE
+float
+glm_hermite(float s, float p0, float t0, float t1, float p1) {
+  float ss, d, a, b, c, e, f;
+
+  ss = s  * s;
+  a  = ss + ss;
+  c  = a  + ss;
+  b  = a  * s;
+  d  = s  * ss;
+  f  = d  - ss;
+  e  = b  - c;
+
+  return p0 * (e + 1.0f) + t0 * (f - ss + s) + t1 * f - p1 * e;
+}
+
+/*!
+ * @brief iterative way to solve cubic equation
+ *
+ * @param[in]  prm  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);
+}
+
+#endif /* cglm_bezier_h */
diff --git a/include/cglm/call.h b/include/cglm/call.h
@@ -27,6 +27,8 @@ extern "C" {
 #include "call/project.h"
 #include "call/sphere.h"
 #include "call/ease.h"
+#include "call/curve.h"
+#include "call/bezier.h"
 
 #ifdef __cplusplus
 }
-Original file line number
+Diff line change
@@ Expand Up / @@ -69,3 +69,4 @@ win/cglm_test_* @@
     win/x64
     win/x85
     win/Debug
+    cglm-test-ios*