Chiebot-Mirror
/
FFmpeg
mirror of https://github.com/FFmpeg/FFmpeg.git
8
0
Fork 0
Code Issues Projects Releases Wiki Activity

avfilter/vf_curves: add PCHIP interpolator and interp option

Browse Source 
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>
pull/388/head
Takeshi (Kesh) Ikuma 2 years ago committed by Paul B Mahol
parent 993e3a6b54
commit 65f96a965a
2 changed files with 262 additions and 25 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 23
      doc/filters.texi
  2. 264
      libavfilter/vf_curves.c

23
doc/filters.texi
Unescape View File

@ -10404,11 +10404,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:
@ -10452,6 +10456,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

264
libavfilter/vf_curves.c
Unescape View File

@ -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,239 @@ end:
av_free(h);
av_free(r);
return ret;
}
#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); \
#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;
}
/**
* 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)
{
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; /* input index/value */
int ret = 0;
/* no change for n = 0 or 1 */
if (n == 0) {
/* no points, no change */
for (int 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 (int i = 0; i < lut_size; i++) y[i] = yval;
return 0;
}
xi = av_calloc(3*n + 2*(n-1), sizeof(double)); /* output values at interval 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 (int 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 (int 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 (int 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 double yi0 = fi[0];
const uint16_t yval = CLIP(yi0);
for (; x < xi0; x++) {
y[x] = yval;
av_log(log_ctx, AV_LOG_TRACE, "f(%f)=%f -> y[%d]=%d\n", xi0, yi0, x, y[x]);
}
av_log(log_ctx, AV_LOG_DEBUG, "Interval -1: [0, %d] -> %d\n", x - 1, yval);
}
/* for each interval */
for (int 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);
av_log(log_ctx, AV_LOG_TRACE, "f(%f)=%f -> y[%d]=%d\n", xx, yy, x, y[x]);
}
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 && x < lut_size) {
/* above the last endpoint, use the last endpoint value*/
const double xi1 = xi[n - 1];
const double yi1 = fi[n - 1];
const uint16_t yval = CLIP(yi1);
av_log(log_ctx, AV_LOG_DEBUG, "Interval %d: [%d, %d] -> %d\n",
n-1, x, lut_size - 1, yval);
for (; x && x < lut_size; x++) { /* loop until int overflow */
y[x] = yval;
av_log(log_ctx, AV_LOG_TRACE, "f(%f)=%f -> y[%d]=%d\n", xi1, yi1, 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 +879,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 +959,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);

Write Preview
Loading…
Cancel
Save