From patchwork Mon Oct 3 01:44:03 2022 Content-Type: text/plain; charset="utf-8" MIME-Version: 1.0 Content-Transfer-Encoding: 7bit X-Patchwork-Submitter: "Takeshi (Kesh) Ikuma" X-Patchwork-Id: 38531 Delivered-To: ffmpegpatchwork2@gmail.com Received: by 2002:a05:6a20:3b1c:b0:96:9ee8:5cfd with SMTP id c28csp1044019pzh; Sun, 2 Oct 2022 18:44:22 -0700 (PDT) X-Google-Smtp-Source: AMsMyM6WZ1wAb63di4a56DIemD3ldDMRT8xuUetq2Q6oxphjXouOsD0+g9D90o8oT4zhYYVGWN9p X-Received: by 2002:aa7:dc10:0:b0:440:b446:c0cc with SMTP id b16-20020aa7dc10000000b00440b446c0ccmr16845221edu.34.1664761461882; Sun, 02 Oct 2022 18:44:21 -0700 (PDT) ARC-Seal: i=1; a=rsa-sha256; t=1664761461; cv=none; d=google.com; s=arc-20160816; b=UEaC1EC1vi90nesWBNRZT3ORyzpoNkxmktalYLGLhafIEkowmmEV+g7N8/ofa9YpXu ismY/rXMB7880Hn9fWaZKks0xvKKF7igq4GRWpoXMSCpxMR5oTrMOHXm5KkhajWI9O/4 cGjpgwZ68OrwRga2mWiqSj95UzmefC8m/9Uz+a7jmPS4MYu4JfuOXwN7yHreexYSmrSd ajsU0W6WA19DhPZDpW+sooMaUtIInrPORVPITY2Ry9xH0WRl3XAjc/wjq4/zRIyIWS45 66XuaaRTTa9rtjx97LAUceSe1OmhHyvU7SYS308iI11vSCdmLT95ogsQkIk6HsehzVyf YUeg== ARC-Message-Signature: i=1; a=rsa-sha256; c=relaxed/relaxed; d=google.com; s=arc-20160816; h=sender:errors-to:content-transfer-encoding:cc:reply-to :list-subscribe:list-help:list-post:list-archive:list-unsubscribe :list-id:precedence:subject:mime-version:message-id:date:to:from :dkim-signature:delivered-to; bh=xB45uESYk2X15M/pHJVTG+T4EfClMBVvjl467YiWfhw=; b=bo9n1Y5gZAZbyE7SjCj5vPK65ykRUeEQuNP8C8DYBgxw8oVZNFdGnrgE9+ZaYN3x9M PXoj+x0LVgpSkBuylp2Ph+egQWlu1dqapYU4Zgxhdy+AiaTzbRdtGUkVVhcN044/Nxa2 bDdAz/7fHzRQoaT/TEj/3TiqGQKcan8zsd2MmPcD7Ht/J7OCWi9GzbmOeM0rpDQ/McyQ QQgGLeEJoa5aEtxQCCCYvmP1yM2NGwFfvRk8zk84TVUVqY/AAVFqSHBDdvjPo4/r+RXm ytcblbbgQqmWwdB99eC6/RwpL5JVGOlrmUf21BRQeiMKSJwaYGbnyOkvZ5ZkHaEh5uSv o4lA== ARC-Authentication-Results: i=1; mx.google.com; dkim=neutral (body hash did not verify) header.i=@gmail.com header.s=20210112 header.b=cl5j6c5C; spf=pass (google.com: domain of ffmpeg-devel-bounces@ffmpeg.org designates 79.124.17.100 as permitted sender) smtp.mailfrom=ffmpeg-devel-bounces@ffmpeg.org; dmarc=fail (p=NONE sp=QUARANTINE dis=NONE) header.from=gmail.com Return-Path: Received: from ffbox0-bg.mplayerhq.hu (ffbox0-bg.ffmpeg.org. [79.124.17.100]) by mx.google.com with ESMTP id q7-20020a056402518700b00457c4773001si8012597edd.131.2022.10.02.18.44.21; Sun, 02 Oct 2022 18:44:21 -0700 (PDT) Received-SPF: pass (google.com: domain of ffmpeg-devel-bounces@ffmpeg.org designates 79.124.17.100 as permitted sender) client-ip=79.124.17.100; Authentication-Results: mx.google.com; dkim=neutral (body hash did not verify) header.i=@gmail.com header.s=20210112 header.b=cl5j6c5C; spf=pass (google.com: domain of ffmpeg-devel-bounces@ffmpeg.org designates 79.124.17.100 as permitted sender) smtp.mailfrom=ffmpeg-devel-bounces@ffmpeg.org; dmarc=fail (p=NONE sp=QUARANTINE dis=NONE) header.from=gmail.com Received: from [127.0.1.1] (localhost [127.0.0.1]) by ffbox0-bg.mplayerhq.hu (Postfix) with ESMTP id C594968ABFD; Mon, 3 Oct 2022 04:44:17 +0300 (EEST) X-Original-To: ffmpeg-devel@ffmpeg.org Delivered-To: ffmpeg-devel@ffmpeg.org 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-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" Errors-To: ffmpeg-devel-bounces@ffmpeg.org Sender: "ffmpeg-devel" X-TUID: 08nKC+l/7Aw5 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);