From patchwork Sun Sep 18 11:36:29 2016 Content-Type: text/plain; charset="utf-8" MIME-Version: 1.0 Content-Transfer-Encoding: 8bit X-Patchwork-Submitter: =?utf-8?b?Q2zDqW1lbnQgQsWTc2No?= X-Patchwork-Id: 613 Delivered-To: ffmpegpatchwork@gmail.com Received: by 10.103.140.134 with SMTP id o128csp531910vsd; Sun, 18 Sep 2016 04:36:05 -0700 (PDT) X-Received: by 10.28.25.71 with SMTP id 68mr4928575wmz.91.1474198564949; Sun, 18 Sep 2016 04:36:04 -0700 (PDT) Return-Path: Received: from ffbox0-bg.mplayerhq.hu (ffbox0-bg.ffmpeg.org. [79.124.17.100]) by mx.google.com with ESMTP id 16si13074858wmh.7.2016.09.18.04.36.01; Sun, 18 Sep 2016 04:36:04 -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; 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 Received: from [127.0.1.1] (localhost [127.0.0.1]) by ffbox0-bg.mplayerhq.hu (Postfix) with ESMTP id 0F83A68A071; Sun, 18 Sep 2016 14:35:45 +0300 (EEST) X-Original-To: ffmpeg-devel@ffmpeg.org Delivered-To: ffmpeg-devel@ffmpeg.org Received: from kimiko.pkh.me (pkh.me [5.196.73.185]) by ffbox0-bg.mplayerhq.hu (Postfix) with ESMTPS id A133F689F4D for ; Sun, 18 Sep 2016 14:35:38 +0300 (EEST) Received: from localhost (kimiko.pkh.me [local]) by kimiko.pkh.me (OpenSMTPD) with ESMTPA id c12a1562; Sun, 18 Sep 2016 11:36:30 +0000 (UTC) From: =?UTF-8?q?Cl=C3=A9ment=20B=C5=93sch?= To: ffmpeg-devel@ffmpeg.org Date: Sun, 18 Sep 2016 13:36:29 +0200 Message-Id: <20160918113629.16877-1-u@pkh.me> X-Mailer: git-send-email 2.9.3 MIME-Version: 1.0 Subject: [FFmpeg-devel] [PATCH] avfilter/WIP: add nlmeans 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: =?UTF-8?q?Cl=C3=A9ment=20B=C5=93sch?= Errors-To: ffmpeg-devel-bounces@ffmpeg.org Sender: "ffmpeg-devel" Fixes Ticket #4910 TODO (before 1st pushed version): - add chroma paramters (for patch size and research window) - doc - lavfi minor bump, Changelog - GBRP TODO++ (after 1st version): - SIMD for compute_safe_ssd_integral_image - SIMD for final weighted averaging - smart parameters for block and research size according to sigma like suggested by ipol? - temporal support --- libavfilter/Makefile | 3 +- libavfilter/allfilters.c | 1 + libavfilter/tests/integral.c | 89 ++++++++ libavfilter/vf_nlmeans.c | 527 +++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 619 insertions(+), 1 deletion(-) create mode 100644 libavfilter/tests/integral.c create mode 100644 libavfilter/vf_nlmeans.c diff --git a/libavfilter/Makefile b/libavfilter/Makefile index 5cd10fa..57a38d3 100644 --- a/libavfilter/Makefile +++ b/libavfilter/Makefile @@ -220,6 +220,7 @@ OBJS-$(CONFIG_METADATA_FILTER) += f_metadata.o OBJS-$(CONFIG_MINTERPOLATE_FILTER) += vf_minterpolate.o motion_estimation.o OBJS-$(CONFIG_MPDECIMATE_FILTER) += vf_mpdecimate.o OBJS-$(CONFIG_NEGATE_FILTER) += vf_lut.o +OBJS-$(CONFIG_NLMEANS_FILTER) += vf_nlmeans.o OBJS-$(CONFIG_NNEDI_FILTER) += vf_nnedi.o OBJS-$(CONFIG_NOFORMAT_FILTER) += vf_format.o OBJS-$(CONFIG_NOISE_FILTER) += vf_noise.o @@ -354,7 +355,7 @@ SKIPHEADERS-$(CONFIG_OPENCL) += opencl_internal.h deshake_opencl OBJS-$(CONFIG_SHARED) += log2_tab.o TOOLS = graph2dot -TESTPROGS = drawutils filtfmts formats +TESTPROGS = drawutils filtfmts formats integral TOOLS-$(CONFIG_LIBZMQ) += zmqsend diff --git a/libavfilter/allfilters.c b/libavfilter/allfilters.c index 47d95f5..8f542fd 100644 --- a/libavfilter/allfilters.c +++ b/libavfilter/allfilters.c @@ -237,6 +237,7 @@ void avfilter_register_all(void) REGISTER_FILTER(MINTERPOLATE, minterpolate, vf); REGISTER_FILTER(MPDECIMATE, mpdecimate, vf); REGISTER_FILTER(NEGATE, negate, vf); + REGISTER_FILTER(NLMEANS, nlmeans, vf); REGISTER_FILTER(NNEDI, nnedi, vf); REGISTER_FILTER(NOFORMAT, noformat, vf); REGISTER_FILTER(NOISE, noise, vf); diff --git a/libavfilter/tests/integral.c b/libavfilter/tests/integral.c new file mode 100644 index 0000000..7690254 --- /dev/null +++ b/libavfilter/tests/integral.c @@ -0,0 +1,89 @@ +/* + * 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 + */ + +#include "libavfilter/vf_nlmeans.c" + +int main(void) +{ + int ret = 0, xoff, yoff; + + // arbitrary test source of size 6x4 and linesize=8 + const int w = 6, h = 5, lz = 8; + static const uint8_t src[] = { + 0xb0, 0x71, 0xfb, 0xd8, 0x01, 0xd9, /***/ 0x01, 0x02, + 0x51, 0x8e, 0x41, 0x0f, 0x84, 0x58, /***/ 0x03, 0x04, + 0xc7, 0x8d, 0x07, 0x70, 0x5c, 0x47, /***/ 0x05, 0x06, + 0x09, 0x4e, 0xfc, 0x74, 0x8f, 0x9a, /***/ 0x07, 0x08, + 0x60, 0x8e, 0x20, 0xaa, 0x95, 0x7d, /***/ 0x09, 0x0a, + }; + + const int r = 3; + const int ii_w = w+r*2, ii_h = h+r*2; + + // align to 4 the linesize, "+1" is for the space of the left 0-column + const int ii_lz_32 = ((ii_w + 1) + 3) & ~3; + + // "+1" is for the space of the top 0-line + uint32_t *ii = calloc(ii_h + 1, ii_lz_32 * sizeof(*ii)); + uint32_t *ii2 = calloc(ii_h + 1, ii_lz_32 * sizeof(*ii2)); + + uint32_t *ii_start = ii + ii_lz_32 + 1; // skip top 0-line and left 0-column + uint32_t *ii_start2 = ii2 + ii_lz_32 + 1; // skip top 0-line and left 0-column + + for (yoff = -r; yoff <= r; yoff++) { + for (xoff = -r; xoff <= r; xoff++) { + int x, y; + + printf("xoff=%d yoff=%d\n", xoff, yoff); + + compute_ssd_integral_image(ii_start, ii_lz_32, + src, lz, xoff, yoff, r, w, h); + + for (y = 0; y < ii_h; y++) { + for (x = 0; x < ii_w; x++) + printf(" %7x", ii_start[y*ii_lz_32 + x]); + printf("\n"); + } + printf("---------------\n"); + + compute_unsafe_ssd_integral_image(ii_start2, ii_lz_32, + 0, 0, + src, lz, + xoff, yoff, r, w, h, + ii_w, ii_h); + + for (y = 0; y < ii_h; y++) { + for (x = 0; x < ii_w; x++) + printf(" %7x", ii_start2[y*ii_lz_32 + x]); + printf("\n"); + } + printf("---------------\n"); + + if (memcmp(ii, ii2, (ii_h+1) * ii_lz_32 * sizeof(*ii))) { + printf("Integral mismatch\n"); + ret = 1; + goto end; + } + } + } + +end: + free(ii); + free(ii2); + return ret; +} diff --git a/libavfilter/vf_nlmeans.c b/libavfilter/vf_nlmeans.c new file mode 100644 index 0000000..f541a80 --- /dev/null +++ b/libavfilter/vf_nlmeans.c @@ -0,0 +1,527 @@ +/* + * Copyright (c) 2016 Clément Bœsch + * + * 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 + */ + +#include "libavutil/avassert.h" +#include "libavutil/opt.h" +#include "libavutil/pixdesc.h" +#include "avfilter.h" +#include "formats.h" +#include "internal.h" +#include "video.h" + +struct weighted_avg { + double total_weight; + double sum; +}; + +#define WEIGHT_LUT_NBITS 9 +#define WEIGHT_LUT_SIZE (1<= 0); av_assert1(startx_safe - s1x < w); + av_assert1(starty_safe - s1y >= 0); av_assert1(starty_safe - s1y < h); + av_assert1(startx_safe - s2x >= 0); av_assert1(startx_safe - s2x < w); + av_assert1(starty_safe - s2y >= 0); av_assert1(starty_safe - s2y < h); + compute_safe_ssd_integral_image_c(ii + starty_safe*ii_linesize_32 + startx_safe, ii_linesize_32, + src + (starty_safe - s1y) * linesize + (startx_safe - s1x), linesize, + src + (starty_safe - s2y) * linesize + (startx_safe - s2x), linesize, + endx_safe - startx_safe, endy_safe - starty_safe); + + // right part of the integral + compute_unsafe_ssd_integral_image(ii, ii_linesize_32, + endx_safe, starty_safe, + src, linesize, + offx, offy, e, w, h, + ii_w - endx_safe, endy_safe - starty_safe); + + // bottom part where only one of s1 and s2 is still readable, or none at all + compute_unsafe_ssd_integral_image(ii, ii_linesize_32, + 0, endy_safe, + src, linesize, + offx, offy, e, w, h, + ii_w, ii_h - endy_safe); +} + +static int config_input(AVFilterLink *inlink) +{ + AVFilterContext *ctx = inlink->dst; + NLMeansContext *s = ctx->priv; + const AVPixFmtDescriptor *desc = av_pix_fmt_desc_get(inlink->format); + const int e = s->research_hsize + s->patch_hsize; + + s->chroma_w = FF_CEIL_RSHIFT(inlink->w, desc->log2_chroma_w); + s->chroma_h = FF_CEIL_RSHIFT(inlink->h, desc->log2_chroma_h); + s->nb_planes = av_pix_fmt_count_planes(inlink->format); + + /* Allocate the integral image with extra edges of thickness "e" + * + * +_+-------------------------------+ + * |0|0000000000000000000000000000000| + * +-x-------------------------------+ + * |0|\ ^ | + * |0| ii | e | + * |0| v | + * |0| +-----------------------+ | + * |0| | | | + * |0|<->| | | + * |0| e | | | + * |0| | | | + * |0| +-----------------------+ | + * |0| | + * |0| | + * |0| | + * +-+-------------------------------+ + */ + s->ii_w = inlink->w + e*2; + s->ii_h = inlink->h + e*2; + + // align to 4 the linesize, "+1" is for the space of the left 0-column + s->ii_lz_32 = FFALIGN(s->ii_w + 1, 4); + + // "+1" is for the space of the top 0-line + s->ii_orig = av_mallocz_array(s->ii_h + 1, s->ii_lz_32 * sizeof(*s->ii_orig)); + if (!s->ii_orig) + return AVERROR(ENOMEM); + + // skip top 0-line and left 0-column + s->ii = s->ii_orig + s->ii_lz_32 + 1; + + // allocate weighted average for every pixel + s->wa = av_malloc_array(inlink->w, inlink->h * sizeof(*s->wa)); + if (!s->wa) + return AVERROR(ENOMEM); + + return 0; +} + +struct thread_data { + const uint8_t *src; + int src_linesize; + int startx, starty; + int endx, endy; + const uint32_t *ii_start; + int wa_linesize; +}; + +static int nlmeans_slice(AVFilterContext *ctx, void *arg, int jobnr, int nb_jobs) +{ + int x, y; + NLMeansContext *s = ctx->priv; + const struct thread_data *td = arg; + const uint8_t *src = td->src; + const int src_linesize = td->src_linesize; + const int process_h = td->endy - td->starty; + const int slice_start = (process_h * jobnr ) / nb_jobs; + const int slice_end = (process_h * (jobnr+1)) / nb_jobs; + const int starty = td->starty + slice_start; + const int endy = td->starty + slice_end; + const int p = s->patch_hsize; + + for (y = starty; y < endy; y++) { + for (x = td->startx; x < td->endx; x++) { + const int patch_diff_sq = get_ssd_patch(td->ii_start, s->ii_lz_32, x, y, p); + if (patch_diff_sq < s->max_meaningful_diff) { + struct weighted_avg *wa = &s->wa[y*td->wa_linesize + x]; + const int weight_lut_idx = patch_diff_sq * s->pdiff_lut_scale; + const double weight = s->weight_lut[weight_lut_idx]; // exp(-patch_diff_sq * s->pdiff_scale) + wa->total_weight += weight; + wa->sum += weight * src[y*src_linesize + x]; + } + } + } + return 0; +} + +static int nlmeans_plane(AVFilterContext *ctx, int w, int h, + uint8_t *dst, int dst_linesize, + const uint8_t *src, int src_linesize) +{ + int x, y; + int offx, offy; + NLMeansContext *s = ctx->priv; + const int r = s->research_hsize; + /* patches center points cover the whole research window so the patches + * themselves overflow the research window */ + const int e = r + s->patch_hsize; + const uint32_t *centered_ii = s->ii + e*s->ii_lz_32 + e; + + memset(s->wa, 0, w * h * sizeof(*s->wa)); + + for (offy = -r; offy <= r; offy++) { + for (offx = -r; offx <= r; offx++) { + if (offx || offy) { + struct thread_data td = { + .src = src + offy*src_linesize + offx, + .src_linesize = src_linesize, + .startx = FFMAX(0, -offx), + .starty = FFMAX(0, -offy), + .endx = FFMIN(w, w - offx), + .endy = FFMIN(h, h - offy), + .ii_start = centered_ii + offy*s->ii_lz_32 + offx, + .wa_linesize = w, + }; + + compute_ssd_integral_image(s->ii, s->ii_lz_32, + src, src_linesize, + offx, offy, e, w, h); + ctx->internal->execute(ctx, nlmeans_slice, &td, NULL, + FFMIN(td.endy - td.starty, ff_filter_get_nb_threads(ctx))); + } + } + } + for (y = 0; y < h; y++) { + for (x = 0; x < w; x++) { + struct weighted_avg *wa = &s->wa[y*w + x]; + + // Also weight the centered pixel + wa->total_weight += 1.0; + wa->sum += 1.0 * src[y*src_linesize + x]; + + dst[y*dst_linesize + x] = av_clip_uint8(wa->sum / wa->total_weight); + } + } + return 0; +} + +static int filter_frame(AVFilterLink *inlink, AVFrame *in) +{ + int p; + AVFilterContext *ctx = inlink->dst; + NLMeansContext *s = ctx->priv; + AVFilterLink *outlink = ctx->outputs[0]; + + AVFrame *out = ff_get_video_buffer(outlink, outlink->w, outlink->h); + if (!out) { + av_frame_free(&in); + return AVERROR(ENOMEM); + } + av_frame_copy_props(out, in); + + for (p = 0; p < s->nb_planes; p++) { + const int w = p ? s->chroma_w : inlink->w; + const int h = p ? s->chroma_h : inlink->h; + nlmeans_plane(ctx, w, h, + out->data[p], out->linesize[p], + in->data[p], in->linesize[p]); + } + + av_frame_free(&in); + return ff_filter_frame(outlink, out); +} + +static av_cold int init(AVFilterContext *ctx) +{ + int i; + NLMeansContext *s = ctx->priv; + const double h = s->sigma * 10.; + + if (!(s->research_size & 1)) { + s->research_size |= 1; + av_log(ctx, AV_LOG_WARNING, "Research window size must be odd, setting it to %d\n", + s->research_size); + } + + if (!(s->patch_size & 1)) { + s->patch_size |= 1; + av_log(ctx, AV_LOG_WARNING, "Patch size must be odd, setting it to %d\n", + s->patch_size); + } + + s->pdiff_scale = 1. / (h * h); + s->max_meaningful_diff = -log(1/255.) / s->pdiff_scale; + s->pdiff_lut_scale = 1./s->max_meaningful_diff * WEIGHT_LUT_SIZE; + av_assert0((s->max_meaningful_diff - 1) * s->pdiff_lut_scale < FF_ARRAY_ELEMS(s->weight_lut)); + for (i = 0; i < WEIGHT_LUT_SIZE; i++) + s->weight_lut[i] = exp(-i / s->pdiff_lut_scale * s->pdiff_scale); + + s->research_hsize = s->research_size / 2; + s->patch_hsize = s->patch_size / 2; + + av_log(ctx, AV_LOG_INFO, "Research window: %dx%d, patch size: %dx%d\n", + s->research_size, s->research_size, s->patch_size, s->patch_size); + + return 0; +} + +static av_cold void uninit(AVFilterContext *ctx) +{ + NLMeansContext *s = ctx->priv; + av_freep(&s->ii_orig); +} + +static const AVFilterPad nlmeans_inputs[] = { + { + .name = "default", + .type = AVMEDIA_TYPE_VIDEO, + .config_props = config_input, + .filter_frame = filter_frame, + }, + { NULL } +}; + +static const AVFilterPad nlmeans_outputs[] = { + { + .name = "default", + .type = AVMEDIA_TYPE_VIDEO, + }, + { NULL } +}; + +AVFilter ff_vf_nlmeans = { + .name = "nlmeans", + .description = NULL_IF_CONFIG_SMALL("Non-locals mean denoiser."), + .priv_size = sizeof(NLMeansContext), + .init = init, + .uninit = uninit, + .query_formats = query_formats, + .inputs = nlmeans_inputs, + .outputs = nlmeans_outputs, + .priv_class = &nlmeans_class, + .flags = AVFILTER_FLAG_SUPPORT_TIMELINE_GENERIC | AVFILTER_FLAG_SLICE_THREADS, +};