Re: [FFmpeg-devel] [PATCH] avfilter/vf_curves: add PCHIP interpolator and interp option

[Paul B Mahol <onemda@gmail.com>]

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".
>
_______________________________________________
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".