| Message ID | 20221003014404.749-1-tikuma@hotmail.com |
|---|---|
| State | New |
| Headers | show |
| Series | [FFmpeg-devel] avfilter/vf_curves: add PCHIP interpolator and interp option | expand |
| Context | Check | Description |
|---|---|---|
| andriy/make_x86 | success | Make finished |
| andriy/make_fate_x86 | success | Make fate finished |
On 10/3/22, Takeshi (Kesh) Ikuma <tikuma@gmail.com> wrote: > summary: This patch modifies the `curves` filter with new `interp` option > to let user pick the existing natural cubic spline interpolation > and the new PCHIP interapolation. > > reason: The natural cubic spline does not impose monotonicity between > the keypoints. As such, the fitted curve may vary wildly against > user's intension. The PCHIP interpolation is not as smooth as > the natural spline but guarantees the monotonicity. Providing > both options enhances users experience (e.g., reduces the number > of keypoints to realize the desired curve). See the related bug > report for the example of an ill-interpolated curve. > > alternate solution: > Both Photoshop and GIMP appear to use monotonic interpolation in > their curve tools, which were the models for this filter. As > such, an alternate solution is to drop the natural spline and > go without the `interp` option. > > related bug report: https://trac.ffmpeg.org/ticket/9947 (filed by myself) > > Signed-off-by: Takeshi (Kesh) Ikuma <tikuma@hotmail.com> > > modified: doc/filters.texi > modified: libavfilter/vf_curves.c > --- > doc/filters.texi | 23 +++- > libavfilter/vf_curves.c | 255 ++++++++++++++++++++++++++++++++++++---- > 2 files changed, 253 insertions(+), 25 deletions(-) > Looks like this new interpolator implementationi hangs with 16bit pixel formats and with some presets: ffplay input.video -vf format=yuv422p16,curves=vintage:interp=1,format=yuv422p16 > diff --git a/doc/filters.texi b/doc/filters.texi > index d0f718678c..08a79644e1 100644 > --- a/doc/filters.texi > +++ b/doc/filters.texi > @@ -10389,11 +10389,15 @@ By default, a component curve is defined by the > two points @var{(0;0)} and > "adjusted" to its own value, which means no change to the image. > > The filter allows you to redefine these two points and add some more. A > new > -curve (using a natural cubic spline interpolation) will be define to pass > -smoothly through all these new coordinates. The new defined points needs to > be > -strictly increasing over the x-axis, and their @var{x} and @var{y} values > must > -be in the @var{[0;1]} interval. If the computed curves happened to go > outside > -the vector spaces, the values will be clipped accordingly. > +curve will be define to pass smoothly through all these new coordinates. > The > +new defined points needs to be strictly increasing over the x-axis, and > their > +@var{x} and @var{y} values must be in the @var{[0;1]} interval. The curve > is > +formed by using a natural or monotonic cubic spline interpolation, > depending > +on the @var{interp} option (default: @code{natural}). The @code{natural} > +spline produces a smoother curve in general while the monotonic > (@code{pchip}) > +spline guarantees the transitions between the specified points to be > +monotonic. If the computed curves happened to go outside the vector > spaces, > +the values will be clipped accordingly. > > The filter accepts the following options: > > @@ -10437,6 +10441,15 @@ options. In this case, the unset component(s) will > fallback on this > Specify a Photoshop curves file (@code{.acv}) to import the settings from. > @item plot > Save Gnuplot script of the curves in specified file. > +@item interp > +Specify the kind of interpolation. Available algorithms are: > +@table @samp > +@item natural > +Natural cubic spline using a piece-wise cubic polynomial that is twice > continuously differentiable. > +@item pchip > +Monotonic cubic spline using a piecewise cubic Hermite interpolating > polynomial (PCHIP). > +@end table > + > @end table > > To avoid some filtergraph syntax conflicts, each key points list need to > be > diff --git a/libavfilter/vf_curves.c b/libavfilter/vf_curves.c > index 498b06f6e5..d0efa380e1 100644 > --- a/libavfilter/vf_curves.c > +++ b/libavfilter/vf_curves.c > @@ -58,6 +58,12 @@ enum preset { > NB_PRESETS, > }; > > +enum interp { > + INTERP_NATURAL, > + INTERP_PCHIP, > + NB_INTERPS, > +}; > + > typedef struct CurvesContext { > const AVClass *class; > int preset; > @@ -73,6 +79,7 @@ typedef struct CurvesContext { > int is_16bit; > int depth; > int parsed_psfile; > + int interp; > > int (*filter_slice)(AVFilterContext *ctx, void *arg, int jobnr, int > nb_jobs); > } CurvesContext; > @@ -107,6 +114,9 @@ static const AVOption curves_options[] = { > { "all", "set points coordinates for all components", > OFFSET(comp_points_str_all), AV_OPT_TYPE_STRING, {.str=NULL}, .flags = FLAGS > }, > { "psfile", "set Photoshop curves file name", OFFSET(psfile), > AV_OPT_TYPE_STRING, {.str=NULL}, .flags = FLAGS }, > { "plot", "save Gnuplot script of the curves in specified file", > OFFSET(plot_filename), AV_OPT_TYPE_STRING, {.str=NULL}, .flags = FLAGS }, > + { "interp", "specify the kind of interpolation", OFFSET(interp), > AV_OPT_TYPE_INT, {.i64=INTERP_NATURAL}, INTERP_NATURAL, NB_INTERPS-1, FLAGS, > "interp_name" }, > + { "natural", "natural cubic spline", 0, AV_OPT_TYPE_CONST, > {.i64=INTERP_NATURAL}, 0, 0, FLAGS, "interp_name" }, > + { "pchip", "monotonically cubic interpolation", 0, > AV_OPT_TYPE_CONST, {.i64=INTERP_PCHIP}, 0, 0, FLAGS, "interp_name" }, > { NULL } > }; > > @@ -336,21 +346,230 @@ end: > av_free(h); > av_free(r); > return ret; > + > +} > + > +#define SIGN(x) (x>0.0?1:x<0.0?-1:0) > + > +/** > + * Evalaute the derivative of an edge endpoint > + * > + * @param h0 input interval of the interval closest to the edge > + * @param h1 input interval of the interval next to the closest > + * @param m0 linear slope of the interval closest to the edge > + * @param m1 linear slope of the intervalnext to the closest > + * @return edge endpoint derivative > + * > + * Based on scipy.interpolate._edge_case() > + * > https://github.com/scipy/scipy/blob/2e5883ef7af4f5ed4a5b80a1759a45e43163bf3f/scipy/interpolate/_cubic.py#L239 > + * which is a python implementation of the special case endpoints, as > suggested in > + * Cleve Moler, Numerical Computing with MATLAB, Chap 3.6 (pchiptx.m) > +*/ > +static double pchip_edge_case(double h0, double h1, double m0, double m1) > +{ > + int mask, mask2; > + double d; > + > + d = ((2 * h0 + h1) * m0 - h0 * m1) / (h0 + h1); > + > + mask = SIGN(d) != SIGN(m0); > + mask2 = (SIGN(m0) != SIGN(m1)) && (fabs(d) > 3. * fabs(m0)); > + > + if (mask) d = 0.0; > + else if (mask2) d = 3.0 * m0; > + > + return d; > } > > -#define DECLARE_INTERPOLATE_FUNC(nbits) > \ > -static int interpolate##nbits(void *log_ctx, uint16_t *y, > \ > - const struct keypoint *points) > \ > -{ > \ > - return interpolate(log_ctx, y, points, nbits); > \ > +/** > + * Evalaute the piecewise polynomial derivatives at endpoints > + * > + * @param n input interval of the interval closest to the edge > + * @param hk input intervals > + * @param mk linear slopes over intervals > + * @param dk endpoint derivatives (output) > + * @return 0 success > + * > + * Based on scipy.interpolate._find_derivatives() > + * > https://github.com/scipy/scipy/blob/2e5883ef7af4f5ed4a5b80a1759a45e43163bf3f/scipy/interpolate/_cubic.py#L254 > +*/ > + > +static int pchip_find_derivatives(const int n, const double *hk, const > double *mk, double *dk) > +{ > + int ret = 0; > + const int m = n - 1; > + int8_t *smk; > + > + smk = av_malloc(n); > + if (!smk) { > + ret = AVERROR(ENOMEM); > + goto end; > + } > + > + /* smk = sgn(mk) */ > + for (int i = 0; i < n; ++i) smk[i] = SIGN(mk[i]); > + > + /* check the strict monotonicity */ > + for (int i = 0; i < m; ++i) { > + int8_t condition = (smk[i + 1] != smk[i]) || (mk[i + 1] == 0) || > (mk[i] == 0); > + if (condition) { > + dk[i + 1] = 0.0; > + } else { > + double w1 = 2 * hk[i + 1] + hk[i]; > + double w2 = hk[i + 1] + 2 * hk[i]; > + dk[i + 1] = (w1 + w2) / (w1 / mk[i] + w2 / mk[i + 1]); > + } > + } > + > + dk[0] = pchip_edge_case(hk[0], hk[1], mk[0], mk[1]); > + dk[n] = pchip_edge_case(hk[n - 1], hk[n - 2], mk[n - 1], mk[n - 2]); > + > +end: > + av_free(smk); > + > + return ret; > +} > + > +/** > + * Evalaute half of the cubic hermite interpolation expression, wrt one > interval endpoint > + * > + * @param x normalized input value at the endpoint > + * @param f output value at the endpoint > + * @param d derivative at the endpoint: normalized to the interval, and > properly sign adjusted > + * @return half of the interpolated value > +*/ > +static inline double interp_cubic_hermite_half(const double x, const double > f, > + const double d) > +{ > + double x2 = x * x, x3 = x2 * x; > + return f * (3.0 * x2 - 2.0 * x3) + d * (x3 - x2); > +} > + > +/** > + * Prepare the lookup table by piecewise monotonic cubic interpolation > (PCHIP) > + * > + * @param log_ctx for logging > + * @param y output lookup table (output) > + * @param points user-defined control points/endpoints > + * @param nbits bitdepth > + * @return 0 success > + * > + * References: > + * [1] F. N. Fritsch and J. Butland, A method for constructing local > monotone piecewise > + * cubic interpolants, SIAM J. Sci. Comput., 5(2), 300-304 (1984). > DOI:10.1137/0905021. > + * [2] scipy.interpolate: > https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.PchipInterpolator.html > +*/ > +static inline int interpolate_pchip(void *log_ctx, uint16_t *y, > + const struct keypoint *points, int > nbits) > +{ > + int i, ret = 0; > + const struct keypoint *point = points; > + const int lut_size = 1<<nbits; > + const int n = get_nb_points(points); // number of endpoints > + double *xi, *fi, *di, *hi, *mi; > + const int scale = lut_size - 1; // white value > + uint16_t x, x0; /* input index/value */ > + > + /* no change for n = 0 or 1 */ > + if (n == 0) { > + /* no points, no change */ > + for (i = 0; i < lut_size; ++i) y[i] = i; > + return 0; > + } > + > + if (n == 1) { > + /* 1 point - 1 color everywhere */ > + const uint16_t yval = CLIP(point->y * scale); > + for (i = 0; i < lut_size; ++i) y[i] = yval; > + return 0; > + } > + > + xi = av_calloc(3*n + 2*(n-1), sizeof(double)); /* output values at > inteval endpoints */ > + > + if (!xi) { > + ret = AVERROR(ENOMEM); > + goto end; > + } > + > + fi = xi + n; /* output values at inteval endpoints */ > + di = fi + n; /* output slope wrt normalized input at interval > endpoints */ > + hi = di + n; /* interval widths */ > + mi = hi + n - 1; /* linear slope over intervals */ > + > + /* scale endpoints and store them in a contiguous memory block */ > + for (i = 0; i < n; ++i) { > + xi[i] = point->x * scale; > + fi[i] = point->y * scale; > + point = point->next; > + } > + > + /* h(i) = x(i+1) - x(i); mi(i) = (f(i+1)-f(i))/h(i) */ > + for (i = 0; i < n - 1; ++i) { > + const double val = (xi[i+1]-xi[i]); > + hi[i] = val; > + mi[i] = (fi[i+1]-fi[i]) / val; > + } > + > + if (n == 2) { > + /* edge case, use linear interpolation */ > + const double m = mi[0], b = fi[0] - xi[0]*m; > + for (i = 0; i < lut_size; ++i) y[i] = CLIP((i*m + b)); > + goto end; > + } > + > + /* compute the derivatives at the endpoints*/ > + ret = pchip_find_derivatives(n-1,hi,mi,di); > + if (ret) goto end; > + > + /* interpolate/extrapolate */ > + x = 0; > + if (xi[0] > 0) { > + /* below first endpoint, use the first endpoint value*/ > + const double xi0 = xi[0]; > + const uint16_t yval = CLIP(fi[0]); > + for (; x < xi0; ++x) y[x] = yval; > + av_log(log_ctx, AV_LOG_DEBUG, "Interval -1: [0, %d] -> %d\n", x - > 1, yval); > + } > + > + /* for each interval */ > + for (i = 0, x0 = x; i < n-1; ++i, x0 = x) { > + > + const double xi0 = xi[i]; /* start-of-interval input value */ > + const double xi1 = xi[i + 1]; /* end-of-interval input value */ > + const double h = hi[i]; /* interval width */ > + const double f0 = fi[i]; /* start-of-interval output value */ > + const double f1 = fi[i + 1]; /* end-of-interval output value */ > + const double d0 = di[i]; /* start-of-interval derivative */ > + const double d1 = di[i + 1]; /* end-of-interval derivative */ > + > + /* fill the lut over the interval */ > + for (; x < xi1; ++x) { /* safe not to check j < lut_size */ > + const double xx = (x - xi0) / h; /* normalize input */ > + const double yy = interp_cubic_hermite_half(1 - xx, f0, -h * > d0) > + + interp_cubic_hermite_half(xx, f1, h * d1); > + y[x] = CLIP(yy); > + } > + > + if (x > x0) > + av_log(log_ctx, AV_LOG_DEBUG, "Interval %d: [%d, %d] -> [%d, > %d]\n", > + i, x0, x-1, y[x0], > y[x-1]); > + else > + av_log(log_ctx, AV_LOG_DEBUG, "Interval %d: empty\n", i); > + } > + > + if (x < lut_size) { > + /* above the last endpoint, use the last endpoint value*/ > + const uint16_t yval = CLIP(fi[n - 1]); > + av_log(log_ctx, AV_LOG_DEBUG, "Interval %d: [%d, %d] -> %d\n", > + n, x, lut_size - 1, yval); > + for (; x < lut_size; ++x) y[x] = yval; > + } > + > +end: > + av_free(xi); > + return ret; > } > > -DECLARE_INTERPOLATE_FUNC(8) > -DECLARE_INTERPOLATE_FUNC(9) > -DECLARE_INTERPOLATE_FUNC(10) > -DECLARE_INTERPOLATE_FUNC(12) > -DECLARE_INTERPOLATE_FUNC(14) > -DECLARE_INTERPOLATE_FUNC(16) > > static int parse_psfile(AVFilterContext *ctx, const char *fname) > { > @@ -651,14 +870,10 @@ static int config_input(AVFilterLink *inlink) > ret = parse_points_str(ctx, comp_points + i, > curves->comp_points_str[i], curves->lut_size); > if (ret < 0) > return ret; > - switch (curves->depth) { > - case 8: ret = interpolate8 (ctx, curves->graph[i], > comp_points[i]); break; > - case 9: ret = interpolate9 (ctx, curves->graph[i], > comp_points[i]); break; > - case 10: ret = interpolate10(ctx, curves->graph[i], > comp_points[i]); break; > - case 12: ret = interpolate12(ctx, curves->graph[i], > comp_points[i]); break; > - case 14: ret = interpolate14(ctx, curves->graph[i], > comp_points[i]); break; > - case 16: ret = interpolate16(ctx, curves->graph[i], > comp_points[i]); break; > - } > + if (curves->interp==INTERP_PCHIP) > + ret = interpolate_pchip (ctx, curves->graph[i], comp_points[i], > curves->depth); > + else > + ret = interpolate (ctx, curves->graph[i], comp_points[i], > curves->depth); > if (ret < 0) > return ret; > } > @@ -735,7 +950,7 @@ static int process_command(AVFilterContext *ctx, const > char *cmd, const char *ar > > if (!strcmp(cmd, "plot")) { > curves->saved_plot = 0; > - } else if (!strcmp(cmd, "all") || !strcmp(cmd, "preset") || > !strcmp(cmd, "psfile")) { > + } else if (!strcmp(cmd, "all") || !strcmp(cmd, "preset") || > !strcmp(cmd, "psfile") || !strcmp(cmd, "interp")) { > if (!strcmp(cmd, "psfile")) > curves->parsed_psfile = 0; > av_freep(&curves->comp_points_str_all); > -- > 2.25.1 > > _______________________________________________ > ffmpeg-devel mailing list > ffmpeg-devel@ffmpeg.org > https://ffmpeg.org/mailman/listinfo/ffmpeg-devel > > To unsubscribe, visit link above, or email > ffmpeg-devel-request@ffmpeg.org with subject "unsubscribe". >
diff --git a/doc/filters.texi b/doc/filters.texi index d0f718678c..08a79644e1 100644 --- a/doc/filters.texi +++ b/doc/filters.texi @@ -10389,11 +10389,15 @@ By default, a component curve is defined by the two points @var{(0;0)} and "adjusted" to its own value, which means no change to the image. The filter allows you to redefine these two points and add some more. A new -curve (using a natural cubic spline interpolation) will be define to pass -smoothly through all these new coordinates. The new defined points needs to be -strictly increasing over the x-axis, and their @var{x} and @var{y} values must -be in the @var{[0;1]} interval. If the computed curves happened to go outside -the vector spaces, the values will be clipped accordingly. +curve will be define to pass smoothly through all these new coordinates. The +new defined points needs to be strictly increasing over the x-axis, and their +@var{x} and @var{y} values must be in the @var{[0;1]} interval. The curve is +formed by using a natural or monotonic cubic spline interpolation, depending +on the @var{interp} option (default: @code{natural}). The @code{natural} +spline produces a smoother curve in general while the monotonic (@code{pchip}) +spline guarantees the transitions between the specified points to be +monotonic. If the computed curves happened to go outside the vector spaces, +the values will be clipped accordingly. The filter accepts the following options: @@ -10437,6 +10441,15 @@ options. In this case, the unset component(s) will fallback on this Specify a Photoshop curves file (@code{.acv}) to import the settings from. @item plot Save Gnuplot script of the curves in specified file. +@item interp +Specify the kind of interpolation. Available algorithms are: +@table @samp +@item natural +Natural cubic spline using a piece-wise cubic polynomial that is twice continuously differentiable. +@item pchip +Monotonic cubic spline using a piecewise cubic Hermite interpolating polynomial (PCHIP). +@end table + @end table To avoid some filtergraph syntax conflicts, each key points list need to be diff --git a/libavfilter/vf_curves.c b/libavfilter/vf_curves.c index 498b06f6e5..d0efa380e1 100644 --- a/libavfilter/vf_curves.c +++ b/libavfilter/vf_curves.c @@ -58,6 +58,12 @@ enum preset { NB_PRESETS, }; +enum interp { + INTERP_NATURAL, + INTERP_PCHIP, + NB_INTERPS, +}; + typedef struct CurvesContext { const AVClass *class; int preset; @@ -73,6 +79,7 @@ typedef struct CurvesContext { int is_16bit; int depth; int parsed_psfile; + int interp; int (*filter_slice)(AVFilterContext *ctx, void *arg, int jobnr, int nb_jobs); } CurvesContext; @@ -107,6 +114,9 @@ static const AVOption curves_options[] = { { "all", "set points coordinates for all components", OFFSET(comp_points_str_all), AV_OPT_TYPE_STRING, {.str=NULL}, .flags = FLAGS }, { "psfile", "set Photoshop curves file name", OFFSET(psfile), AV_OPT_TYPE_STRING, {.str=NULL}, .flags = FLAGS }, { "plot", "save Gnuplot script of the curves in specified file", OFFSET(plot_filename), AV_OPT_TYPE_STRING, {.str=NULL}, .flags = FLAGS }, + { "interp", "specify the kind of interpolation", OFFSET(interp), AV_OPT_TYPE_INT, {.i64=INTERP_NATURAL}, INTERP_NATURAL, NB_INTERPS-1, FLAGS, "interp_name" }, + { "natural", "natural cubic spline", 0, AV_OPT_TYPE_CONST, {.i64=INTERP_NATURAL}, 0, 0, FLAGS, "interp_name" }, + { "pchip", "monotonically cubic interpolation", 0, AV_OPT_TYPE_CONST, {.i64=INTERP_PCHIP}, 0, 0, FLAGS, "interp_name" }, { NULL } }; @@ -336,21 +346,230 @@ end: av_free(h); av_free(r); return ret; + +} + +#define SIGN(x) (x>0.0?1:x<0.0?-1:0) + +/** + * Evalaute the derivative of an edge endpoint + * + * @param h0 input interval of the interval closest to the edge + * @param h1 input interval of the interval next to the closest + * @param m0 linear slope of the interval closest to the edge + * @param m1 linear slope of the intervalnext to the closest + * @return edge endpoint derivative + * + * Based on scipy.interpolate._edge_case() + * https://github.com/scipy/scipy/blob/2e5883ef7af4f5ed4a5b80a1759a45e43163bf3f/scipy/interpolate/_cubic.py#L239 + * which is a python implementation of the special case endpoints, as suggested in + * Cleve Moler, Numerical Computing with MATLAB, Chap 3.6 (pchiptx.m) +*/ +static double pchip_edge_case(double h0, double h1, double m0, double m1) +{ + int mask, mask2; + double d; + + d = ((2 * h0 + h1) * m0 - h0 * m1) / (h0 + h1); + + mask = SIGN(d) != SIGN(m0); + mask2 = (SIGN(m0) != SIGN(m1)) && (fabs(d) > 3. * fabs(m0)); + + if (mask) d = 0.0; + else if (mask2) d = 3.0 * m0; + + return d; } -#define DECLARE_INTERPOLATE_FUNC(nbits) \ -static int interpolate##nbits(void *log_ctx, uint16_t *y, \ - const struct keypoint *points) \ -{ \ - return interpolate(log_ctx, y, points, nbits); \ +/** + * Evalaute the piecewise polynomial derivatives at endpoints + * + * @param n input interval of the interval closest to the edge + * @param hk input intervals + * @param mk linear slopes over intervals + * @param dk endpoint derivatives (output) + * @return 0 success + * + * Based on scipy.interpolate._find_derivatives() + * https://github.com/scipy/scipy/blob/2e5883ef7af4f5ed4a5b80a1759a45e43163bf3f/scipy/interpolate/_cubic.py#L254 +*/ + +static int pchip_find_derivatives(const int n, const double *hk, const double *mk, double *dk) +{ + int ret = 0; + const int m = n - 1; + int8_t *smk; + + smk = av_malloc(n); + if (!smk) { + ret = AVERROR(ENOMEM); + goto end; + } + + /* smk = sgn(mk) */ + for (int i = 0; i < n; ++i) smk[i] = SIGN(mk[i]); + + /* check the strict monotonicity */ + for (int i = 0; i < m; ++i) { + int8_t condition = (smk[i + 1] != smk[i]) || (mk[i + 1] == 0) || (mk[i] == 0); + if (condition) { + dk[i + 1] = 0.0; + } else { + double w1 = 2 * hk[i + 1] + hk[i]; + double w2 = hk[i + 1] + 2 * hk[i]; + dk[i + 1] = (w1 + w2) / (w1 / mk[i] + w2 / mk[i + 1]); + } + } + + dk[0] = pchip_edge_case(hk[0], hk[1], mk[0], mk[1]); + dk[n] = pchip_edge_case(hk[n - 1], hk[n - 2], mk[n - 1], mk[n - 2]); + +end: + av_free(smk); + + return ret; +} + +/** + * Evalaute half of the cubic hermite interpolation expression, wrt one interval endpoint + * + * @param x normalized input value at the endpoint + * @param f output value at the endpoint + * @param d derivative at the endpoint: normalized to the interval, and properly sign adjusted + * @return half of the interpolated value +*/ +static inline double interp_cubic_hermite_half(const double x, const double f, + const double d) +{ + double x2 = x * x, x3 = x2 * x; + return f * (3.0 * x2 - 2.0 * x3) + d * (x3 - x2); +} + +/** + * Prepare the lookup table by piecewise monotonic cubic interpolation (PCHIP) + * + * @param log_ctx for logging + * @param y output lookup table (output) + * @param points user-defined control points/endpoints + * @param nbits bitdepth + * @return 0 success + * + * References: + * [1] F. N. Fritsch and J. Butland, A method for constructing local monotone piecewise + * cubic interpolants, SIAM J. Sci. Comput., 5(2), 300-304 (1984). DOI:10.1137/0905021. + * [2] scipy.interpolate: https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.PchipInterpolator.html +*/ +static inline int interpolate_pchip(void *log_ctx, uint16_t *y, + const struct keypoint *points, int nbits) +{ + int i, ret = 0; + const struct keypoint *point = points; + const int lut_size = 1<<nbits; + const int n = get_nb_points(points); // number of endpoints + double *xi, *fi, *di, *hi, *mi; + const int scale = lut_size - 1; // white value + uint16_t x, x0; /* input index/value */ + + /* no change for n = 0 or 1 */ + if (n == 0) { + /* no points, no change */ + for (i = 0; i < lut_size; ++i) y[i] = i; + return 0; + } + + if (n == 1) { + /* 1 point - 1 color everywhere */ + const uint16_t yval = CLIP(point->y * scale); + for (i = 0; i < lut_size; ++i) y[i] = yval; + return 0; + } + + xi = av_calloc(3*n + 2*(n-1), sizeof(double)); /* output values at inteval endpoints */ + + if (!xi) { + ret = AVERROR(ENOMEM); + goto end; + } + + fi = xi + n; /* output values at inteval endpoints */ + di = fi + n; /* output slope wrt normalized input at interval endpoints */ + hi = di + n; /* interval widths */ + mi = hi + n - 1; /* linear slope over intervals */ + + /* scale endpoints and store them in a contiguous memory block */ + for (i = 0; i < n; ++i) { + xi[i] = point->x * scale; + fi[i] = point->y * scale; + point = point->next; + } + + /* h(i) = x(i+1) - x(i); mi(i) = (f(i+1)-f(i))/h(i) */ + for (i = 0; i < n - 1; ++i) { + const double val = (xi[i+1]-xi[i]); + hi[i] = val; + mi[i] = (fi[i+1]-fi[i]) / val; + } + + if (n == 2) { + /* edge case, use linear interpolation */ + const double m = mi[0], b = fi[0] - xi[0]*m; + for (i = 0; i < lut_size; ++i) y[i] = CLIP((i*m + b)); + goto end; + } + + /* compute the derivatives at the endpoints*/ + ret = pchip_find_derivatives(n-1,hi,mi,di); + if (ret) goto end; + + /* interpolate/extrapolate */ + x = 0; + if (xi[0] > 0) { + /* below first endpoint, use the first endpoint value*/ + const double xi0 = xi[0]; + const uint16_t yval = CLIP(fi[0]); + for (; x < xi0; ++x) y[x] = yval; + av_log(log_ctx, AV_LOG_DEBUG, "Interval -1: [0, %d] -> %d\n", x - 1, yval); + } + + /* for each interval */ + for (i = 0, x0 = x; i < n-1; ++i, x0 = x) { + + const double xi0 = xi[i]; /* start-of-interval input value */ + const double xi1 = xi[i + 1]; /* end-of-interval input value */ + const double h = hi[i]; /* interval width */ + const double f0 = fi[i]; /* start-of-interval output value */ + const double f1 = fi[i + 1]; /* end-of-interval output value */ + const double d0 = di[i]; /* start-of-interval derivative */ + const double d1 = di[i + 1]; /* end-of-interval derivative */ + + /* fill the lut over the interval */ + for (; x < xi1; ++x) { /* safe not to check j < lut_size */ + const double xx = (x - xi0) / h; /* normalize input */ + const double yy = interp_cubic_hermite_half(1 - xx, f0, -h * d0) + + interp_cubic_hermite_half(xx, f1, h * d1); + y[x] = CLIP(yy); + } + + if (x > x0) + av_log(log_ctx, AV_LOG_DEBUG, "Interval %d: [%d, %d] -> [%d, %d]\n", + i, x0, x-1, y[x0], y[x-1]); + else + av_log(log_ctx, AV_LOG_DEBUG, "Interval %d: empty\n", i); + } + + if (x < lut_size) { + /* above the last endpoint, use the last endpoint value*/ + const uint16_t yval = CLIP(fi[n - 1]); + av_log(log_ctx, AV_LOG_DEBUG, "Interval %d: [%d, %d] -> %d\n", + n, x, lut_size - 1, yval); + for (; x < lut_size; ++x) y[x] = yval; + } + +end: + av_free(xi); + return ret; } -DECLARE_INTERPOLATE_FUNC(8) -DECLARE_INTERPOLATE_FUNC(9) -DECLARE_INTERPOLATE_FUNC(10) -DECLARE_INTERPOLATE_FUNC(12) -DECLARE_INTERPOLATE_FUNC(14) -DECLARE_INTERPOLATE_FUNC(16) static int parse_psfile(AVFilterContext *ctx, const char *fname) { @@ -651,14 +870,10 @@ static int config_input(AVFilterLink *inlink) ret = parse_points_str(ctx, comp_points + i, curves->comp_points_str[i], curves->lut_size); if (ret < 0) return ret; - switch (curves->depth) { - case 8: ret = interpolate8 (ctx, curves->graph[i], comp_points[i]); break; - case 9: ret = interpolate9 (ctx, curves->graph[i], comp_points[i]); break; - case 10: ret = interpolate10(ctx, curves->graph[i], comp_points[i]); break; - case 12: ret = interpolate12(ctx, curves->graph[i], comp_points[i]); break; - case 14: ret = interpolate14(ctx, curves->graph[i], comp_points[i]); break; - case 16: ret = interpolate16(ctx, curves->graph[i], comp_points[i]); break; - } + if (curves->interp==INTERP_PCHIP) + ret = interpolate_pchip (ctx, curves->graph[i], comp_points[i], curves->depth); + else + ret = interpolate (ctx, curves->graph[i], comp_points[i], curves->depth); if (ret < 0) return ret; } @@ -735,7 +950,7 @@ static int process_command(AVFilterContext *ctx, const char *cmd, const char *ar if (!strcmp(cmd, "plot")) { curves->saved_plot = 0; - } else if (!strcmp(cmd, "all") || !strcmp(cmd, "preset") || !strcmp(cmd, "psfile")) { + } else if (!strcmp(cmd, "all") || !strcmp(cmd, "preset") || !strcmp(cmd, "psfile") || !strcmp(cmd, "interp")) { if (!strcmp(cmd, "psfile")) curves->parsed_psfile = 0; av_freep(&curves->comp_points_str_all);
summary: This patch modifies the `curves` filter with new `interp` option to let user pick the existing natural cubic spline interpolation and the new PCHIP interapolation. reason: The natural cubic spline does not impose monotonicity between the keypoints. As such, the fitted curve may vary wildly against user's intension. The PCHIP interpolation is not as smooth as the natural spline but guarantees the monotonicity. Providing both options enhances users experience (e.g., reduces the number of keypoints to realize the desired curve). See the related bug report for the example of an ill-interpolated curve. alternate solution: Both Photoshop and GIMP appear to use monotonic interpolation in their curve tools, which were the models for this filter. As such, an alternate solution is to drop the natural spline and go without the `interp` option. related bug report: https://trac.ffmpeg.org/ticket/9947 (filed by myself) Signed-off-by: Takeshi (Kesh) Ikuma <tikuma@hotmail.com> modified: doc/filters.texi modified: libavfilter/vf_curves.c --- doc/filters.texi | 23 +++- libavfilter/vf_curves.c | 255 ++++++++++++++++++++++++++++++++++++---- 2 files changed, 253 insertions(+), 25 deletions(-)