git.sesse.net Git - ffmpeg/blob - libavfilter/vf_dctdnoiz.c

   1 /*
   2  * Copyright (c) 2013 Clément Bœsch
   3  *
   4  * This file is part of FFmpeg.
   5  *
   6  * FFmpeg is free software; you can redistribute it and/or
   7  * modify it under the terms of the GNU Lesser General Public
   8  * License as published by the Free Software Foundation; either
   9  * version 2.1 of the License, or (at your option) any later version.
  10  *
  11  * FFmpeg is distributed in the hope that it will be useful,
  12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14  * Lesser General Public License for more details.
  15  *
  16  * You should have received a copy of the GNU Lesser General Public
  17  * License along with FFmpeg; if not, write to the Free Software
  18  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
  19  */
  20
  21 /**
  22  * A simple, relatively efficient and extremely slow DCT image denoiser.
  23  * @see http://www.ipol.im/pub/art/2011/ys-dct/
  24  */
  25
  26 #include "libavcodec/avfft.h"
  27 #include "libavutil/eval.h"
  28 #include "libavutil/opt.h"
  29 #include "drawutils.h"
  30 #include "internal.h"
  31
  32 #define NBITS 4
  33 #define BSIZE (1<<(NBITS))
  34
  35 static const char *const var_names[] = { "c", NULL };
  36 enum { VAR_C, VAR_VARS_NB };
  37
  38 typedef struct {
  39     const AVClass *class;
  40
  41     /* coefficient factor expression */
  42     char *expr_str;
  43     AVExpr *expr;
  44     double var_values[VAR_VARS_NB];
  45
  46     int pr_width, pr_height;    // width and height to process
  47     float sigma;                // used when no expression are st
  48     float th;                   // threshold (3*sigma)
  49     float color_dct[3][3];      // 3x3 DCT for color decorrelation
  50     float *cbuf[2][3];          // two planar rgb color buffers
  51     float *weights;             // dct coeff are cumulated with overlapping; these values are used for averaging
  52     int p_linesize;             // line sizes for color and weights
  53     int overlap;                // number of block overlapping pixels
  54     int step;                   // block step increment (BSIZE - overlap)
  55     DCTContext *dct, *idct;     // DCT and inverse DCT contexts
  56     float *block, *tmp_block;   // two BSIZE x BSIZE block buffers
  57 } DCTdnoizContext;
  58
  59 #define OFFSET(x) offsetof(DCTdnoizContext, x)
  60 #define FLAGS AV_OPT_FLAG_FILTERING_PARAM|AV_OPT_FLAG_VIDEO_PARAM
  61 static const AVOption dctdnoiz_options[] = {
  62     { "sigma",   "set noise sigma constant",               OFFSET(sigma),    AV_OPT_TYPE_FLOAT,  {.dbl=0},            0, 999,          .flags = FLAGS },
  63     { "s",       "set noise sigma constant",               OFFSET(sigma),    AV_OPT_TYPE_FLOAT,  {.dbl=0},            0, 999,          .flags = FLAGS },
  64     { "overlap", "set number of block overlapping pixels", OFFSET(overlap),  AV_OPT_TYPE_INT,    {.i64=(1<<NBITS)-1}, 0, (1<<NBITS)-1, .flags = FLAGS },
  65     { "expr",    "set coefficient factor expression",      OFFSET(expr_str), AV_OPT_TYPE_STRING, {.str=NULL},                          .flags = FLAGS },
  66     { "e",       "set coefficient factor expression",      OFFSET(expr_str), AV_OPT_TYPE_STRING, {.str=NULL},                          .flags = FLAGS },
  67     { NULL }
  68 };
  69
  70 AVFILTER_DEFINE_CLASS(dctdnoiz);
  71
  72 static float *dct_block(DCTdnoizContext *ctx, const float *src, int src_linesize)
  73 {
  74     int x, y;
  75     float *column;
  76
  77     for (y = 0; y < BSIZE; y++) {
  78         float *line = ctx->block;
  79
  80         memcpy(line, src, BSIZE * sizeof(*line));
  81         src += src_linesize;
  82         av_dct_calc(ctx->dct, line);
  83
  84         column = ctx->tmp_block + y;
  85         column[0] = line[0] * (1. / sqrt(BSIZE));
  86         column += BSIZE;
  87         for (x = 1; x < BSIZE; x++) {
  88             *column = line[x] * sqrt(2. / BSIZE);
  89             column += BSIZE;
  90         }
  91     }
  92
  93     column = ctx->tmp_block;
  94     for (x = 0; x < BSIZE; x++) {
  95         av_dct_calc(ctx->dct, column);
  96         column[0] *= 1. / sqrt(BSIZE);
  97         for (y = 1; y < BSIZE; y++)
  98             column[y] *= sqrt(2. / BSIZE);
  99         column += BSIZE;
 100     }
 101
 102     for (y = 0; y < BSIZE; y++)
 103         for (x = 0; x < BSIZE; x++)
 104             ctx->block[y*BSIZE + x] = ctx->tmp_block[x*BSIZE + y];
 105
 106     return ctx->block;
 107 }
 108
 109 static void idct_block(DCTdnoizContext *ctx, float *dst, int dst_linesize)
 110 {
 111     int x, y;
 112     float *block = ctx->block;
 113     float *tmp = ctx->tmp_block;
 114
 115     for (y = 0; y < BSIZE; y++) {
 116         block[0] *= sqrt(BSIZE);
 117         for (x = 1; x < BSIZE; x++)
 118             block[x] *= 1./sqrt(2. / BSIZE);
 119         av_dct_calc(ctx->idct, block);
 120         block += BSIZE;
 121     }
 122
 123     block = ctx->block;
 124     for (y = 0; y < BSIZE; y++) {
 125         tmp[0] = block[y] * sqrt(BSIZE);
 126         for (x = 1; x < BSIZE; x++)
 127             tmp[x] = block[x*BSIZE + y] * (1./sqrt(2. / BSIZE));
 128         av_dct_calc(ctx->idct, tmp);
 129         for (x = 0; x < BSIZE; x++)
 130             dst[x*dst_linesize + y] += tmp[x];
 131     }
 132 }
 133
 134 static int config_input(AVFilterLink *inlink)
 135 {
 136     AVFilterContext *ctx = inlink->dst;
 137     DCTdnoizContext *s = ctx->priv;
 138     int i, x, y, bx, by, linesize, *iweights;
 139     const float dct_3x3[3][3] = {
 140         { 1./sqrt(3),  1./sqrt(3),  1./sqrt(3) },
 141         { 1./sqrt(2),           0, -1./sqrt(2) },
 142         { 1./sqrt(6), -2./sqrt(6),  1./sqrt(6) },
 143     };
 144     uint8_t rgba_map[4];
 145
 146     ff_fill_rgba_map(rgba_map, inlink->format);
 147     for (y = 0; y < 3; y++)
 148         for (x = 0; x < 3; x++)
 149             s->color_dct[y][x] = dct_3x3[rgba_map[y]][rgba_map[x]];
 150
 151     s->pr_width  = inlink->w - (inlink->w - BSIZE) % s->step;
 152     s->pr_height = inlink->h - (inlink->h - BSIZE) % s->step;
 153     if (s->pr_width != inlink->w)
 154         av_log(ctx, AV_LOG_WARNING, "The last %d horizontal pixels won't be denoised\n",
 155                inlink->w - s->pr_width);
 156     if (s->pr_height != inlink->h)
 157         av_log(ctx, AV_LOG_WARNING, "The last %d vertical pixels won't be denoised\n",
 158                inlink->h - s->pr_height);
 159
 160     s->p_linesize = linesize = FFALIGN(s->pr_width, 32);
 161     for (i = 0; i < 2; i++) {
 162         s->cbuf[i][0] = av_malloc(linesize * s->pr_height * sizeof(*s->cbuf[i][0]));
 163         s->cbuf[i][1] = av_malloc(linesize * s->pr_height * sizeof(*s->cbuf[i][1]));
 164         s->cbuf[i][2] = av_malloc(linesize * s->pr_height * sizeof(*s->cbuf[i][2]));
 165         if (!s->cbuf[i][0] || !s->cbuf[i][1] || !s->cbuf[i][2])
 166             return AVERROR(ENOMEM);
 167     }
 168
 169     s->weights = av_malloc(s->pr_height * linesize * sizeof(*s->weights));
 170     if (!s->weights)
 171         return AVERROR(ENOMEM);
 172     iweights = av_calloc(s->pr_height, linesize * sizeof(*iweights));
 173     if (!iweights)
 174         return AVERROR(ENOMEM);
 175     for (y = 0; y < s->pr_height - BSIZE + 1; y += s->step)
 176         for (x = 0; x < s->pr_width - BSIZE + 1; x += s->step)
 177             for (by = 0; by < BSIZE; by++)
 178                 for (bx = 0; bx < BSIZE; bx++)
 179                     iweights[(y + by)*linesize + x + bx]++;
 180     for (y = 0; y < s->pr_height; y++)
 181         for (x = 0; x < s->pr_width; x++)
 182             s->weights[y*linesize + x] = 1. / iweights[y*linesize + x];
 183     av_free(iweights);
 184
 185     return 0;
 186 }
 187
 188 static av_cold int init(AVFilterContext *ctx)
 189 {
 190     DCTdnoizContext *s = ctx->priv;
 191
 192     if (s->expr_str) {
 193         int ret = av_expr_parse(&s->expr, s->expr_str, var_names,
 194                                 NULL, NULL, NULL, NULL, 0, ctx);
 195         if (ret < 0)
 196             return ret;
 197     }
 198
 199     s->th   = s->sigma * 3.;
 200     s->step = BSIZE - s->overlap;
 201     s->dct  = av_dct_init(NBITS, DCT_II);
 202     s->idct = av_dct_init(NBITS, DCT_III);
 203     s->block     = av_malloc(BSIZE * BSIZE * sizeof(*s->block));
 204     s->tmp_block = av_malloc(BSIZE * BSIZE * sizeof(*s->tmp_block));
 205
 206     if (!s->dct || !s->idct || !s->tmp_block || !s->block)
 207         return AVERROR(ENOMEM);
 208
 209     return 0;
 210 }
 211
 212 static int query_formats(AVFilterContext *ctx)
 213 {
 214     static const enum AVPixelFormat pix_fmts[] = {
 215         AV_PIX_FMT_BGR24, AV_PIX_FMT_RGB24,
 216         AV_PIX_FMT_NONE
 217     };
 218     ff_set_common_formats(ctx, ff_make_format_list(pix_fmts));
 219     return 0;
 220 }
 221
 222 static void color_decorrelation(float dct3ch[3][3], float **dst, int dst_linesize,
 223                                 const uint8_t *src, int src_linesize, int w, int h)
 224 {
 225     int x, y;
 226     float *dstp_r = dst[0];
 227     float *dstp_g = dst[1];
 228     float *dstp_b = dst[2];
 229
 230     for (y = 0; y < h; y++) {
 231         const uint8_t *srcp = src;
 232
 233         for (x = 0; x < w; x++) {
 234             dstp_r[x] = srcp[0] * dct3ch[0][0] + srcp[1] * dct3ch[0][1] + srcp[2] * dct3ch[0][2];
 235             dstp_g[x] = srcp[0] * dct3ch[1][0] + srcp[1] * dct3ch[1][1] + srcp[2] * dct3ch[1][2];
 236             dstp_b[x] = srcp[0] * dct3ch[2][0] + srcp[1] * dct3ch[2][1] + srcp[2] * dct3ch[2][2];
 237             srcp += 3;
 238         }
 239         src += src_linesize;
 240         dstp_r += dst_linesize;
 241         dstp_g += dst_linesize;
 242         dstp_b += dst_linesize;
 243     }
 244 }
 245
 246 static void color_correlation(float dct3ch[3][3], uint8_t *dst, int dst_linesize,
 247                               float **src, int src_linesize, int w, int h)
 248 {
 249     int x, y;
 250     const float *src_r = src[0];
 251     const float *src_g = src[1];
 252     const float *src_b = src[2];
 253
 254     for (y = 0; y < h; y++) {
 255         uint8_t *dstp = dst;
 256
 257         for (x = 0; x < w; x++) {
 258             dstp[0] = av_clip_uint8(src_r[x] * dct3ch[0][0] + src_g[x] * dct3ch[1][0] + src_b[x] * dct3ch[2][0]);
 259             dstp[1] = av_clip_uint8(src_r[x] * dct3ch[0][1] + src_g[x] * dct3ch[1][1] + src_b[x] * dct3ch[2][1]);
 260             dstp[2] = av_clip_uint8(src_r[x] * dct3ch[0][2] + src_g[x] * dct3ch[1][2] + src_b[x] * dct3ch[2][2]);
 261             dstp += 3;
 262         }
 263         dst += dst_linesize;
 264         src_r += src_linesize;
 265         src_g += src_linesize;
 266         src_b += src_linesize;
 267     }
 268 }
 269
 270 static void filter_plane(AVFilterContext *ctx,
 271                          float *dst, int dst_linesize,
 272                          const float *src, int src_linesize,
 273                          int w, int h)
 274 {
 275     int x, y, bx, by;
 276     DCTdnoizContext *s = ctx->priv;
 277     float *dst0 = dst;
 278     const float *weights = s->weights;
 279
 280     // reset block sums
 281     memset(dst, 0, h * dst_linesize * sizeof(*dst));
 282
 283     // block dct sums
 284     for (y = 0; y < h - BSIZE + 1; y += s->step) {
 285         for (x = 0; x < w - BSIZE + 1; x += s->step) {
 286             float *ftb = dct_block(s, src + x, src_linesize);
 287
 288             if (s->expr) {
 289                 for (by = 0; by < BSIZE; by++) {
 290                     for (bx = 0; bx < BSIZE; bx++) {
 291                         s->var_values[VAR_C] = FFABS(*ftb);
 292                         *ftb++ *= av_expr_eval(s->expr, s->var_values, s);
 293                     }
 294                 }
 295             } else {
 296                 for (by = 0; by < BSIZE; by++) {
 297                     for (bx = 0; bx < BSIZE; bx++) {
 298                         if (FFABS(*ftb) < s->th)
 299                             *ftb = 0;
 300                         ftb++;
 301                     }
 302                 }
 303             }
 304             idct_block(s, dst + x, dst_linesize);
 305         }
 306         src += s->step * src_linesize;
 307         dst += s->step * dst_linesize;
 308     }
 309
 310     // average blocks
 311     dst = dst0;
 312     for (y = 0; y < h; y++) {
 313         for (x = 0; x < w; x++)
 314             dst[x] *= weights[x];
 315         dst += dst_linesize;
 316         weights += dst_linesize;
 317     }
 318 }
 319
 320 static int filter_frame(AVFilterLink *inlink, AVFrame *in)
 321 {
 322     AVFilterContext *ctx = inlink->dst;
 323     DCTdnoizContext *s = ctx->priv;
 324     AVFilterLink *outlink = inlink->dst->outputs[0];
 325     int direct, plane;
 326     AVFrame *out;
 327
 328     if (av_frame_is_writable(in)) {
 329         direct = 1;
 330         out = in;
 331     } else {
 332         direct = 0;
 333         out = ff_get_video_buffer(outlink, outlink->w, outlink->h);
 334         if (!out) {
 335             av_frame_free(&in);
 336             return AVERROR(ENOMEM);
 337         }
 338         av_frame_copy_props(out, in);
 339     }
 340
 341     color_decorrelation(s->color_dct, s->cbuf[0], s->p_linesize,
 342                         in->data[0], in->linesize[0], s->pr_width, s->pr_height);
 343     for (plane = 0; plane < 3; plane++)
 344         filter_plane(ctx, s->cbuf[1][plane], s->p_linesize,
 345                           s->cbuf[0][plane], s->p_linesize,
 346                           s->pr_width, s->pr_height);
 347     color_correlation(s->color_dct, out->data[0], out->linesize[0],
 348                       s->cbuf[1], s->p_linesize, s->pr_width, s->pr_height);
 349
 350     if (!direct) {
 351         int y;
 352         uint8_t *dst = out->data[0];
 353         const uint8_t *src = in->data[0];
 354         const int dst_linesize = out->linesize[0];
 355         const int src_linesize = in->linesize[0];
 356         const int hpad = (inlink->w - s->pr_width) * 3;
 357         const int vpad = (inlink->h - s->pr_height);
 358
 359         if (hpad) {
 360             uint8_t       *dstp = dst + s->pr_width * 3;
 361             const uint8_t *srcp = src + s->pr_width * 3;
 362
 363             for (y = 0; y < s->pr_height; y++) {
 364                 memcpy(dstp, srcp, hpad);
 365                 dstp += dst_linesize;
 366                 srcp += src_linesize;
 367             }
 368         }
 369         if (vpad) {
 370             uint8_t       *dstp = dst + s->pr_height * dst_linesize;
 371             const uint8_t *srcp = src + s->pr_height * src_linesize;
 372
 373             for (y = 0; y < vpad; y++) {
 374                 memcpy(dstp, srcp, inlink->w * 3);
 375                 dstp += dst_linesize;
 376                 srcp += src_linesize;
 377             }
 378         }
 379
 380         av_frame_free(&in);
 381     }
 382
 383     return ff_filter_frame(outlink, out);
 384 }
 385
 386 static av_cold void uninit(AVFilterContext *ctx)
 387 {
 388     int i;
 389     DCTdnoizContext *s = ctx->priv;
 390
 391     av_dct_end(s->dct);
 392     av_dct_end(s->idct);
 393     av_free(s->block);
 394     av_free(s->tmp_block);
 395     av_free(s->weights);
 396     for (i = 0; i < 2; i++) {
 397         av_free(s->cbuf[i][0]);
 398         av_free(s->cbuf[i][1]);
 399         av_free(s->cbuf[i][2]);
 400     }
 401     av_expr_free(s->expr);
 402 }
 403
 404 static const AVFilterPad dctdnoiz_inputs[] = {
 405     {
 406         .name         = "default",
 407         .type         = AVMEDIA_TYPE_VIDEO,
 408         .filter_frame = filter_frame,
 409         .config_props = config_input,
 410     },
 411     { NULL }
 412 };
 413
 414 static const AVFilterPad dctdnoiz_outputs[] = {
 415     {
 416         .name = "default",
 417         .type = AVMEDIA_TYPE_VIDEO,
 418     },
 419     { NULL }
 420 };
 421
 422 AVFilter ff_vf_dctdnoiz = {
 423     .name          = "dctdnoiz",
 424     .description   = NULL_IF_CONFIG_SMALL("Denoise frames using 2D DCT."),
 425     .priv_size     = sizeof(DCTdnoizContext),
 426     .init          = init,
 427     .uninit        = uninit,
 428     .query_formats = query_formats,
 429     .inputs        = dctdnoiz_inputs,
 430     .outputs       = dctdnoiz_outputs,
 431     .priv_class    = &dctdnoiz_class,
 432     .flags         = AVFILTER_FLAG_SUPPORT_TIMELINE_GENERIC,
 433 };