curve: cubic hermite intrpolation

include/cglm/bezier.h

@@ -8,13 +8,17 @@
 #ifndef cglm_bezier_h
 #define cglm_bezier_h
 
-#define GLM_BEZIER_MAT_INIT  {{-1.0f,  3.0f, -3.0f, 1.0f},                    \
+#define GLM_BEZIER_MAT_INIT  {{-1.0f,  3.0f, -3.0f,  1.0f},                   \
-                              { 3.0f, -6.0f,  3.0f, 0.0f},                    \
+                              { 3.0f, -6.0f,  3.0f,  0.0f},                   \
-                              {-3.0f,  3.0f,  0.0f, 0.0f},                    \
+                              {-3.0f,  3.0f,  0.0f,  0.0f},                   \
-                              { 1.0f,  0.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_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
@@ -54,10 +58,46 @@ glm_bezier(float s, float p0, float c0, float c1, float p1) {
 }
 
 /*!
- * @brief iterative way to solve cubic equation
+ * @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 B(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
@@ -112,7 +152,7 @@ glm_decasteljau(float prm, float p0, float c0, float c1, float p1) {
 /*!
  * @brief solve cubic bezier equation
  *
- * @param[in]  s    parameter between 0 and 1
+ * @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

include/cglm/call/bezier.h

@@ -17,6 +17,10 @@ CGLM_EXPORT
 float
 glmc_bezier(float s, float p0, float c0, float c1, float p1);
 
+CGLM_EXPORT
+float
+glmc_hermite(float s, float p0, float t0, float t1, float p1);
+
 CGLM_EXPORT
 float
 glmc_decasteljau(float prm, float p0, float c0, float c1, float p1);

makefile.am

@@ -58,7 +58,8 @@ cglm_HEADERS = include/cglm/version.h \
                include/cglm/project.h \
                include/cglm/sphere.h \
                include/cglm/ease.h \
-               include/cglm/curve.h
+               include/cglm/curve.h \
+               include/cglm/bezier.h
 
 cglm_calldir=$(includedir)/cglm/call
 cglm_call_HEADERS = include/cglm/call/mat4.h \
@@ -76,7 +77,8 @@ cglm_call_HEADERS = include/cglm/call/mat4.h \
                     include/cglm/call/project.h \
                     include/cglm/call/sphere.h \
                     include/cglm/call/ease.h \
-                    include/cglm/call/curve.h
+                    include/cglm/call/curve.h \
+                    include/cglm/call/bezier.h
 
 cglm_simddir=$(includedir)/cglm/simd
 cglm_simd_HEADERS = include/cglm/simd/intrin.h \
@@ -112,7 +114,8 @@ libcglm_la_SOURCES=\
     src/project.c \
     src/sphere.c \
     src/ease.c \
-    src/curve.c
+    src/curve.c \
+    src/bezier.c
 
 test_tests_SOURCES=\
     test/src/test_common.c \
@@ -126,7 +129,8 @@ test_tests_SOURCES=\
     test/src/test_vec4.c \
     test/src/test_vec3.c \
     test/src/test_mat3.c \
-    test/src/test_affine.c
+    test/src/test_affine.c \
+    test/src/test_bezier.c
 
 all-local:
 	sh ./post-build.sh

src/bezier.c

@@ -14,6 +14,12 @@ glmc_bezier(float s, float p0, float c0, float c1, float p1) {
   return glm_bezier(s, p0, c0, c1, p1);
 }
 
+CGLM_EXPORT
+float
+glmc_hermite(float s, float p0, float t0, float t1, float p1) {
+  return glm_hermite(s, p0, t0, t1, p1);
+}
+
 CGLM_EXPORT
 float
 glmc_decasteljau(float prm, float p0, float c0, float c1, float p1) {

test/src/test_bezier.c

@@ -21,6 +21,20 @@ test_bezier_plain(float s, float p0, float c0, float c1, float p1) {
   return p0 * xxx + 3.0f * (c0 * s * xx + c1 * ss * x) + p1 * sss;
 }
 
+CGLM_INLINE
+float
+test_hermite_plain(float s, float p0, float t0, float t1, float p1) {
+  float ss, sss;
+
+  ss  = s  * s;
+  sss = ss * s;
+
+  return p0 * (2.0f * sss - 3.0f * ss + 1.0f)
+       + t0 * (sss - 2.0f * ss + s)
+       + p1 * (-2.0f * sss + 3.0f * ss)
+       + t1 * (sss - ss);
+}
+
 void
 test_bezier(void **state) {
   float s, p0, p1, c0, c1, smc, Bs, Bs_plain;
@@ -31,6 +45,7 @@ test_bezier(void **state) {
   c0       = test_rand();
   c1       = test_rand();
 
+  /* test cubic bezier */
   smc      = glm_smc(s, GLM_BEZIER_MAT, (vec4){p0, c0, c1, p1});
   Bs       = glm_bezier(s, p0, c0, c1, p1);
   Bs_plain = test_bezier_plain(s, p0, c0, c1, p1);
@@ -38,4 +53,13 @@ test_bezier(void **state) {
   assert_true(glm_eq(Bs,  Bs_plain));
   assert_true(glm_eq(smc, Bs_plain));
   assert_true(glm_eq(Bs,  smc));
+
+  /* test cubic hermite */
+  smc      = glm_smc(s, GLM_HERMITE_MAT, (vec4){p0, p1, c0, c1});
+  Bs       = glm_hermite(s, p0, c0, c1, p1);
+  Bs_plain = test_hermite_plain(s, p0, c0, c1, p1);
+
+  assert_true(glm_eq(Bs,  Bs_plain));
+  assert_true(glm_eq(smc, Bs_plain));
+  assert_true(glm_eq(Bs,  smc));
 }

win/cglm.vcxproj

@@ -20,6 +20,7 @@
   </ItemGroup>
   <ItemGroup>
     <ClCompile Include="..\src\affine.c" />
+    <ClCompile Include="..\src\bezier.c" />
     <ClCompile Include="..\src\box.c" />
     <ClCompile Include="..\src\cam.c" />
     <ClCompile Include="..\src\curve.c" />
@@ -40,9 +41,11 @@
   <ItemGroup>
     <ClInclude Include="..\include\cglm\affine-mat.h" />
     <ClInclude Include="..\include\cglm\affine.h" />
+    <ClInclude Include="..\include\cglm\bezier.h" />
     <ClInclude Include="..\include\cglm\box.h" />
     <ClInclude Include="..\include\cglm\call.h" />
     <ClInclude Include="..\include\cglm\call\affine.h" />
+    <ClInclude Include="..\include\cglm\call\bezier.h" />
     <ClInclude Include="..\include\cglm\call\box.h" />
     <ClInclude Include="..\include\cglm\call\cam.h" />
     <ClInclude Include="..\include\cglm\call\curve.h" />

win/cglm.vcxproj.filters

@@ -87,6 +87,9 @@
     <ClCompile Include="..\src\curve.c">
       <Filter>src</Filter>
     </ClCompile>
+    <ClCompile Include="..\src\bezier.c">
+      <Filter>src</Filter>
+    </ClCompile>
   </ItemGroup>
   <ItemGroup>
     <ClInclude Include="..\src\config.h">
@@ -248,5 +251,11 @@
     <ClInclude Include="..\include\cglm\curve.h">
       <Filter>include\cglm</Filter>
     </ClInclude>
+    <ClInclude Include="..\include\cglm\bezier.h">
+      <Filter>include\cglm</Filter>
+    </ClInclude>
+    <ClInclude Include="..\include\cglm\call\bezier.h">
+      <Filter>include\cglm\call</Filter>
+    </ClInclude>
   </ItemGroup>
 </Project>