Insert keyframes while preserving shape of curve

Apply the De Casteljau algorithm to split the Bèzier curve at the X
coordinate where the new key is inserted, and uses the result to update
both the newly inserted and surrounding handles.

For curves that use Auto keyframes this has been largely addressed by
the new algorithm from D2884. This commit extends this to non-auto
handles.

This code is heavily based on D3172 by Alexander Gavrilov (@angavrilov).

Manifest Task: https://developer.blender.org/T81353

source/blender/blenkernel/BKE_fcurve.h

@@ -267,6 +267,18 @@ typedef enum eFCU_Cycle_Type {
 
 eFCU_Cycle_Type BKE_fcurve_get_cycle_type(struct FCurve *fcu);
 
+/** Adjust Bezier handles of all three given BezTriples, so that `bezt` can be inserted between
+ * `prev` and `next` without changing the resulting curve shape.
+ *
+ * \param r_pdelta: return Y difference between `bezt` and the original curve value at its X
+ * position.
+ * \return Whether the split was succesful.
+ */
+bool BKE_bezt_subdivide_handles(struct BezTriple *bezt,
+                                struct BezTriple *prev,
+                                struct BezTriple *next,
+                                float *r_pdelta);
+
 /* -------- Curve Sanity --------  */
 
 void calchandles_fcurve(struct FCurve *fcu);

source/blender/blenkernel/intern/fcurve.c

@@ -1352,17 +1352,14 @@ void correct_bezpart(const float v1[2], float v2[2], float v3[2], const float v4
   }
 }
 
-/* find root ('zero') */
-static int findzero(float x, float q0, float q1, float q2, float q3, float *o)
+/** Find roots of cubic equation (c0 x³ + c1 x² + c2 x + c3)
+ * \return number of roots in `o`.
+ * NOTE: it is up to the caller to allocate enough memory for `o`. */
+static int solve_cubic(double c0, double c1, double c2, double c3, float *o)
 {
-  double c0, c1, c2, c3, a, b, c, p, q, d, t, phi;
+  double a, b, c, p, q, d, t, phi;
   int nr = 0;
 
-  c0 = q0 - x;
-  c1 = 3.0f * (q1 - q0);
-  c2 = 3.0f * (q0 - 2.0f * q1 + q2);
-  c3 = q3 - q0 + 3.0f * (q1 - q2);
-
   if (c3 != 0.0) {
     a = c2 / c3;
     b = c1 / c3;
@@ -1469,6 +1466,17 @@ static int findzero(float x, float q0, float q1, float q2, float q3, float *o)
   return 0;
 }
 
+/* Find root(s) ('zero') of a Bezier curve. */
+static int findzero(float x, float q0, float q1, float q2, float q3, float *o)
+{
+  const double c0 = q0 - x;
+  const double c1 = 3.0f * (q1 - q0);
+  const double c2 = 3.0f * (q0 - 2.0f * q1 + q2);
+  const double c3 = q3 - q0 + 3.0f * (q1 - q2);
+
+  return solve_cubic(c0, c1, c2, c3, o);
+}
+
 static void berekeny(float f1, float f2, float f3, float f4, float *o, int b)
 {
   float t, c0, c1, c2, c3;
@@ -1485,6 +1493,68 @@ static void berekeny(float f1, float f2, float f3, float f4, float *o, int b)
   }
 }
 
+/* Recompute handles to neatly subdivide the prev-next range at bezt. */
+bool BKE_bezt_subdivide_handles(struct BezTriple *bezt,
+                                struct BezTriple *prev,
+                                struct BezTriple *next,
+                                float *r_pdelta)
+{
+  /* The four points that make up this section of the Bezier curve. */
+  const float *prev_coords = prev->vec[1];
+  float *prev_handle_right = prev->vec[2];
+  float *next_handle_left = next->vec[0];
+  const float *next_coords = next->vec[1];
+
+  float *new_handle_left = bezt->vec[0];
+  const float *new_coords = bezt->vec[1];
+  float *new_handle_right = bezt->vec[2];
+
+  if (new_coords[0] <= prev_coords[0] || new_coords[0] >= next_coords[0]) {
+    /* The new keyframe is outside the (prev_coords, next_coords) range. */
+    return false;
+  }
+
+  /* Apply evaluation-time limits and compute the effective curve. */
+  correct_bezpart(prev_coords, prev_handle_right, next_handle_left, next_coords);
+  float roots[4];
+  if (!findzero(new_coords[0],
+                prev_coords[0],
+                prev_handle_right[0],
+                next_handle_left[0],
+                next_coords[0],
+                roots)) {
+    return false;
+  }
+
+  const float t = roots[0]; /* Percentage of the curve at which the split should occur. */
+  if (t <= 0.0f || t >= 1.0f) {
+    /* The split would occur outside the curve, which isn't possible. */
+    return false;
+  }
+
+  /* De Casteljau split, requires three iterations of splitting.
+   * See https://pomax.github.io/bezierinfo/#decasteljau */
+  float split1[3][2], split2[2][2], split3[2];
+  interp_v2_v2v2(split1[0], prev_coords, prev_handle_right, t);
+  interp_v2_v2v2(split1[1], prev_handle_right, next_handle_left, t);
+  interp_v2_v2v2(split1[2], next_handle_left, next_coords, t);
+  interp_v2_v2v2(split2[0], split1[0], split1[1], t);
+  interp_v2_v2v2(split2[1], split1[1], split1[2], t);
+  interp_v2_v2v2(split3, split2[0], split2[1], t);
+
+  /* Update the existing handles. */
+  copy_v2_v2(prev_handle_right, split1[0]);
+  copy_v2_v2(next_handle_left, split1[2]);
+
+  float diff_coords[2];
+  sub_v2_v2v2(diff_coords, new_coords, split3);
+  add_v2_v2v2(new_handle_left, split2[0], diff_coords);
+  add_v2_v2v2(new_handle_right, split2[1], diff_coords);
+
+  *r_pdelta = diff_coords[1];
+  return true;
+}
+
 /** \} */
 
 /* -------------------------------------------------------------------- */

