From c29b532a942233b5c804e620d475b3b888fe6e72 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Cl=C3=A9ment=20B=C5=93sch?= Date: Wed, 30 Sep 2015 20:29:30 +0200 Subject: [PATCH] lavfi: add nlmeans filter Fixes Ticket #4910 --- Changelog | 1 + doc/filters.texi | 35 +++ libavfilter/Makefile | 3 +- libavfilter/allfilters.c | 1 + libavfilter/tests/integral.c | 92 ++++++ libavfilter/version.h | 2 +- libavfilter/vf_nlmeans.c | 551 +++++++++++++++++++++++++++++++++++ 7 files changed, 683 insertions(+), 2 deletions(-) create mode 100644 libavfilter/tests/integral.c create mode 100644 libavfilter/vf_nlmeans.c diff --git a/Changelog b/Changelog index 2d0a449c13..a6916bef74 100644 --- a/Changelog +++ b/Changelog @@ -31,6 +31,7 @@ version : - MediaCodec HEVC decoding - TrueHD encoder - Meridian Lossless Packing (MLP) encoder +- Non-Local Means (nlmeans) denoising filter version 3.1: diff --git a/doc/filters.texi b/doc/filters.texi index 070e57d40f..7aa541f42b 100644 --- a/doc/filters.texi +++ b/doc/filters.texi @@ -9695,6 +9695,41 @@ Negate input video. It accepts an integer in input; if non-zero it negates the alpha component (if available). The default value in input is 0. +@section nlmeans + +Denoise frames using Non-Local Means algorithm. + +Each pixel is adjusted by looking for other pixels with similar contexts. This +context similarity is defined by comparing their surrounding patches of size +@option{p}x@option{p}. Patches are searched in an area of @option{r}x@option{r} +around the pixel. + +Note that the research area defines centers for patches, which means some +patches will be made of pixels outside that research area. + +The filter accepts the following options. + +@table @option +@item s +Set denoising strength. + +@item p +Set patch size. + +@item pc +Same as @option{p} but for chroma planes. + +The default value is @var{0} and means automatic. + +@item r +Set research size. + +@item rc +Same as @option{r} but for chroma planes. + +The default value is @var{0} and means automatic. +@end table + @section nnedi Deinterlace video using neural network edge directed interpolation. diff --git a/libavfilter/Makefile b/libavfilter/Makefile index 5cd10faabb..57a38d3754 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 47d95f5f6a..8f542fddb6 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 0000000000..8a2286b54d --- /dev/null +++ b/libavfilter/tests/integral.c @@ -0,0 +1,92 @@ +/* + * 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" + +static void display_integral(const uint32_t *ii, int w, int h, int lz_32) +{ + int x, y; + + for (y = 0; y < h; y++) { + for (x = 0; x < w; x++) + printf(" %7x", ii[y*lz_32 + x]); + printf("\n"); + } + printf("---------------\n"); +} + +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 e = 3; + const int ii_w = w+e*2, ii_h = h+e*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 = av_mallocz_array(ii_h + 1, ii_lz_32 * sizeof(*ii)); + uint32_t *ii2 = av_mallocz_array(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 + + if (!ii || !ii2) + return -1; + + for (yoff = -e; yoff <= e; yoff++) { + for (xoff = -e; xoff <= e; 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, e, w, h); + display_integral(ii_start, ii_w, ii_h, ii_lz_32); + + compute_unsafe_ssd_integral_image(ii_start2, ii_lz_32, + 0, 0, + src, lz, + xoff, yoff, e, w, h, + ii_w, ii_h); + display_integral(ii_start2, ii_w, ii_h, ii_lz_32); + + if (memcmp(ii, ii2, (ii_h+1) * ii_lz_32 * sizeof(*ii))) { + printf("Integral mismatch\n"); + ret = 1; + goto end; + } + } + } + +end: + av_freep(&ii); + av_freep(&ii2); + return ret; +} diff --git a/libavfilter/version.h b/libavfilter/version.h index 1a7a1b5686..93d249be10 100644 --- a/libavfilter/version.h +++ b/libavfilter/version.h @@ -30,7 +30,7 @@ #include "libavutil/version.h" #define LIBAVFILTER_VERSION_MAJOR 6 -#define LIBAVFILTER_VERSION_MINOR 62 +#define LIBAVFILTER_VERSION_MINOR 63 #define LIBAVFILTER_VERSION_MICRO 100 #define LIBAVFILTER_VERSION_INT AV_VERSION_INT(LIBAVFILTER_VERSION_MAJOR, \ diff --git a/libavfilter/vf_nlmeans.c b/libavfilter/vf_nlmeans.c new file mode 100644 index 0000000000..2487813504 --- /dev/null +++ b/libavfilter/vf_nlmeans.c @@ -0,0 +1,551 @@ +/* + * 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 + */ + +/** + * @todo + * - SIMD for compute_safe_ssd_integral_image + * - SIMD for final weighted averaging + * - better automatic defaults? see "Parameters" @ http://www.ipol.im/pub/art/2011/bcm_nlm/ + * - temporal support (probably doesn't need any displacement according to + * "Denoising image sequences does not require motion estimation") + * - Bayer pixel format support for at least raw photos? (DNG support would be + * handy here) + * - FATE test (probably needs visual threshold test mechanism due to the use + * of floats) + */ + +#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 = FFMAX(s->research_hsize, s->research_hsize_uv) + + FFMAX(s->patch_hsize, s->patch_hsize_uv); + + 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_linesize = inlink->w; + s->wa = av_malloc_array(s->wa_linesize, 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 p; +}; + +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; + + for (y = starty; y < endy; y++) { + for (x = td->startx; x < td->endx; x++) { + const int patch_diff_sq = get_integral_patch_value(td->ii_start, s->ii_lz_32, x, y, td->p); + if (patch_diff_sq < s->max_meaningful_diff) { + struct weighted_avg *wa = &s->wa[y*s->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, int p, int r, + uint8_t *dst, int dst_linesize, + const uint8_t *src, int src_linesize) +{ + int x, y; + int offx, offy; + NLMeansContext *s = ctx->priv; + /* patches center points cover the whole research window so the patches + * themselves overflow the research window */ + const int e = r + p; + /* focus an integral pointer on the centered image (s1) */ + const uint32_t *centered_ii = s->ii + e*s->ii_lz_32 + e; + + memset(s->wa, 0, s->wa_linesize * 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, + .p = p, + }; + + 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*s->wa_linesize + 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 i; + 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 (i = 0; i < s->nb_planes; i++) { + const int w = i ? s->chroma_w : inlink->w; + const int h = i ? s->chroma_h : inlink->h; + const int p = i ? s->patch_hsize_uv : s->patch_hsize; + const int r = i ? s->research_hsize_uv : s->research_hsize; + nlmeans_plane(ctx, w, h, p, r, + out->data[i], out->linesize[i], + in->data[i], in->linesize[i]); + } + + av_frame_free(&in); + return ff_filter_frame(outlink, out); +} + +#define CHECK_ODD_FIELD(field, name) do { \ + if (!(s->field & 1)) { \ + s->field |= 1; \ + av_log(ctx, AV_LOG_WARNING, name " size must be odd, " \ + "setting it to %d\n", s->field); \ + } \ +} while (0) + +static av_cold int init(AVFilterContext *ctx) +{ + int i; + NLMeansContext *s = ctx->priv; + const double h = s->sigma * 10.; + + 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); + + CHECK_ODD_FIELD(research_size, "Luma research window"); + CHECK_ODD_FIELD(patch_size, "Luma patch"); + + if (!s->research_size_uv) s->research_size_uv = s->research_size; + if (!s->patch_size_uv) s->patch_size_uv = s->patch_size; + + CHECK_ODD_FIELD(research_size_uv, "Chroma research window"); + CHECK_ODD_FIELD(patch_size_uv, "Chroma patch"); + + s->research_hsize = s->research_size / 2; + s->research_hsize_uv = s->research_size_uv / 2; + s->patch_hsize = s->patch_size / 2; + s->patch_hsize_uv = s->patch_size_uv / 2; + + av_log(ctx, AV_LOG_INFO, "Research window: %dx%d / %dx%d, patch size: %dx%d / %dx%d\n", + s->research_size, s->research_size, s->research_size_uv, s->research_size_uv, + s->patch_size, s->patch_size, s->patch_size_uv, s->patch_size_uv); + + return 0; +} + +static av_cold void uninit(AVFilterContext *ctx) +{ + NLMeansContext *s = ctx->priv; + av_freep(&s->ii_orig); + av_freep(&s->wa); +} + +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-local means 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, +};