From c29b532a942233b5c804e620d475b3b888fe6e72 Mon Sep 17 00:00:00 2001 From: Clément Bœsch Date: Wed, 30 Sep 2015 20:29:30 +0200 Subject: lavfi: add nlmeans filter Fixes Ticket #4910 --- libavfilter/vf_nlmeans.c | 551 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 551 insertions(+) create mode 100644 libavfilter/vf_nlmeans.c (limited to 'libavfilter/vf_nlmeans.c') 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, +}; -- cgit v1.2.3