source/blender/blenkernel/intern/fcurve_test.cc

@@ -210,4 +210,67 @@ TEST(evaluate_fcurve, ExtrapolationBezierKeys)
   BKE_fcurve_free(fcu);
 }
 
+TEST(fcurve_subdivide, BKE_bezt_subdivide_handles)
+{
+  FCurve *fcu = BKE_fcurve_create();
+
+  /* Insert two keyframes and set handles to something non-default. */
+  EXPECT_EQ(insert_vert_fcurve(fcu, 1.0f, 0.0f, BEZT_KEYTYPE_KEYFRAME, INSERTKEY_NO_USERPREF), 0);
+  EXPECT_EQ(insert_vert_fcurve(fcu, 13.0f, 2.0f, BEZT_KEYTYPE_KEYFRAME, INSERTKEY_NO_USERPREF), 1);
+
+  fcu->bezt[0].h1 = fcu->bezt[0].h2 = HD_FREE;
+  fcu->bezt[0].vec[0][0] = -5.0f;
+  fcu->bezt[0].vec[0][1] = 0.0f;
+  fcu->bezt[0].vec[2][0] = 2.0f;
+  fcu->bezt[0].vec[2][1] = 4.0f;
+
+  fcu->bezt[1].h1 = fcu->bezt[1].h2 = HD_FREE;
+  fcu->bezt[1].vec[0][0] = 13.0f;
+  fcu->bezt[1].vec[0][1] = -2.0f;
+  fcu->bezt[1].vec[2][0] = 16.0f;
+  fcu->bezt[1].vec[2][1] = -3.0f;
+
+  /* Create new keyframe point with defaults from insert_vert_fcurve(). */
+  BezTriple beztr;
+  const float x = 7.375f;  // at this X-coord, the FCurve should evaluate to 1.000f.
+  const float y = 1.000f;
+  beztr.vec[0][0] = x - 1.0f;
+  beztr.vec[0][1] = y;
+  beztr.vec[1][0] = x;
+  beztr.vec[1][1] = y;
+  beztr.vec[2][0] = x + 1.0f;
+  beztr.vec[2][1] = y;
+  beztr.h1 = beztr.h2 = HD_AUTO_ANIM;
+  beztr.ipo = BEZT_IPO_BEZ;
+
+  /* This should update the existing handles as well as the new BezTriple. */
+  float y_delta;
+  BKE_bezt_subdivide_handles(&beztr, &fcu->bezt[0], &fcu->bezt[1], &y_delta);
+
+  EXPECT_FLOAT_EQ(y_delta, 0.0f);
+
+  EXPECT_FLOAT_EQ(fcu->bezt[0].vec[0][0], -5.0f);  // Left handle should not be touched.
+  EXPECT_FLOAT_EQ(fcu->bezt[0].vec[0][1], 0.0f);
+  EXPECT_FLOAT_EQ(fcu->bezt[0].vec[1][0], 1.0f);  // Coordinates should not be touched.
+  EXPECT_FLOAT_EQ(fcu->bezt[0].vec[1][1], 0.0f);
+  EXPECT_FLOAT_EQ(fcu->bezt[0].vec[2][0], 1.5f);  // Right handle should be updated.
+  EXPECT_FLOAT_EQ(fcu->bezt[0].vec[2][1], 2.0f);
+
+  EXPECT_FLOAT_EQ(fcu->bezt[1].vec[0][0], 13.0f);  // Left handle should be updated.
+  EXPECT_FLOAT_EQ(fcu->bezt[1].vec[0][1], 0.0f);
+  EXPECT_FLOAT_EQ(fcu->bezt[1].vec[1][0], 13.0f);  // Coordinates should not be touched.
+  EXPECT_FLOAT_EQ(fcu->bezt[1].vec[1][1], 2.0f);
+  EXPECT_FLOAT_EQ(fcu->bezt[1].vec[2][0], 16.0f);  // Right handle should not be touched
+  EXPECT_FLOAT_EQ(fcu->bezt[1].vec[2][1], -3.0f);
+
+  EXPECT_FLOAT_EQ(beztr.vec[0][0], 4.5f);  // Left handle should be updated.
+  EXPECT_FLOAT_EQ(beztr.vec[0][1], 1.5f);
+  EXPECT_FLOAT_EQ(beztr.vec[1][0], 7.375f);  // Coordinates should not be touched.
+  EXPECT_FLOAT_EQ(beztr.vec[1][1], 1.0f);
+  EXPECT_FLOAT_EQ(beztr.vec[2][0], 10.250);  // Right handle should be updated.
+  EXPECT_FLOAT_EQ(beztr.vec[2][1], 0.5);
+
+  BKE_fcurve_free(fcu);
+}
+
 }  // namespace blender::bke::tests

