FFmpeg
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Groups Pages
vf_dctdnoiz.c
Go to the documentation of this file.
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;
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
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 
213 {
214  static const enum AVPixelFormat pix_fmts[] = {
217  };
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 
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 
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,
429  .inputs = dctdnoiz_inputs,
430  .outputs = dctdnoiz_outputs,
431  .priv_class = &dctdnoiz_class,
433 };