Merge commit 'cebdedca57d95834a8f8098c7b6a322a1163e26b'
[ffmpeg.git] / 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         for (x = 0; x < BSIZE; x++) {
86             *line *= x == 0 ? 1. / sqrt(BSIZE) : sqrt(2. / BSIZE);
87             *column = *line++;
88             column += BSIZE;
89         }
90     }
91
92     column = ctx->tmp_block;
93     for (x = 0; x < BSIZE; x++) {
94         av_dct_calc(ctx->dct, column);
95         for (y = 0; y < BSIZE; y++)
96             column[y] *= y == 0 ? 1. / sqrt(BSIZE) : sqrt(2. / BSIZE);
97         column += BSIZE;
98     }
99
100     for (y = 0; y < BSIZE; y++)
101         for (x = 0; x < BSIZE; x++)
102             ctx->block[y*BSIZE + x] = ctx->tmp_block[x*BSIZE + y];
103
104     return ctx->block;
105 }
106
107 static void idct_block(DCTdnoizContext *ctx, float *dst, int dst_linesize)
108 {
109     int x, y;
110     float *block = ctx->block;
111     float *tmp = ctx->tmp_block;
112
113     for (y = 0; y < BSIZE; y++) {
114         for (x = 0; x < BSIZE; x++)
115             block[x] *= x == 0 ? sqrt(BSIZE) : 1./sqrt(2. / BSIZE);
116         av_dct_calc(ctx->idct, block);
117         block += BSIZE;
118     }
119
120     block = ctx->block;
121     for (y = 0; y < BSIZE; y++) {
122         for (x = 0; x < BSIZE; x++) {
123             tmp[x] = block[x*BSIZE + y];
124             tmp[x] *= x == 0 ? sqrt(BSIZE) : 1./sqrt(2. / BSIZE);
125         }
126         av_dct_calc(ctx->idct, tmp);
127         for (x = 0; x < BSIZE; x++)
128             dst[x*dst_linesize + y] += tmp[x];
129     }
130 }
131
132 static int config_input(AVFilterLink *inlink)
133 {
134     AVFilterContext *ctx = inlink->dst;
135     DCTdnoizContext *s = ctx->priv;
136     int i, x, y, bx, by, linesize, *iweights;
137     const float dct_3x3[3][3] = {
138         { 1./sqrt(3),  1./sqrt(3),  1./sqrt(3) },
139         { 1./sqrt(2),           0, -1./sqrt(2) },
140         { 1./sqrt(6), -2./sqrt(6),  1./sqrt(6) },
141     };
142     uint8_t rgba_map[4];
143
144     ff_fill_rgba_map(rgba_map, inlink->format);
145     for (y = 0; y < 3; y++)
146         for (x = 0; x < 3; x++)
147             s->color_dct[y][x] = dct_3x3[rgba_map[y]][rgba_map[x]];
148
149     s->pr_width  = inlink->w - (inlink->w - BSIZE) % s->step;
150     s->pr_height = inlink->h - (inlink->h - BSIZE) % s->step;
151     if (s->pr_width != inlink->w)
152         av_log(ctx, AV_LOG_WARNING, "The last %d horizontal pixels won't be denoised\n",
153                inlink->w - s->pr_width);
154     if (s->pr_height != inlink->h)
155         av_log(ctx, AV_LOG_WARNING, "The last %d vertical pixels won't be denoised\n",
156                inlink->h - s->pr_height);
157
158     s->p_linesize = linesize = FFALIGN(s->pr_width, 32);
159     for (i = 0; i < 2; i++) {
160         s->cbuf[i][0] = av_malloc(linesize * s->pr_height * sizeof(*s->cbuf[i][0]));
161         s->cbuf[i][1] = av_malloc(linesize * s->pr_height * sizeof(*s->cbuf[i][1]));
162         s->cbuf[i][2] = av_malloc(linesize * s->pr_height * sizeof(*s->cbuf[i][2]));
163         if (!s->cbuf[i][0] || !s->cbuf[i][1] || !s->cbuf[i][2])
164             return AVERROR(ENOMEM);
165     }
166
167     s->weights = av_malloc(s->pr_height * linesize * sizeof(*s->weights));
168     if (!s->weights)
169         return AVERROR(ENOMEM);
170     iweights = av_calloc(s->pr_height, linesize * sizeof(*iweights));
171     if (!iweights)
172         return AVERROR(ENOMEM);
173     for (y = 0; y < s->pr_height - BSIZE + 1; y += s->step)
174         for (x = 0; x < s->pr_width - BSIZE + 1; x += s->step)
175             for (by = 0; by < BSIZE; by++)
176                 for (bx = 0; bx < BSIZE; bx++)
177                     iweights[(y + by)*linesize + x + bx]++;
178     for (y = 0; y < s->pr_height; y++)
179         for (x = 0; x < s->pr_width; x++)
180             s->weights[y*linesize + x] = 1. / iweights[y*linesize + x];
181     av_free(iweights);
182
183     return 0;
184 }
185
186 static av_cold int init(AVFilterContext *ctx)
187 {
188     DCTdnoizContext *s = ctx->priv;
189
190     if (s->expr_str) {
191         int ret = av_expr_parse(&s->expr, s->expr_str, var_names,
192                                 NULL, NULL, NULL, NULL, 0, ctx);
193         if (ret < 0)
194             return ret;
195     }
196
197     s->th   = s->sigma * 3.;
198     s->step = BSIZE - s->overlap;
199     s->dct  = av_dct_init(NBITS, DCT_II);
200     s->idct = av_dct_init(NBITS, DCT_III);
201     s->block     = av_malloc(BSIZE * BSIZE * sizeof(*s->block));
202     s->tmp_block = av_malloc(BSIZE * BSIZE * sizeof(*s->tmp_block));
203
204     if (!s->dct || !s->idct || !s->tmp_block || !s->block)
205         return AVERROR(ENOMEM);
206
207     return 0;
208 }
209
210 static int query_formats(AVFilterContext *ctx)
211 {
212     static const enum AVPixelFormat pix_fmts[] = {
213         AV_PIX_FMT_BGR24, AV_PIX_FMT_RGB24,
214         AV_PIX_FMT_NONE
215     };
216     ff_set_common_formats(ctx, ff_make_format_list(pix_fmts));
217     return 0;
218 }
219
220 static void color_decorrelation(float dct3ch[3][3], float **dst, int dst_linesize,
221                                 const uint8_t *src, int src_linesize, int w, int h)
222 {
223     int x, y;
224     float *dstp_r = dst[0];
225     float *dstp_g = dst[1];
226     float *dstp_b = dst[2];
227
228     for (y = 0; y < h; y++) {
229         const uint8_t *srcp = src;
230
231         for (x = 0; x < w; x++) {
232             dstp_r[x] = srcp[0] * dct3ch[0][0] + srcp[1] * dct3ch[0][1] + srcp[2] * dct3ch[0][2];
233             dstp_g[x] = srcp[0] * dct3ch[1][0] + srcp[1] * dct3ch[1][1] + srcp[2] * dct3ch[1][2];
234             dstp_b[x] = srcp[0] * dct3ch[2][0] + srcp[1] * dct3ch[2][1] + srcp[2] * dct3ch[2][2];
235             srcp += 3;
236         }
237         src += src_linesize;
238         dstp_r += dst_linesize;
239         dstp_g += dst_linesize;
240         dstp_b += dst_linesize;
241     }
242 }
243
244 static void color_correlation(float dct3ch[3][3], uint8_t *dst, int dst_linesize,
245                               float **src, int src_linesize, int w, int h)
246 {
247     int x, y;
248     const float *src_r = src[0];
249     const float *src_g = src[1];
250     const float *src_b = src[2];
251
252     for (y = 0; y < h; y++) {
253         uint8_t *dstp = dst;
254
255         for (x = 0; x < w; x++) {
256             dstp[0] = av_clip_uint8(src_r[x] * dct3ch[0][0] + src_g[x] * dct3ch[1][0] + src_b[x] * dct3ch[2][0]);
257             dstp[1] = av_clip_uint8(src_r[x] * dct3ch[0][1] + src_g[x] * dct3ch[1][1] + src_b[x] * dct3ch[2][1]);
258             dstp[2] = av_clip_uint8(src_r[x] * dct3ch[0][2] + src_g[x] * dct3ch[1][2] + src_b[x] * dct3ch[2][2]);
259             dstp += 3;
260         }
261         dst += dst_linesize;
262         src_r += src_linesize;
263         src_g += src_linesize;
264         src_b += src_linesize;
265     }
266 }
267
268 static void filter_plane(AVFilterContext *ctx,
269                          float *dst, int dst_linesize,
270                          const float *src, int src_linesize,
271                          int w, int h)
272 {
273     int x, y, bx, by;
274     DCTdnoizContext *s = ctx->priv;
275     float *dst0 = dst;
276     const float *weights = s->weights;
277
278     // reset block sums
279     memset(dst, 0, h * dst_linesize * sizeof(*dst));
280
281     // block dct sums
282     for (y = 0; y < h - BSIZE + 1; y += s->step) {
283         for (x = 0; x < w - BSIZE + 1; x += s->step) {
284             float *ftb = dct_block(s, src + x, src_linesize);
285
286             if (s->expr) {
287                 for (by = 0; by < BSIZE; by++) {
288                     for (bx = 0; bx < BSIZE; bx++) {
289                         s->var_values[VAR_C] = FFABS(*ftb);
290                         *ftb++ *= av_expr_eval(s->expr, s->var_values, s);
291                     }
292                 }
293             } else {
294                 for (by = 0; by < BSIZE; by++) {
295                     for (bx = 0; bx < BSIZE; bx++) {
296                         if (FFABS(*ftb) < s->th)
297                             *ftb = 0;
298                         ftb++;
299                     }
300                 }
301             }
302             idct_block(s, dst + x, dst_linesize);
303         }
304         src += s->step * src_linesize;
305         dst += s->step * dst_linesize;
306     }
307
308     // average blocks
309     dst = dst0;
310     for (y = 0; y < h; y++) {
311         for (x = 0; x < w; x++)
312             dst[x] *= weights[x];
313         dst += dst_linesize;
314         weights += dst_linesize;
315     }
316 }
317
318 static int filter_frame(AVFilterLink *inlink, AVFrame *in)
319 {
320     AVFilterContext *ctx = inlink->dst;
321     DCTdnoizContext *s = ctx->priv;
322     AVFilterLink *outlink = inlink->dst->outputs[0];
323     int direct, plane;
324     AVFrame *out;
325
326     if (av_frame_is_writable(in)) {
327         direct = 1;
328         out = in;
329     } else {
330         direct = 0;
331         out = ff_get_video_buffer(outlink, outlink->w, outlink->h);
332         if (!out) {
333             av_frame_free(&in);
334             return AVERROR(ENOMEM);
335         }
336         av_frame_copy_props(out, in);
337     }
338
339     color_decorrelation(s->color_dct, s->cbuf[0], s->p_linesize,
340                         in->data[0], in->linesize[0], s->pr_width, s->pr_height);
341     for (plane = 0; plane < 3; plane++)
342         filter_plane(ctx, s->cbuf[1][plane], s->p_linesize,
343                           s->cbuf[0][plane], s->p_linesize,
344                           s->pr_width, s->pr_height);
345     color_correlation(s->color_dct, out->data[0], out->linesize[0],
346                       s->cbuf[1], s->p_linesize, s->pr_width, s->pr_height);
347
348     if (!direct) {
349         int y;
350         uint8_t *dst = out->data[0];
351         const uint8_t *src = in->data[0];
352         const int dst_linesize = out->linesize[0];
353         const int src_linesize = in->linesize[0];
354         const int hpad = (inlink->w - s->pr_width) * 3;
355         const int vpad = (inlink->h - s->pr_height);
356
357         if (hpad) {
358             uint8_t       *dstp = dst + s->pr_width * 3;
359             const uint8_t *srcp = src + s->pr_width * 3;
360
361             for (y = 0; y < s->pr_height; y++) {
362                 memcpy(dstp, srcp, hpad);
363                 dstp += dst_linesize;
364                 srcp += src_linesize;
365             }
366         }
367         if (vpad) {
368             uint8_t       *dstp = dst + s->pr_height * dst_linesize;
369             const uint8_t *srcp = src + s->pr_height * src_linesize;
370
371             for (y = 0; y < vpad; y++) {
372                 memcpy(dstp, srcp, inlink->w * 3);
373                 dstp += dst_linesize;
374                 srcp += src_linesize;
375             }
376         }
377
378         av_frame_free(&in);
379     }
380
381     return ff_filter_frame(outlink, out);
382 }
383
384 static av_cold void uninit(AVFilterContext *ctx)
385 {
386     int i;
387     DCTdnoizContext *s = ctx->priv;
388
389     av_dct_end(s->dct);
390     av_dct_end(s->idct);
391     av_free(s->block);
392     av_free(s->tmp_block);
393     av_free(s->weights);
394     for (i = 0; i < 2; i++) {
395         av_free(s->cbuf[i][0]);
396         av_free(s->cbuf[i][1]);
397         av_free(s->cbuf[i][2]);
398     }
399     av_expr_free(s->expr);
400 }
401
402 static const AVFilterPad dctdnoiz_inputs[] = {
403     {
404         .name         = "default",
405         .type         = AVMEDIA_TYPE_VIDEO,
406         .filter_frame = filter_frame,
407         .config_props = config_input,
408     },
409     { NULL }
410 };
411
412 static const AVFilterPad dctdnoiz_outputs[] = {
413     {
414         .name = "default",
415         .type = AVMEDIA_TYPE_VIDEO,
416     },
417     { NULL }
418 };
419
420 AVFilter avfilter_vf_dctdnoiz = {
421     .name          = "dctdnoiz",
422     .description   = NULL_IF_CONFIG_SMALL("Denoise frames using 2D DCT."),
423     .priv_size     = sizeof(DCTdnoizContext),
424     .init          = init,
425     .uninit        = uninit,
426     .query_formats = query_formats,
427     .inputs        = dctdnoiz_inputs,
428     .outputs       = dctdnoiz_outputs,
429     .priv_class    = &dctdnoiz_class,
430     .flags         = AVFILTER_FLAG_SUPPORT_TIMELINE_GENERIC,
431 };