From patchwork Fri Jul 28 10:53:05 2017 Content-Type: text/plain; charset="utf-8" MIME-Version: 1.0 Content-Transfer-Encoding: 7bit X-Patchwork-Submitter: Ashish Singh X-Patchwork-Id: 4485 Delivered-To: ffmpegpatchwork@gmail.com Received: by 10.103.1.76 with SMTP id 73csp239975vsb; Fri, 28 Jul 2017 04:22:57 -0700 (PDT) X-Received: by 10.28.30.138 with SMTP id e132mr5547267wme.43.1501240976121; Fri, 28 Jul 2017 04:22:56 -0700 (PDT) ARC-Seal: i=1; a=rsa-sha256; t=1501240976; cv=none; d=google.com; s=arc-20160816; b=i7Csq1nWwniev3fIl+wkNGjKJ9qIhjDsP89164K/UutvAswj4aBDUmB7wvIbE9CGO0 AOs1g+58CCbNb4d79saiZGQox7Ol4flnn2wTUezojvo0aUU5f4Bvhj47a4RZMgyFBHlv z+YSONzTQ2Qb9NpMy1FOYq6VyYDfwTLq9GrRUAZ+/mMUnWuoJDeEGdO5YpArVil3D1oG GJL4TFp85DL8u5PN5Yl6WMYhnvyQzQOO01NoQdzXoBt2Neb8ac2Ae3UihTN48Ys0Bh8j lCMFfmERcaHHP7m0CqpxhtYiE7QAYIrjUJY1wuZ4HE/XY8KK8S/A2AvPMtR9KGgfvULr F6ZQ== ARC-Message-Signature: i=1; a=rsa-sha256; c=relaxed/relaxed; d=google.com; s=arc-20160816; h=sender:errors-to:content-transfer-encoding:mime-version:cc:reply-to :list-subscribe:list-help:list-post:list-archive:list-unsubscribe :list-id:precedence:subject:references:in-reply-to:message-id:date :to:from:dkim-signature:delivered-to:arc-authentication-results; bh=jnIWLXHzGcDrBGmhQmzOZA1h2Do4wOfRT1WeKhr+grU=; b=RSb1oUbpoJw/blpl4/qunq+fTScn52bhN80uAZ8QLviBSnnptUBV3//R+g3CNULps+ f+rJ7+5KQvtk28+rlfm/5k1hlpK+zquQpuvJ6lNZ8Ma3D1qeh1xDV27wlFcBW4EiYClb Dwa/EJRBIPBasi6JTxTMHZAyxHXO0h82kgAnXqZ8yN6w7MtGwkS558Sf/LN4NqKBbG+v ApbfwFCaHwob3BrhuGWnYa0Y/uWWVsifhJTXdxZRDPJbKg1Mcu/fozQjiYqN5fA4n+r9 0SuCZI2OGGhFGV2zQhnTZY81W+heHaCsmeu/abjox56n/xCIgpigllqsG26gY5du1cP9 Jy+w== ARC-Authentication-Results: i=1; mx.google.com; dkim=neutral (body hash did not verify) header.i=@gmail.com header.s=20161025 header.b=mBGdin21; 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=NONE 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 h3si7465806wmd.176.2017.07.28.04.22.55; Fri, 28 Jul 2017 04:22:56 -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=20161025 header.b=mBGdin21; 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=NONE 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 B2B49689A46; Fri, 28 Jul 2017 14:22:51 +0300 (EEST) X-Original-To: ffmpeg-devel@ffmpeg.org Delivered-To: ffmpeg-devel@ffmpeg.org Received: from mail-pf0-f195.google.com (mail-pf0-f195.google.com [209.85.192.195]) by ffbox0-bg.mplayerhq.hu (Postfix) with ESMTPS id E7260689827 for ; Fri, 28 Jul 2017 14:22:44 +0300 (EEST) Received: by mail-pf0-f195.google.com with SMTP id e3so2947344pfc.5 for ; Fri, 28 Jul 2017 04:22:46 -0700 (PDT) DKIM-Signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=20161025; h=from:to:cc:subject:date:message-id:in-reply-to:references; bh=d/fY3fKFU8Jcz3NH6YJ19aZd8PpwM4VcI3b+MV90sYQ=; b=mBGdin21PeDCroGv0njzKpuqfTtwV2TWRTeE4g9RvcivDMMXQ9JxxxuhjttInSqU3j AyuvwjSCiXiRaS8DLwC591G/RwxgdhYxw49wRVef1VN0GLy9djPsQRKZneRP5sqk8eju LwCB9DNM2SA2dCAI1B8cwrCqqw4XrPCDBdhIaz0a9dDfDH3m/skRsNSCElJVOVO8use1 aXHnDO0TPcqPmkG3kij+EtSXcl7Xl3jQsfrJ1jsILJBPjB/bLCXpqXG7Gv/IRt36FLZ4 15D/YopBTLxPe1yfO9dgLaJ4480BIZC6Krq+y4JcHDRlPc+735Pn5spV3FEgXKLt6TeG rSkQ== X-Google-DKIM-Signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=1e100.net; s=20161025; h=x-gm-message-state:from:to:cc:subject:date:message-id:in-reply-to :references; bh=d/fY3fKFU8Jcz3NH6YJ19aZd8PpwM4VcI3b+MV90sYQ=; b=rbga0I1cNjhsDap9bKv76yISdcD/Wy04JelsXVK4LXsJor/q2bg0XMCZRkz25Mxkcs bfQzi7qX0tZstUb/t/3rwIopiGZb6Zu3bCPcsbYXmY6DxBF9XYcRJRLarH66AAmKM7JU To+xLcsgNkrujWQfrJyByGebX4SLlTCsCxL0e1BVqV87Qd8pi54wCXD1g6TJvw0z17CI uJHD2jpcCmeXxEpwx/8rJQQdPipWBfeoFmGiUYpDpHruHfYNQCP3vkklw8hzIaZgUgSh eUflLZe2cK1PREqR1xG1YishIXWVPloFW6z3+oU8mcm81893/6fn+k1BFYCwDmQXMpKq /gmg== X-Gm-Message-State: AIVw113Ljowl6uNxrNek+vFKVusueKUoNdf8ClhYIp+0rIRAn0IDlHbv QQIp2cfYccxyTi5D X-Received: by 10.99.156.1 with SMTP id f1mr256211pge.327.1501239293105; Fri, 28 Jul 2017 03:54:53 -0700 (PDT) Received: from localhost.localdomain ([117.244.58.121]) by smtp.gmail.com with ESMTPSA id c63sm46007583pfk.79.2017.07.28.03.54.49 (version=TLS1_2 cipher=ECDHE-RSA-AES128-SHA bits=128/128); Fri, 28 Jul 2017 03:54:52 -0700 (PDT) From: Ashish Pratap Singh To: ffmpeg-devel@ffmpeg.org Date: Fri, 28 Jul 2017 16:23:05 +0530 Message-Id: <1501239185-4593-1-git-send-email-ashk43712@gmail.com> X-Mailer: git-send-email 2.7.4 In-Reply-To: <1499626548-3238-1-git-send-email-ashk43712@gmail.com> References: <1499626548-3238-1-git-send-email-ashk43712@gmail.com> Subject: [FFmpeg-devel] [PATCH] avfilter: add ADM filter X-BeenThere: ffmpeg-devel@ffmpeg.org X-Mailman-Version: 2.1.20 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: Ashish Singh MIME-Version: 1.0 Errors-To: ffmpeg-devel-bounces@ffmpeg.org Sender: "ffmpeg-devel" From: Ashish Singh Signed-off-by: Ashish Singh --- Changelog | 1 + doc/filters.texi | 19 ++ libavfilter/Makefile | 1 + libavfilter/adm.h | 88 +++++ libavfilter/allfilters.c | 1 + libavfilter/vf_adm.c | 819 +++++++++++++++++++++++++++++++++++++++++++++++ 6 files changed, 929 insertions(+) create mode 100644 libavfilter/adm.h create mode 100644 libavfilter/vf_adm.c diff --git a/Changelog b/Changelog index 187ae79..5a9b211 100644 --- a/Changelog +++ b/Changelog @@ -29,6 +29,7 @@ version : - limiter video filter - libvmaf video filter - Dolby E decoder and SMPTE 337M demuxer +- adm video filter version 3.3: - CrystalHD decoder moved to new decode API diff --git a/doc/filters.texi b/doc/filters.texi index 2324b96..86cf865 100644 --- a/doc/filters.texi +++ b/doc/filters.texi @@ -4604,6 +4604,25 @@ build. Below is a description of the currently available video filters. +@section adm + +Obtain the average ADM/DLM (Detail Loss Metric) between two input videos. + +This filter takes two input videos. + +Both input videos must have the same resolution and pixel format for +this filter to work correctly. Also it assumes that both inputs +have the same number of frames, which are compared one by one. + +The obtained average ADM is printed through the logging system. + +In the below example the input file @file{main.mpg} being processed is compared +with the reference file @file{ref.mpg}. + +@example +ffmpeg -i main.mpg -i ref.mpg -lavfi adm -f null - +@end example + @section alphaextract Extract the alpha component from the input as a grayscale video. This diff --git a/libavfilter/Makefile b/libavfilter/Makefile index ee16361..bca3eff 100644 --- a/libavfilter/Makefile +++ b/libavfilter/Makefile @@ -126,6 +126,7 @@ OBJS-$(CONFIG_SINE_FILTER) += asrc_sine.o OBJS-$(CONFIG_ANULLSINK_FILTER) += asink_anullsink.o # video filters +OBJS-$(CONFIG_ADM_FILTER) += vf_adm.o dualinput.o framesync.o OBJS-$(CONFIG_ALPHAEXTRACT_FILTER) += vf_extractplanes.o OBJS-$(CONFIG_ALPHAMERGE_FILTER) += vf_alphamerge.o OBJS-$(CONFIG_ASS_FILTER) += vf_subtitles.o diff --git a/libavfilter/adm.h b/libavfilter/adm.h new file mode 100644 index 0000000..8b2db8b --- /dev/null +++ b/libavfilter/adm.h @@ -0,0 +1,88 @@ +/* + * Copyright (c) 2017 Ronald S. Bultje + * Copyright (c) 2017 Ashish Pratap Singh + * + * This file is part of FFmpeg. + * + * FFmpeg is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * FFmpeg is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with FFmpeg; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA + */ + +#ifndef AVFILTER_ADM_H +#define AVFILTER_ADM_H + +#define VIEW_DIST 3.0f +#define REF_DISPLAY_HEIGHT 1080 +/** Percentage of frame to discard on all 4 sides */ +#define ADM_BORDER_FACTOR (0.1) + +typedef struct adm_dwt_band_t { + float *band_a; /** Low-pass V + low-pass H. */ + float *band_v; /** Low-pass V + high-pass H. */ + float *band_h; /** High-pass V + low-pass H. */ + float *band_d; /** High-pass V + high-pass H. */ +} adm_dwt_band_t; + +/** + * The following dwt visibility threshold parameters are taken from + * "Visibility of Wavelet Quantization Noise" + * by A. B. Watson, G. Y. Yang, J. A. Solomon and J. Villasenor + * IEEE Trans. on Image Processing, Vol. 6, No 8, Aug. 1997 + * Page 1170, formula (7) and corresponding Table IV + * Table IV has 2 entries for Cb and Cr thresholds + * Chose those corresponding to subject "sfl" since they are lower + * These thresholds were obtained and modeled for the 7-9 biorthogonal wavelet basis + */ +struct dwt_model_params { + float a; + float k; + float f0; + float g[4]; +}; + +static const float dwt2_db2_coeffs_lo[4] = { 0.482962913144690, 0.836516303737469, 0.224143868041857, -0.129409522550921 }; +static const float dwt2_db2_coeffs_hi[4] = { -0.129409522550921, -0.224143868041857, 0.836516303737469, -0.482962913144690 }; + +/** 0 -> Y, 1 -> Cb, 2 -> Cr */ +static const struct dwt_model_params dwt_7_9_YCbCr_threshold[3] = { + { .a = 0.495, .k = 0.466, .f0 = 0.401, .g = { 1.501, 1.0, 0.534, 1.0} }, + { .a = 1.633, .k = 0.353, .f0 = 0.209, .g = { 1.520, 1.0, 0.502, 1.0} }, + { .a = 0.944, .k = 0.521, .f0 = 0.404, .g = { 1.868, 1.0, 0.516, 1.0} } +}; + +/** + * The following dwt basis function amplitudes, A(lambda,theta), are taken from + * "Visibility of Wavelet Quantization Noise" + * by A. B. Watson, G. Y. Yang, J. A. Solomon and J. Villasenor + * IEEE Trans. on Image Processing, Vol. 6, No 8, Aug. 1997 + * Page 1172, Table V + * The table has been transposed, i.e. it can be used directly to obtain A[lambda][theta] + * These amplitudes were calculated for the 7-9 biorthogonal wavelet basis + */ +static const float dwt_7_9_basis_function_amplitudes[6][4] = { + { 0.62171, 0.67234, 0.72709, 0.67234 }, + { 0.34537, 0.41317, 0.49428, 0.41317 }, + { 0.18004, 0.22727, 0.28688, 0.22727 }, + { 0.091401, 0.11792, 0.15214, 0.11792 }, + { 0.045943, 0.059758, 0.077727, 0.059758 }, + { 0.023013, 0.030018, 0.039156, 0.030018 } +}; + +/** function to compute adm score */ +int compute_adm2(const float *ref, const float *main, int w, int h, + int ref_stride, int main_stride, double *score, + double *score_num, double *score_den, double *scores, + float *data_buf, float *temp_lo, float* temp_hi); + +#endif /* AVFILTER_ADM_H */ diff --git a/libavfilter/allfilters.c b/libavfilter/allfilters.c index b1c2d11..bbd7db8 100644 --- a/libavfilter/allfilters.c +++ b/libavfilter/allfilters.c @@ -138,6 +138,7 @@ static void register_all(void) REGISTER_FILTER(ANULLSINK, anullsink, asink); + REGISTER_FILTER(ADM, adm, vf); REGISTER_FILTER(ALPHAEXTRACT, alphaextract, vf); REGISTER_FILTER(ALPHAMERGE, alphamerge, vf); REGISTER_FILTER(ASS, ass, vf); diff --git a/libavfilter/vf_adm.c b/libavfilter/vf_adm.c new file mode 100644 index 0000000..0c3effc --- /dev/null +++ b/libavfilter/vf_adm.c @@ -0,0 +1,819 @@ +/* + * Copyright (c) 2017 Ronald S. Bultje + * Copyright (c) 2017 Ashish Pratap Singh + * + * This file is part of FFmpeg. + * + * FFmpeg is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * FFmpeg is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with FFmpeg; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA + */ + +/** + * @file + * Calculate the ADM between two input videos. + */ + +#include "libavutil/avstring.h" +#include "libavutil/opt.h" +#include "libavutil/pixdesc.h" +#include "avfilter.h" +#include "dualinput.h" +#include "drawutils.h" +#include "formats.h" +#include "internal.h" +#include "adm.h" +#include "video.h" +#include + +typedef struct ADMContext { + const AVClass *class; + FFDualInputContext dinput; + const AVPixFmtDescriptor *desc; + int width; + int height; + float *ref_data; + float *main_data; + float *data_buf; + float *temp_lo; + float *temp_hi; + double adm_sum; + uint64_t nb_frames; +} ADMContext; + +static const AVOption adm_options[] = { + { NULL } +}; + +AVFILTER_DEFINE_CLASS(adm); + +#define MAX_ALIGN 32 +#define ALIGN_CEIL(x) ((x) + ((x) % MAX_ALIGN ? MAX_ALIGN - (x) % MAX_ALIGN : 0)) +#define OPT_RANGE_PIXEL_OFFSET (-128) +#define FLAGS AV_OPT_FLAG_FILTERING_PARAM|AV_OPT_FLAG_VIDEO_PARAM + +static float rcp(float x) +{ + float xi = _mm_cvtss_f32(_mm_rcp_ss(_mm_load_ss(&x))); + return xi + xi * (1.0f - x * xi); +} + +#define DIVS(n, d) ((n) * rcp(d)) + +/** + * lambda = 0 (finest scale), 1, 2, 3 (coarsest scale); + * theta = 0 (ll), 1 (lh - vertical), 2 (hh - diagonal), 3(hl - horizontal). + */ +av_always_inline static float dwt_quant_step(const struct dwt_model_params *params, + int lambda, int theta) +{ + /** Formula (1), page 1165 - display visual resolution (DVR), + * in pixels/degree of visual angle. This should be 56.55 + */ + float r = VIEW_DIST * REF_DISPLAY_HEIGHT * M_PI / 180.0; + + /** Formula (9), page 1171 */ + float temp = log10(pow(2.0, lambda+1) * params->f0 * params->g[theta] / r); + float Q = 2.0 * params->a * pow(10.0, params->k * temp * temp) / + dwt_7_9_basis_function_amplitudes[lambda][theta]; + + return Q; +} + +static float get_cube(float val) +{ + return val * val * val; +} + +static float adm_sum_cube(const float *x, int w, int h, int stride, + double border_factor) +{ + int px_stride = stride / sizeof(float); + int left = w * border_factor - 0.5; + int top = h * border_factor - 0.5; + int right = w - left; + int bottom = h - top; + + int i, j; + + float val; + float sum = 0; + + for (i = top; i < bottom; i++) { + for (j = left; j < right; j++) { + val = fabsf(x[i * px_stride + j]); + sum += get_cube(val); + } + } + + return powf(sum, 1.0f / 3.0f) + powf((bottom - top) * (right - left) / + 32.0f, 1.0f / 3.0f); +} + +static void adm_decouple(const adm_dwt_band_t *ref, const adm_dwt_band_t *main, + const adm_dwt_band_t *r, const adm_dwt_band_t *a, + int w, int h, int ref_stride, int main_stride, + int r_stride, int a_stride) +{ + const float cos_1deg_sq = cos(1.0 * M_PI / 180.0) * cos(1.0 * M_PI / 180.0); + const float eps = 1e-30; + + int ref_px_stride = ref_stride / sizeof(float); + int main_px_stride = main_stride / sizeof(float); + int r_px_stride = r_stride / sizeof(float); + int a_px_stride = a_stride / sizeof(float); + + float oh, ov, od, th, tv, td; + float kh, kv, kd, tmph, tmpv, tmpd; + float ot_dp, o_mag_sq, t_mag_sq; + int angle_flag; + int i, j; + + for (i = 0; i < h; i++) { + for (j = 0; j < w; j++) { + oh = ref->band_h[i * ref_px_stride + j]; + ov = ref->band_v[i * ref_px_stride + j]; + od = ref->band_d[i * ref_px_stride + j]; + th = main->band_h[i * main_px_stride + j]; + tv = main->band_v[i * main_px_stride + j]; + td = main->band_d[i * main_px_stride + j]; + + kh = DIVS(th, oh + eps); + kv = DIVS(tv, ov + eps); + kd = DIVS(td, od + eps); + + kh = kh < 0.0f ? 0.0f : (kh > 1.0f ? 1.0f : kh); + + + kv = kv < 0.0f ? 0.0f : (kv > 1.0f ? 1.0f : kv); + kd = kd < 0.0f ? 0.0f : (kd > 1.0f ? 1.0f : kd); + + tmph = kh * oh; + tmpv = kv * ov; + tmpd = kd * od; + /** Determine if angle between (oh,ov) and (th,tv) is less than 1 degree. + * Given that u is the angle (oh,ov) and v is the angle (th,tv), this can + * be done by testing the inequvality. + * + * { (u.v.) >= 0 } AND { (u.v)^2 >= cos(1deg)^2 * ||u||^2 * ||v||^2 } + * + * Proof: + * + * cos(theta) = (u.v) / (||u|| * ||v||) + * + * IF u.v >= 0 THEN + * cos(theta)^2 = (u.v)^2 / (||u||^2 * ||v||^2) + * (u.v)^2 = cos(theta)^2 * ||u||^2 * ||v||^2 + * + * IF |theta| < 1deg THEN + * (u.v)^2 >= cos(1deg)^2 * ||u||^2 * ||v||^2 + * END + * ELSE + * |theta| > 90deg + * END + */ + ot_dp = oh * th + ov * tv; + o_mag_sq = oh * oh + ov * ov; + t_mag_sq = th * th + tv * tv; + + angle_flag = (ot_dp >= 0.0f) && (ot_dp * ot_dp >= cos_1deg_sq * + o_mag_sq * t_mag_sq); + + if (angle_flag) { + tmph = th; + tmpv = tv; + tmpd = td; + } + + r->band_h[i * r_px_stride + j] = tmph; + r->band_v[i * r_px_stride + j] = tmpv; + r->band_d[i * r_px_stride + j] = tmpd; + + a->band_h[i * a_px_stride + j] = th - tmph; + a->band_v[i * a_px_stride + j] = tv - tmpv; + a->band_d[i * a_px_stride + j] = td - tmpd; + } + } +} + +static void adm_csf(const adm_dwt_band_t *src, const adm_dwt_band_t *dst, + int orig_h, int scale, int w, int h, int src_stride, + int dst_stride) +{ + const float *src_angles[3] = { src->band_h, src->band_v, src->band_d }; + float *dst_angles[3] = { dst->band_h, dst->band_v, dst->band_d }; + + const float *src_ptr; + float *dst_ptr; + + int src_px_stride = src_stride / sizeof(float); + int dst_px_stride = dst_stride / sizeof(float); + + float factor1 = dwt_quant_step(&dwt_7_9_YCbCr_threshold[0], scale, 1); + float factor2 = dwt_quant_step(&dwt_7_9_YCbCr_threshold[0], scale, 2); + float rfactor[3] = {1.0f / factor1, 1.0f / factor1, 1.0f / factor2}; + + int i, j, theta; + + for (theta = 0; theta < 3; theta++) { + src_ptr = src_angles[theta]; + dst_ptr = dst_angles[theta]; + + for (i = 0; i < h; i++) { + for (j = 0; j < w; j++) { + dst_ptr[i * dst_px_stride + j] = rfactor[theta] * + src_ptr[i * src_px_stride + j]; + } + } + } +} + +static void adm_cm_thresh(const adm_dwt_band_t *src, float *dst, int w, int h, + int src_stride, int dst_stride) +{ + const float *angles[3] = { src->band_h, src->band_v, src->band_d }; + const float *src_ptr; + + int src_px_stride = src_stride / sizeof(float); + int dst_px_stride = dst_stride / sizeof(float); + + float filt_coeff, img_coeff; + + int theta, i, j, filt_i, filt_j, src_i, src_j; + + for (i = 0; i < h; i++) { + + for (j = 0; j < w; j++) { + dst[i * dst_px_stride + j] = 0; + } + + for (theta = 0; theta < 3; ++theta) { + src_ptr = angles[theta]; + + for (j = 0; j < w; j++) { + float sum = 0; + + for (filt_i = 0; filt_i < 3; filt_i++) { + for (filt_j = 0; filt_j < 3; filt_j++) { + filt_coeff = (filt_i == 1 && filt_j == 1) ? 1.0f / 15.0f : 1.0f / + 30.0f; + + src_i = i - 1 + filt_i; + src_j = j - 1 + filt_j; + + src_i = FFABS(src_i); + if (src_i >= h) { + src_i = 2 * h - src_i - 1; + } + src_j = FFABS(src_j); + if (src_j >= w) { + src_j = 2 * w - src_j - 1; + } + img_coeff = fabsf(src_ptr[src_i * src_px_stride + src_j]); + + sum += filt_coeff * img_coeff; + } + } + + dst[i * dst_px_stride + j] += sum; + } + } + } +} + +static void adm_cm(const adm_dwt_band_t *src, const adm_dwt_band_t *dst, + const float *thresh, int w, int h, int src_stride, + int dst_stride, int thresh_stride) +{ + int src_px_stride = src_stride / sizeof(float); + int dst_px_stride = dst_stride / sizeof(float); + int thresh_px_stride = thresh_stride / sizeof(float); + + float xh, xv, xd, thr; + + int i, j; + + for (i = 0; i < h; i++) { + for (j = 0; j < w; j++) { + xh = src->band_h[i * src_px_stride + j]; + xv = src->band_v[i * src_px_stride + j]; + xd = src->band_d[i * src_px_stride + j]; + thr = thresh[i * thresh_px_stride + j]; + + xh = fabsf(xh) - thr; + xv = fabsf(xv) - thr; + xd = fabsf(xd) - thr; + + xh = xh < 0.0f ? 0.0f : xh; + xv = xv < 0.0f ? 0.0f : xv; + xd = xd < 0.0f ? 0.0f : xd; + + dst->band_h[i * dst_px_stride + j] = xh; + dst->band_v[i * dst_px_stride + j] = xv; + dst->band_d[i * dst_px_stride + j] = xd; + } + } +} + +static void adm_dwt2(const float *src, const adm_dwt_band_t *dst, int w, int h, + int src_stride, int dst_stride, float *temp_lo, + float* temp_hi) +{ + const float *filter_lo = dwt2_db2_coeffs_lo; + const float *filter_hi = dwt2_db2_coeffs_hi; + int filt_w = sizeof(dwt2_db2_coeffs_lo) / sizeof(float); + + int src_px_stride = src_stride / sizeof(float); + int dst_px_stride = dst_stride / sizeof(float); + + float filt_coeff_lo, filt_coeff_hi, img_coeff; + + int i, j, filt_i, filt_j, src_i, src_j; + + for (i = 0; i < (h + 1) / 2; i++) { + /** Vertical pass. */ + for (j = 0; j < w; j++) { + float sum_lo = 0; + float sum_hi = 0; + + for (filt_i = 0; filt_i < filt_w; filt_i++) { + filt_coeff_lo = filter_lo[filt_i]; + filt_coeff_hi = filter_hi[filt_i]; + + src_i = 2 * i - 1 + filt_i; + + src_i = FFABS(src_i); + if (src_i >= h) { + src_i = 2 * h - src_i - 1; + } + + img_coeff = src[src_i * src_px_stride + j]; + + sum_lo += filt_coeff_lo * img_coeff; + sum_hi += filt_coeff_hi * img_coeff; + } + + temp_lo[j] = sum_lo; + temp_hi[j] = sum_hi; + } + + /** Horizontal pass (lo). */ + for (j = 0; j < (w + 1) / 2; j++) { + float sum_lo = 0; + float sum_hi = 0; + + for (filt_j = 0; filt_j < filt_w; filt_j++) { + filt_coeff_lo = filter_lo[filt_j]; + filt_coeff_hi = filter_hi[filt_j]; + + src_j = 2 * j - 1 + filt_j; + + src_j = FFABS(src_j); + if (src_j >= w) { + src_j = 2 * w - src_j - 1; + } + + img_coeff = temp_lo[src_j]; + + sum_lo += filt_coeff_lo * img_coeff; + sum_hi += filt_coeff_hi * img_coeff; + } + + dst->band_a[i * dst_px_stride + j] = sum_lo; + dst->band_v[i * dst_px_stride + j] = sum_hi; + } + + /** Horizontal pass (hi). */ + for (j = 0; j < (w + 1) / 2; j++) { + float sum_lo = 0; + float sum_hi = 0; + + for (filt_j = 0; filt_j < filt_w; filt_j++) { + filt_coeff_lo = filter_lo[filt_j]; + filt_coeff_hi = filter_hi[filt_j]; + + src_j = 2 * j - 1 + filt_j; + + src_j = FFABS(src_j); + if (src_j >= w) { + src_j = 2 * w - src_j - 1; + } + + img_coeff = temp_hi[src_j]; + + sum_lo += filt_coeff_lo * img_coeff; + sum_hi += filt_coeff_hi * img_coeff; + } + + dst->band_h[i * dst_px_stride + j] = sum_lo; + dst->band_d[i * dst_px_stride + j] = sum_hi; + } + } + +} + +static void adm_buffer_copy(const void *src, void *dst, int linewidth, int h, + int src_stride, int dst_stride) +{ + const char *src_p = src; + char *dst_p = dst; + int i; + + for (i = 0; i < h; i++) { + memcpy(dst_p, src_p, linewidth); + src_p += src_stride; + dst_p += dst_stride; + } +} + +static char *init_dwt_band(adm_dwt_band_t *band, char *data_top, size_t buf_sz) +{ + band->band_a = (float *) data_top; + data_top += buf_sz; + band->band_h = (float *) data_top; + data_top += buf_sz; + band->band_v = (float *) data_top; + data_top += buf_sz; + band->band_d = (float *) data_top; + data_top += buf_sz; + return data_top; +} + +int compute_adm2(const float *ref, const float *main, int w, int h, + int ref_stride, int main_stride, double *score, + double *score_num, double *score_den, double *scores, + float *data_buf, float *temp_lo, float* temp_hi) +{ + double numden_limit = 1e-2 * (w * h) / (1920.0 * 1080.0); + + char *data_top; + + float *ref_scale; + float *main_scale; + + adm_dwt_band_t ref_dwt2; + adm_dwt_band_t main_dwt2; + + adm_dwt_band_t decouple_r; + adm_dwt_band_t decouple_a; + + adm_dwt_band_t csf_o; + adm_dwt_band_t csf_r; + adm_dwt_band_t csf_a; + + float *mta; + + adm_dwt_band_t cm_r; + + const float *curr_ref_scale = (float *) ref; + const float *curr_main_scale = (float *) main; + int curr_ref_stride = ref_stride; + int curr_main_stride = main_stride; + + int orig_h = h; + + int buf_stride = ALIGN_CEIL(((w + 1) / 2) * sizeof(float)); + size_t buf_sz = (size_t)buf_stride * ((h + 1) / 2); + + double num = 0; + double den = 0; + + int scale; + int ret = 1; + + data_top = (char *) (data_buf); + + ref_scale = (float *) data_top; + data_top += buf_sz; + main_scale = (float *) data_top; + data_top += buf_sz; + + data_top = init_dwt_band(&ref_dwt2, data_top, buf_sz); + data_top = init_dwt_band(&main_dwt2, data_top, buf_sz); + data_top = init_dwt_band(&decouple_r, data_top, buf_sz); + data_top = init_dwt_band(&decouple_a, data_top, buf_sz); + data_top = init_dwt_band(&csf_o, data_top, buf_sz); + data_top = init_dwt_band(&csf_r, data_top, buf_sz); + data_top = init_dwt_band(&csf_a, data_top, buf_sz); + + mta = (float *) data_top; + data_top += buf_sz; + + data_top = init_dwt_band(&cm_r, data_top, buf_sz); + + for (scale = 0; scale < 4; scale++) { + float num_scale = 0.0; + float den_scale = 0.0; + + adm_dwt2(curr_ref_scale, &ref_dwt2, w, h, curr_ref_stride, buf_stride, + temp_lo, temp_hi); + adm_dwt2(curr_main_scale, &main_dwt2, w, h, curr_main_stride, buf_stride, + temp_lo, temp_hi); + + w = (w + 1) / 2; + h = (h + 1) / 2; + + adm_decouple(&ref_dwt2, &main_dwt2, &decouple_r, &decouple_a, w, h, + buf_stride, buf_stride, buf_stride, buf_stride); + + adm_csf(&ref_dwt2, &csf_o, orig_h, scale, w, h, buf_stride, buf_stride); + adm_csf(&decouple_r, &csf_r, orig_h, scale, w, h, buf_stride, buf_stride); + adm_csf(&decouple_a, &csf_a, orig_h, scale, w, h, buf_stride, buf_stride); + + adm_cm_thresh(&csf_a, mta, w, h, buf_stride, buf_stride); + adm_cm(&csf_r, &cm_r, mta, w, h, buf_stride, buf_stride, buf_stride); + + num_scale += adm_sum_cube(cm_r.band_h, w, h, buf_stride, ADM_BORDER_FACTOR); + num_scale += adm_sum_cube(cm_r.band_v, w, h, buf_stride, ADM_BORDER_FACTOR); + num_scale += adm_sum_cube(cm_r.band_d, w, h, buf_stride, ADM_BORDER_FACTOR); + + den_scale += adm_sum_cube(csf_o.band_h, w, h, buf_stride, ADM_BORDER_FACTOR); + den_scale += adm_sum_cube(csf_o.band_v, w, h, buf_stride, ADM_BORDER_FACTOR); + den_scale += adm_sum_cube(csf_o.band_d, w, h, buf_stride, ADM_BORDER_FACTOR); + + num += num_scale; + den += den_scale; + + adm_buffer_copy(ref_dwt2.band_a, ref_scale, w * sizeof(float), h, + buf_stride, buf_stride); + adm_buffer_copy(main_dwt2.band_a, main_scale, w * sizeof(float), h, + buf_stride, buf_stride); + + curr_ref_scale = ref_scale; + curr_main_scale = main_scale; + curr_ref_stride = buf_stride; + curr_main_stride = buf_stride; + + scores[2*scale+0] = num_scale; + scores[2*scale+1] = den_scale; + } + + num = num < numden_limit ? 0 : num; + den = den < numden_limit ? 0 : den; + + if (den == 0.0) { + *score = 1.0f; + } else { + *score = num / den; + } + *score_num = num; + *score_den = den; + + ret = 0; + + return ret; +} + +#define offset_fn(type, bits) \ + static void offset_##bits##bit(ADMContext *s, const AVFrame *ref, AVFrame *main, int stride) \ +{ \ + int w = s->width; \ + int h = s->height; \ + int i,j; \ + \ + int ref_stride = ref->linesize[0]; \ + int main_stride = main->linesize[0]; \ + \ + const type *ref_ptr = (const type *) ref->data[0]; \ + const type *main_ptr = (const type *) main->data[0]; \ + \ + float *ref_ptr_data = s->ref_data; \ + float *main_ptr_data = s->main_data; \ + \ + for(i = 0; i < h; i++) { \ + for(j = 0; j < w; j++) { \ + ref_ptr_data[j] = (float) ref_ptr[j] + OPT_RANGE_PIXEL_OFFSET; \ + main_ptr_data[j] = (float) main_ptr[j] + OPT_RANGE_PIXEL_OFFSET; \ + } \ + ref_ptr += ref_stride / sizeof(type); \ + ref_ptr_data += stride / sizeof(float); \ + main_ptr += main_stride / sizeof(type); \ + main_ptr_data += stride / sizeof(float); \ + } \ +} + +offset_fn(uint8_t, 8); +offset_fn(uint16_t, 10); + +static void set_meta(AVDictionary **metadata, const char *key, float d) +{ + char value[128]; + snprintf(value, sizeof(value), "%0.2f", d); + av_dict_set(metadata, key, value, 0); +} + +static AVFrame *do_adm(AVFilterContext *ctx, AVFrame *main, const AVFrame *ref) +{ + ADMContext *s = ctx->priv; + AVDictionary **metadata = &main->metadata; + + double score = 0.0; + double score_num = 0; + double score_den = 0; + double scores[2*4]; + + int w = s->width; + int h = s->height; + + int stride; + + stride = ALIGN_CEIL(w * sizeof(float)); + + /** Offset ref and main pixels by OPT_RANGE_PIXEL_OFFSET */ + if (s->desc->comp[0].depth <= 8) { + offset_8bit(s, ref, main, stride); + } else { + offset_10bit(s, ref, main, stride); + } + + compute_adm2(s->ref_data, s->main_data, w, h, stride, stride, &score, + &score_num, &score_den, scores, s->data_buf, s->temp_lo, + s->temp_hi); + + set_meta(metadata, "lavfi.adm.score", score); + + s->nb_frames++; + + s->adm_sum += score; + + return main; +} + +static av_cold int init(AVFilterContext *ctx) +{ + ADMContext *s = ctx->priv; + + s->dinput.process = do_adm; + + return 0; +} + +static int query_formats(AVFilterContext *ctx) +{ + static const enum AVPixelFormat pix_fmts[] = { + AV_PIX_FMT_YUV444P, AV_PIX_FMT_YUV422P, AV_PIX_FMT_YUV420P, + AV_PIX_FMT_YUV444P10LE, AV_PIX_FMT_YUV422P10LE, AV_PIX_FMT_YUV420P10LE, + AV_PIX_FMT_NONE + }; + + AVFilterFormats *fmts_list = ff_make_format_list(pix_fmts); + if (!fmts_list) + return AVERROR(ENOMEM); + return ff_set_common_formats(ctx, fmts_list); +} + +static int config_input_ref(AVFilterLink *inlink) +{ + AVFilterContext *ctx = inlink->dst; + ADMContext *s = ctx->priv; + int buf_stride; + size_t buf_sz; + int stride; + size_t data_sz; + + if (ctx->inputs[0]->w != ctx->inputs[1]->w || + ctx->inputs[0]->h != ctx->inputs[1]->h) { + av_log(ctx, AV_LOG_ERROR, "Width and height of input videos must be same.\n"); + return AVERROR(EINVAL); + } + if (ctx->inputs[0]->format != ctx->inputs[1]->format) { + av_log(ctx, AV_LOG_ERROR, "Inputs must be of same pixel format.\n"); + return AVERROR(EINVAL); + } + + s->desc = av_pix_fmt_desc_get(inlink->format); + s->width = ctx->inputs[0]->w; + s->height = ctx->inputs[0]->h; + + stride = ALIGN_CEIL(s->width * sizeof(float)); + data_sz = (size_t)stride * s->height; + + if (!(s->ref_data = av_malloc(data_sz))) { + av_log(ctx, AV_LOG_ERROR, "ref data allocation failed.\n"); + return AVERROR(ENOMEM); + } + + if (!(s->main_data = av_malloc(data_sz))) { + av_log(ctx, AV_LOG_ERROR, "main data allocation failed.\n"); + return AVERROR(ENOMEM); + } + + buf_stride = ALIGN_CEIL(((s->width + 1) / 2) * sizeof(float)); + buf_sz = (size_t)buf_stride * ((s->height + 1) / 2); + + if (SIZE_MAX / buf_sz < 35) { + av_log(ctx, AV_LOG_ERROR, "error: SIZE_MAX / buf_sz_one < 35, buf_sz_one = %lu.\n", buf_sz); + return AVERROR(EINVAL); + } + + if (!(s->data_buf = av_malloc(buf_sz * 35))) { + av_log(ctx, AV_LOG_ERROR, "data_buf allocation failed.\n"); + return AVERROR(ENOMEM); + } + + if (!(s->temp_lo = av_malloc(stride))) { + av_log(ctx, AV_LOG_ERROR, "temp lo allocation failed.\n"); + return AVERROR(ENOMEM); + } + + if (!(s->temp_hi = av_malloc(stride))) { + av_log(ctx, AV_LOG_ERROR, "temp hi allocation failed.\n"); + return AVERROR(ENOMEM); + } + + return 0; +} + + +static int config_output(AVFilterLink *outlink) +{ + AVFilterContext *ctx = outlink->src; + ADMContext *s = ctx->priv; + AVFilterLink *mainlink = ctx->inputs[0]; + int ret; + + outlink->w = mainlink->w; + outlink->h = mainlink->h; + outlink->time_base = mainlink->time_base; + outlink->sample_aspect_ratio = mainlink->sample_aspect_ratio; + outlink->frame_rate = mainlink->frame_rate; + if ((ret = ff_dualinput_init(ctx, &s->dinput)) < 0) + return ret; + + return 0; +} + +static int filter_frame(AVFilterLink *inlink, AVFrame *inpicref) +{ + ADMContext *s = inlink->dst->priv; + return ff_dualinput_filter_frame(&s->dinput, inlink, inpicref); +} + +static int request_frame(AVFilterLink *outlink) +{ + ADMContext *s = outlink->src->priv; + return ff_dualinput_request_frame(&s->dinput, outlink); +} + +static av_cold void uninit(AVFilterContext *ctx) +{ + ADMContext *s = ctx->priv; + + if (s->nb_frames > 0) { + av_log(ctx, AV_LOG_INFO, "ADM AVG: %.3f\n", s->adm_sum / s->nb_frames); + } + + av_free(s->ref_data); + av_free(s->main_data); + av_free(s->data_buf); + av_free(s->temp_lo); + av_free(s->temp_hi); + + ff_dualinput_uninit(&s->dinput); +} + +static const AVFilterPad adm_inputs[] = { + { + .name = "main", + .type = AVMEDIA_TYPE_VIDEO, + .filter_frame = filter_frame, + },{ + .name = "reference", + .type = AVMEDIA_TYPE_VIDEO, + .filter_frame = filter_frame, + .config_props = config_input_ref, + }, + { NULL } +}; + +static const AVFilterPad adm_outputs[] = { + { + .name = "default", + .type = AVMEDIA_TYPE_VIDEO, + .config_props = config_output, + .request_frame = request_frame, + }, + { NULL } +}; + +AVFilter ff_vf_adm = { + .name = "adm", + .description = NULL_IF_CONFIG_SMALL("Calculate the ADM between two video streams."), + .init = init, + .uninit = uninit, + .query_formats = query_formats, + .priv_size = sizeof(ADMContext), + .priv_class = &adm_class, + .inputs = adm_inputs, + .outputs = adm_outputs, +};