From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: from ffbox0-bg.mplayerhq.hu (ffbox0-bg.ffmpeg.org [79.124.17.100]) by master.gitmailbox.com (Postfix) with ESMTP id C576E4317E for ; Sat, 22 Oct 2022 07:50:46 +0000 (UTC) Received: from [127.0.1.1] (localhost [127.0.0.1]) by ffbox0-bg.mplayerhq.hu (Postfix) with ESMTP id B122D68BE1A; Sat, 22 Oct 2022 10:50:43 +0300 (EEST) Received: from mail-ua1-f41.google.com (mail-ua1-f41.google.com [209.85.222.41]) by ffbox0-bg.mplayerhq.hu (Postfix) with ESMTPS id 8B61268BCDC for ; Sat, 22 Oct 2022 10:50:37 +0300 (EEST) Received: by mail-ua1-f41.google.com with SMTP id p4so3156597uao.0 for ; Sat, 22 Oct 2022 00:50:37 -0700 (PDT) DKIM-Signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=20210112; h=cc:to:subject:message-id:date:from:references:in-reply-to :mime-version:from:to:cc:subject:date:message-id:reply-to; bh=F5o+GDRXQKUyrhji9V+PY9DeDSZQ4mq6H+/vwZFyo4I=; b=AiXKBr4PAv/L1TV6+i818m9eFbpvafOD2wP3DczWtgoQ8uFlOKh55ddEyzAj2Qrhsp cCZ5eq3oFIj/qfcOvILt3ww2IgjnWiiZczoK3Aq8sf8zL/JMDqEqgDTk/azrBM8efWTD gMo4TuazYnZB6wJLF83Hdw1NdsqBCVhV3Y/cUcO43T0rrqJkKiJIDxN7yATsdbM/er8x tTJO+834wbmEmHV6lrLstkT9oXhZ4IeJeIJf4C/T5uHLw5NJuH/ZvCChXG0dHLDePSbX 6ZgL9vo1t/qvFlV3pGZDZ6hnqXE3uP0dj3hPdDvam03QI94//uAFnNYCFySIJmjIe4Ve vZvA== X-Google-DKIM-Signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=1e100.net; s=20210112; h=cc:to:subject:message-id:date:from:references:in-reply-to :mime-version:x-gm-message-state:from:to:cc:subject:date:message-id :reply-to; bh=F5o+GDRXQKUyrhji9V+PY9DeDSZQ4mq6H+/vwZFyo4I=; b=jIKESnuCnk16wNf5TsmZ6NNVywr2rllQg26OJN9nY72Z0ZR77uiravKSJqfAZOGiQB ErGz5ywrVzdBmvycRMujKjz7bOeckB2gtMP3r3FvQEpvUr0F3Phc3UTOqE121bTodaSD vOKmOTGWi76aUSEEc7zntW90P17bhdavf/QwtNcX9NWYYKCTycJKVUvCxoSQwLbNpDAc jnXaC8LkARo0xKcJMfmhIJz515kR1DKP50pUAjSENYyrG5/YZ5JTA2ebB3IVxyqbHQ4T j2u/b/V7IOt1t+OiqYo93DY+B03ROAb/QMMDGBcQksflOg0WJ6p0ALOBMv1StYDXpRm3 ly2A== X-Gm-Message-State: ACrzQf0bSfoxmTKsQUJJ67/EydPvIjyhQHZ682sPQo/QUN51VBbvMYkw AZD1zCjShNxl0PrJDKtZTNvr29OqmUHEmR3n3NBRP7Px X-Google-Smtp-Source: AMsMyM6Zx1WhHaIAIaKGgIi5myqefOOajW9VykwpyMmcm8bWObNGveD8arWysYrp+yp92q00G3wOvXAo8xfpAZicFOk= X-Received: by 2002:ab0:1e8d:0:b0:3d9:6cb0:4a87 with SMTP id o13-20020ab01e8d000000b003d96cb04a87mr15108668uak.42.1666425035850; Sat, 22 Oct 2022 00:50:35 -0700 (PDT) MIME-Version: 1.0 Received: by 2002:a05:612c:612:b0:314:ac6a:1eb7 with HTTP; Sat, 22 Oct 2022 00:50:34 -0700 (PDT) In-Reply-To: <20221003014404.749-1-tikuma@hotmail.com> References: <20221003014404.749-1-tikuma@hotmail.com> From: Paul B Mahol Date: Sat, 22 Oct 2022 09:50:34 +0200 Message-ID: To: FFmpeg development discussions and patches Subject: Re: [FFmpeg-devel] [PATCH] avfilter/vf_curves: add PCHIP interpolator and interp option X-BeenThere: ffmpeg-devel@ffmpeg.org X-Mailman-Version: 2.1.29 Precedence: list List-Id: FFmpeg development discussions and patches List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , Reply-To: FFmpeg development discussions and patches Cc: "Takeshi \(Kesh\) Ikuma" Content-Type: text/plain; charset="us-ascii" Content-Transfer-Encoding: 7bit Errors-To: ffmpeg-devel-bounces@ffmpeg.org Sender: "ffmpeg-devel" Archived-At: List-Archive: List-Post: On 10/3/22, Takeshi (Kesh) Ikuma 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 > > 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< + 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".