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 41A0044247 for ; Mon, 3 Oct 2022 01:44:19 +0000 (UTC) Received: from [127.0.1.1] (localhost [127.0.0.1]) by ffbox0-bg.mplayerhq.hu (Postfix) with ESMTP id 8AE0868BA58; Mon, 3 Oct 2022 04:44:17 +0300 (EEST) Received: from mail-oo1-f52.google.com (mail-oo1-f52.google.com [209.85.161.52]) by ffbox0-bg.mplayerhq.hu (Postfix) with ESMTPS id 0F67868B7C4 for ; Mon, 3 Oct 2022 04:44:11 +0300 (EEST) Received: by mail-oo1-f52.google.com with SMTP id h1-20020a4aa741000000b004756c611188so5892097oom.4 for ; Sun, 02 Oct 2022 18:44:10 -0700 (PDT) DKIM-Signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=20210112; h=content-transfer-encoding:mime-version:message-id:date:subject:cc :to:from:from:to:cc:subject:date; bh=PtMN6ExMDUCGfQOavqgp5i2CxfeZ70rarmOb9qG8Fb4=; b=cl5j6c5CjciiOs1wv1mBOG5nYLE0EK+OnMArjiXF7C+cEsAY8IwQEHzBJNcLtv2e+V XsmGVavSLGPtUiqgmqk1J08QdIEdbBXqFwedVwx1f6Fvz4BVcRksVTvlI70XsHPXqze3 QALxjgtr8+LkUOYJzxNXb5EwDycycWSlxW9AK24ogPaiHdF6YdcDnexx/8pSm5jaz/r9 RD5AKd/U6rKE7tAdnXUEeR6+8eCvPwf3JwSoV+I/dn7hPf4pqk/3lT2oEWWe7hChlhYr wvR+TESTKqZb8iZoJkYZCnmD13Tvk6T+DqzFiP6Jk+Md0dkZJN2N01dXoXAofenRXPNh cQ1g== X-Google-DKIM-Signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=1e100.net; s=20210112; h=content-transfer-encoding:mime-version:message-id:date:subject:cc :to:from:x-gm-message-state:from:to:cc:subject:date; bh=PtMN6ExMDUCGfQOavqgp5i2CxfeZ70rarmOb9qG8Fb4=; b=IK4GAGb1mTpgkiER8jgfj7Be2YmnIU5ERuXT5Nb3RRaESguL9Koi116/GMKmO8YyjT hgDZ+qrp6i4VHmUyYmxIgR/ZP8aAT2WC6GloFQLBLGudEcCLJPtdOkK5FQRR6xpqvXjf jMMW+6HNuoRxmp2xkfIuRI3DIDkqt1F7pCqvSixy4PECO7Ct8DiqWXGoUyLwb9S+NNRZ w1MS/0C7UrgSyorKytgtIblemxqg9Rxu/exTCsInoN1JJP6buZAcEmGPfpLZyKk8yIA4 CwBmIerVIM16Se/TqIgWcy2bpsHHXG1KmU3JckmbCMDr14I+mB71oZvBd4OqiNpUCh4v cUsA== X-Gm-Message-State: ACrzQf1upLgUM8Pc+WJOS3+feetNWU3vOJBGMk4n4pavmBEewBNDmjsi z2Anuo2NJ7q44UUa3kO6pJsnIqM3Tkc= X-Google-Smtp-Source: AMsMyM6lAA5jXfg7O8mG/dsncOIxlIP4EZER1wtYfGpBCuZC0v1ODn+vskQMgMkQgzhUyvK6/xXULQ== X-Received: by 2002:a4a:4113:0:b0:476:6fbc:f9b7 with SMTP id x19-20020a4a4113000000b004766fbcf9b7mr6905199ooa.47.1664761448753; Sun, 02 Oct 2022 18:44:08 -0700 (PDT) Received: from Excalibur.localdomain ([68.72.114.64]) by smtp.gmail.com with ESMTPSA id er30-20020a056870c89e00b0012b298699dbsm2583965oab.1.2022.10.02.18.44.07 (version=TLS1_3 cipher=TLS_AES_256_GCM_SHA384 bits=256/256); Sun, 02 Oct 2022 18:44:08 -0700 (PDT) From: "Takeshi (Kesh) Ikuma" X-Google-Original-From: "Takeshi (Kesh) Ikuma" To: ffmpeg-devel@ffmpeg.org Date: Sun, 2 Oct 2022 20:44:03 -0500 Message-Id: <20221003014404.749-1-tikuma@hotmail.com> X-Mailer: git-send-email 2.25.1 MIME-Version: 1.0 Subject: [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: 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(-) 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<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".