source/blender/editors/animation/keyframing.c

@@ -483,6 +483,57 @@ int insert_bezt_fcurve(FCurve *fcu, const BezTriple *bezt, eInsertKeyFlags flag)
   return i;
 }
 
+/** Update the FCurve to allow insertion of `bezt` without modifying the curve shape.
+ *
+ * Checks whether it is necessary to apply Bezier subdivision due to involvement of non-auto
+ * handles. If necessary, changes `bezt` handles from Auto to Aligned.
+ *
+ * \param bezt: key being inserted
+ * \param prev: keyframe before that key
+ * \param next: keyframe after that key
+ */
+static void subdivide_nonauto_handles(const FCurve *fcu,
+                                      BezTriple *bezt,
+                                      BezTriple *prev,
+                                      BezTriple *next)
+{
+  if (prev->ipo != BEZT_IPO_BEZ || bezt->ipo != BEZT_IPO_BEZ) {
+    return;
+  }
+
+  /* Don't change Vector handles, or completely auto regions. */
+  const bool bezt_auto = BEZT_IS_AUTOH(bezt) || (bezt->h1 == HD_VECT && bezt->h2 == HD_VECT);
+  const bool prev_auto = BEZT_IS_AUTOH(prev) || (prev->h2 == HD_VECT);
+  const bool next_auto = BEZT_IS_AUTOH(next) || (next->h1 == HD_VECT);
+  if (bezt_auto && prev_auto && next_auto) {
+    return;
+  }
+
+  /* Subdivide the curve. */
+  float delta;
+  if (!BKE_bezt_subdivide_handles(bezt, prev, next, &delta)) {
+    return;
+  }
+
+  /* Decide when to force auto to manual. */
+  if (!BEZT_IS_AUTOH(bezt) || fabsf(delta) >= 0.001f) {
+    return;
+  }
+  if ((prev_auto || next_auto) && fcu->auto_smoothing == FCURVE_SMOOTH_CONT_ACCEL) {
+    const float hx = bezt->vec[1][0] - bezt->vec[0][0];
+    const float dx = bezt->vec[1][0] - prev->vec[1][0];
+
+    /* This mode always uses 1/3 of key distance for handle x size. */
+    const bool auto_works_well = fabsf(hx - dx / 3.0f) < 0.001f;
+    if (auto_works_well) {
+      return;
+    }
+  }
+
+  /* Turn off auto mode. */
+  bezt->h1 = bezt->h2 = HD_ALIGN;
+}
+
 /**
  * This function is a wrapper for #insert_bezt_fcurve(), and should be used when
  * adding a new keyframe to a curve, when the keyframe doesn't exist anywhere else yet.
@@ -562,14 +613,6 @@ int insert_vert_fcurve(
     return -1;
   }
 
-  /* don't recalculate handles if fast is set
-   * - this is a hack to make importers faster
-   * - we may calculate twice (due to autohandle needing to be calculated twice)
-   */
-  if ((flag & INSERTKEY_FAST) == 0) {
-    calchandles_fcurve(fcu);
-  }
-
   /* set handletype and interpolation */
   if ((fcu->totvert > 2) && (flag & INSERTKEY_REPLACE) == 0) {
     BezTriple *bezt = (fcu->bezt + a);
@@ -586,17 +629,21 @@ int insert_vert_fcurve(
       else if (a < fcu->totvert - 1) {
         bezt->ipo = (bezt + 1)->ipo;
       }
-    }
 
-    /* don't recalculate handles if fast is set
-     * - this is a hack to make importers faster
-     * - we may calculate twice (due to autohandle needing to be calculated twice)
-     */
-    if ((flag & INSERTKEY_FAST) == 0) {
-      calchandles_fcurve(fcu);
+      if (0 < a && a < (fcu->totvert - 1) && (flag & INSERTKEY_OVERWRITE_FULL) == 0) {
+        subdivide_nonauto_handles(fcu, bezt, bezt - 1, bezt + 1);
+      }
     }
   }
 
+  /* don't recalculate handles if fast is set
+   * - this is a hack to make importers faster
+   * - we may calculate twice (due to autohandle needing to be calculated twice)
+   */
+  if ((flag & INSERTKEY_FAST) == 0) {
+    calchandles_fcurve(fcu);
+  }
+
   /* return the index at which the keyframe was added */
   return a;
 }