[Bf-blender-cvs] [8861c5fd35f] smooth-curves: Reimplement curve smoothing in a clear way without premature optimizations.
Alexander Gavrilov
noreply at git.blender.org
Sat Apr 29 20:59:38 CEST 2017
Commit: 8861c5fd35f5eae97a36d3d188299cca1eeab02c
Author: Alexander Gavrilov
Date: Sat Apr 29 21:53:07 2017 +0300
Branches: smooth-curves
https://developer.blender.org/rB8861c5fd35f5eae97a36d3d188299cca1eeab02c
Reimplement curve smoothing in a clear way without premature optimizations.
The original code tries too hard to avoid allocating more memory, which
reduces clarity and is not necessary here as there aren't that many
f-curve keyframes in any case. Rewriting the code with a generic solver
somehow fixed the behavior of the curve in the totally unclamped case,
which suggests that there was a bug hidden in the optimizations.
In addition, the new solver allows treating the case when an Auto curve
segment ends at a Vector handle specially as a zero acceleration boundary
condition.
===================================================================
M source/blender/blenkernel/BKE_curve.h
M source/blender/blenkernel/intern/curve.c
M source/blender/blenkernel/intern/fcurve.c
===================================================================
diff --git a/source/blender/blenkernel/BKE_curve.h b/source/blender/blenkernel/BKE_curve.h
index e11e9830160..8f5b91991f4 100644
--- a/source/blender/blenkernel/BKE_curve.h
+++ b/source/blender/blenkernel/BKE_curve.h
@@ -201,7 +201,7 @@ void BKE_nurb_bpoint_calc_normal(struct Nurb *nu, struct BPoint *bp, float r_nor
void BKE_nurb_handle_calc(struct BezTriple *bezt, struct BezTriple *prev, struct BezTriple *next,
const bool is_fcurve);
-bool BKE_nurb_handle_calc_smooth(struct BezTriple *first, int count);
+void BKE_nurb_handle_calc_smooth(struct BezTriple *first, int count);
void BKE_nurb_handle_calc_simple(struct Nurb *nu, struct BezTriple *bezt);
void BKE_nurb_handle_calc_simple_auto(struct Nurb *nu, struct BezTriple *bezt);
diff --git a/source/blender/blenkernel/intern/curve.c b/source/blender/blenkernel/intern/curve.c
index 5ef2350c1e5..e3810160800 100644
--- a/source/blender/blenkernel/intern/curve.c
+++ b/source/blender/blenkernel/intern/curve.c
@@ -3398,16 +3398,57 @@ static void calchandlesNurb_intern(Nurb *nu, bool skip_align)
#define MAX_AUTO_C_VALUE 1000.f
/*
+ * Solves a tridiagonal system of equations:
+ *
+ * a[i] * x[i-1] + b[i] * x[i] + c[i] * x[i+1] = d[i]
+ *
+ * Ignores a[0] and c[count-1]. Uses Thomas algorithm, e.g. see wiki.
+ */
+static void tridiagonal_solve(const float *a, const float *b, const float *c, const float *d, float *x, const int count)
+{
+ float *c1 = (float *)MEM_mallocN(sizeof(float)*count, "tridiagonal_c1");
+ float *d1 = (float *)MEM_mallocN(sizeof(float)*count, "tridiagonal_d1");
+
+ int i;
+ double c_prev, d_prev, x_prev;
+
+ /* forward pass */
+
+ c1[0] = c_prev = ((double)c[0]) / b[0];
+ d1[0] = d_prev = ((double)d[0]) / b[0];
+
+ for (i = 1; i < count; i++)
+ {
+ double denum = b[i] - a[i]*c_prev;
+
+ c1[i] = c_prev = c[i] / denum;
+ d1[i] = d_prev = (d[i] - a[i]*d_prev) / denum;
+ }
+
+ /* back pass */
+
+ x[--i] = x_prev = d_prev;
+
+ while (--i >= 0)
+ {
+ x[i] = x_prev = d1[i] - c1[i] * x_prev;
+ }
+
+ MEM_freeN(c1);
+ MEM_freeN(d1);
+}
+
+/*
* This function computes the handles of a series of auto bezier points
* on the basis of 'no acceleration discontinuities' at the points.
* The first and last bezier points are considered 'fixed' (their handles are not touched)
* The result is the smoothest possible trajectory going through intemediate points.
* The difficulty is that the handles depends on their neighbours.
*
- * The exact solution is a recursive formula that computes the first point handles
- * as the sum of the influence of the nth following points. The influence decreases
- * rapidly (as 1/4^n) so that just a few iterations are necessary. Once the first point
- * is computed, it becomes 'fixed' and the next point is computed.
+ * The exact solution is found by solving a tridiagonal matrix equation formed
+ * by the continuity and boundary conditions. Although theoretically handle position
+ * is affected by all other points of the curve segment, in practice the influence
+ * decreases exponentially with distance.
*
* Note: this algorithm assumes that the handle horizontal size if always 1/3 of the
* of the interval to the next point. This rule ensures linear interpolation of time.
@@ -3424,160 +3465,179 @@ static void calchandlesNurb_intern(Nurb *nu, bool skip_align)
* |-------t1---------t2--------- ~ --------tN-------------------> time (co 0)
*
*
- * Definitions:
- *
- * for n>=2:
- * l(n) = t(n)/t(n-1)
- * b(n) = 2*(1+l(n))
- * a(n) = (y(n)-y(n-1))+l(n)*l(n)*(y(n-1)-y(n-2))
- * p(n) = l(n)*p(n-1) for n>=3
- * = 1 for n=2
- * c(n) = l(n)*(c(n-1)*b(n)-c(n-2)*p(n) for n>=3
- * = b(2) for n=2
- * = 1 for n=1
+ * Mathematical basis:
*
- * Then the nth approximation of right handle height of the first point is:
+ * 1. Handle lengths on either side of each point are connected by a factor
+ * ensuring continuity of the first derivative:
*
- * for n>=2:
- * h(n) = (b(n)*l(n)*c(n-1)*h(n-1)-c(n-2)*h(n-2)*l(n-1)*l(n)*l(n)+((-1)^n)*a(n))/c(n)
- * h(0) = 0
- * h(1) = 0
+ * l[i] = t[i+1]/t[i]
*
- * 1/c(n) is the fraction of influence of the nth points to the first one.
- * You can stop the iterations when c(n) > 1000
- *
- * Example: if there are only 3 points in the series (N=2), the middle point
- * is the auto handle and the optimal right handle height = h(2) = a(2)/c(2) = a(2)/b(2)
+ * 2. The tridiagonal system is formed by the following equation, which is derived
+ * by differentiating the bezier curve and specifies second derivative continuity
+ * at every point:
*
- *
+ * l[i]^2 * h[i-1] + (2*l[i]+2) * h[i] + 1/l[i+1] * h[i+1] = (y[i]-y[i-1])*l[i]^2 + y[i+1]-y[i]
*/
-bool BKE_nurb_handle_calc_smooth(BezTriple *first, int count)
+
+static void nurb_lock_unknown(float *a, float *b, float *c, float *d, int i, float value)
{
- static int mod3[3] = { 1, 2, 0 };
- float *pl;
- float *l;
- float *a;
- float *b;
- float *y;
- float c[3]; /* only need to keep 3 C coef: current and previous 2 */
- float h[3]; /* only need to keey 3 h value: current and previous 2 */
- float t[3];
- float sign, dy, f;
- BezTriple *bezt, *prev, *next;
- int i, j, n, n1, n2, maxpoints;
- bool final = true;
-
- if (count < 3)
- return final;
-
- pl = (float *)MEM_mallocN(sizeof(float)*count, "ha_prod_lamba");
- l = (float *)MEM_mallocN(sizeof(float)*count, "ha_lamba");
- a = (float *)MEM_mallocN(sizeof(float)*count, "ha_a");
- b = (float *)MEM_mallocN(sizeof(float)*count, "ha_b");
- y = (float *)MEM_mallocN(sizeof(float)*count, "ha_y");
-
- bezt = first;
- pl[2] = 1.0f;
- n=0;
- n1=2; /* n-1 mod 3 */
- n2=1; /* n-2 mod 3 */
- for (i=0; i<count; i++) {
- t[n] = bezt->vec[1][0];
- if (i==0) {
- y[i] = bezt->vec[2][1];
- } else if (i==count-1) {
- y[i] = bezt->vec[0][1];
- } else {
- y[i] = bezt->vec[1][1];
- }
- if (i >= 2) {
- l[i] = (t[n]-t[n1])/(t[n1]-t[n2]);
- b[i] = 2.0f*(1.0f+l[i]);
- a[i] = (y[i]-y[i-1]) + l[i]*l[i]*(y[i-1]-y[i-2]);
- if (i>=3)
- pl[i] = pl[i-1]*l[i];
- }
- n = mod3[n];
- n1 = mod3[n1];
- n2 = mod3[n2];
- bezt++;
+ a[i] = c[i] = 0.0f;
+ b[i] = 1.0f;
+ d[i] = value;
+}
+
+static void nurb_clamp_right(float *hmax, float *hmin, float *y, int i)
+{
+ hmax[i] = min_ff(hmax[i], max_ff(y[i+1]-y[i], 0.0f));
+ hmin[i] = max_ff(hmin[i], min_ff(y[i+1]-y[i], 0.0f));
+}
+
+static void nurb_clamp_left(float *hmax, float *hmin, float *y, float *l, int i)
+{
+ hmax[i] = min_ff(hmax[i], max_ff(y[i]-y[i-1], 0.0f)*l[i]);
+ hmin[i] = max_ff(hmin[i], min_ff(y[i]-y[i-1], 0.0f)*l[i]);
+}
+
+void BKE_nurb_handle_calc_smooth(BezTriple *bezt, int count)
+{
+ float *x, *y, *l, *a, *b, *c, *d, *h, *hmax, *hmin;
+ float **arrays[] = { &x, &y, &l, &a, &b, &c, &d, &h, &hmax, &hmin };
+
+ if (count < 2)
+ return;
+
+ if (count == 2)
+ {
+ if (!ELEM(HD_VECT, bezt[0].h2, bezt[1].h1) || (bezt[0].h2 == bezt[1].h1))
+ return;
}
- prev = first;
- bezt = first+1;
- next = first+2;
- dy = 0.f;
- for (i=1; i<count-1; i++) {
- /* compute ideal slope as if the next bezier point was fixed */
- /* adjust a[n] based on the previous computed handle (assumed fixed) */
- a[i+1] -= l[i+1]*l[i+1]*dy;
- c[0] = 0.f;
- c[1] = 1.f;
- c[2] = b[i+1];
- h[0] = 0.f;
- h[1] = 0.f;
- h[2] = a[i+1]/c[2];
- /* evaluate the influence of the remaing bezier points */
- n=0;
- n1=2;
- n2=1;
- sign=-1;
- maxpoints = (i+MAX_AUTO_HANDLE_INFLUENCE+1);
- if (maxpoints > count)
- maxpoints = count;
- for (j=i+2; j<maxpoints; j++) {
- c[n] = l[j]*(c[n1]*b[j]-c[n2]*pl[j]);
- h[n] = (b[j]*l[j]*c[n1]*h[n1]-c[n2]*h[n2]*l[j-1]*l[j]*l[j]+sign*a[j])/c[n];
- sign = -sign;
- n = mod3[n];
- n1 = mod3[n1];
- n2 = mod3[n2];
- if (c[n1] > MAX_AUTO_C_VALUE)
- break;
- }
- /* h[n-1] is the value of the right handle adjusted for 0 accelaration step */
- dy = h[n1];
- if (bezt->h1 == HD_AUTO_ANIM || bezt->h2 == HD_AUTO_ANIM) {
- /* avoid overshoot */
- bool hard_clamped = false;
- if (dy > 0.f) {
- if (dy > (f = (next->vec[1][1]-y[i]))) {
- dy = (f < 0.f) ? 0.f : f;
- hard_clamped = hard_clamped || f >= 0.f;
- }
- if (dy > (f=l[i+1]*(y[i]-prev->vec[1][1]))) {
- dy = (f < 0.f) ? 0.f : f;
- hard_clamped = hard_clamped || f >= 0.f;
- }
- } else {
- if (dy < (f = (next->vec[1][1]-y[i]))) {
- dy = (f > 0.f) ? 0.f : f;
- hard_clamped = hard_clamped || f <= 0.f;
- }
- if (dy < (f = l[i+1]*(y[i]-prev->vec[1][1]))) {
- dy = (f > 0.f) ? 0.f : f;
- hard_clamped = hard_clamped || f <= 0.f;
+
+ /* allocate all */
+ for (int i = 0; i < sizeof(arrays)/sizeof(float*); i++)
+ *arrays[i] = (float *)MEM_mallocN(sizeof(float)*count, "nurb_calc_smooth_tmp");
+
+ /* point locations and ratio of x intervals */
+
+ for (int i = 0; i < count; i++)
+ {
+ x[i] = bezt[i].vec[1][0];
+ y[i] = bezt[i].vec[1][1];
+ hmax[i] = FLT_MAX;
+ hmin[i] = -FLT_MAX;
+ }
+
+ l[0] = l[count-1] = 1.0f;
+
+ for (int i = 1; i < count-1; i++)
+ {
+ l[i] = (x[i+1] - x[i]) / (x[i] - x[i-1]);
+ }
+
+ /* compute handle clamp ranges */
+
+ bool clamped = false;
+
+ for (int i = 0; i < count-1; i++)
+ {
+ if (clamped)
+ nurb_clamp_left(hmax, hmin, y, l, i);
+
+ clamped = ELEM(HD_AUTO_ANIM, bezt[i].h2, bezt[i+1].h1);
+
+ if (clamped)
+ nurb_clamp_right(hmax, hmin, y, i);
+ }
+
+ if (clamped)
+ nurb_clamp_left(hmax, hmin, y, l, count-1);
+
+ /* main tridiagonal system of equations */
+
+ for (int i = 1; i < count-1; i++)
+ {
+ a[i] = l[i]*l[i];
+ b[i] = 2.0f*(l[i] + 1);
+ c[i] = 1.0f/l[i+1];
+ d[i] = (y[i]-y[i-1])*l[i]*l[i] + (y[i+1]-y[i]);
+ }
+
+ /* boundary condition: fixed handles or zero curvature */
+
+ if (bezt[0].h2 == HD_VECT)
+ {
+ a[0] = 0.0f;
+ b[0] = 2.0f;
+ c[0] = 1.0f/l[1];
+ d[0] = y[1]-y[0];
+ }
+ else
+ nurb_lock_unknown(a, b, c, d, 0, bezt[0].vec[2][1] - y[0]);
+
+ if (bezt[count-1].h1 == HD_VECT)
+ {
+ int i = count-1;
+ a[i] = l[i]*l[i];
+ b[i] = 2.0f*l[i];
+ c[i] = 0.0f;
+ d[i] = (y[i]-y[i-1])*l[i]*l[i];
+ }
+ else
+ nurb_lock_unknown(a, b, c, d, count-1, y[count-1] - bezt[count-1].vec[0][1]);
+
+ /* solve and clamp until done */
+
+ bool overshoot;
+
+ do
+ {
+ tridiagonal_solve(a, b,
@@ Diff output truncated at 10240 characters. @@
More information about the Bf-blender-cvs
mailing list