2 * Copyright (c) 2013 Clément Bœsch
4 * This file is part of FFmpeg.
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.
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.
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
22 * A simple, relatively efficient and extremely slow DCT image denoiser.
23 * @see http://www.ipol.im/pub/art/2011/ys-dct/
26 #include "libavcodec/avfft.h"
27 #include "libavutil/eval.h"
28 #include "libavutil/opt.h"
29 #include "drawutils.h"
33 #define BSIZE (1<<(NBITS))
35 static const char *const var_names[] = { "c", NULL };
36 enum { VAR_C, VAR_VARS_NB };
41 /* coefficient factor expression */
44 double var_values[VAR_VARS_NB];
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
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 },
70 AVFILTER_DEFINE_CLASS(dctdnoiz);
72 static float *dct_block(DCTdnoizContext *ctx, const float *src, int src_linesize)
77 for (y = 0; y < BSIZE; y++) {
78 float *line = ctx->block;
80 memcpy(line, src, BSIZE * sizeof(*line));
82 av_dct_calc(ctx->dct, line);
84 column = ctx->tmp_block + y;
85 column[0] = line[0] * (1. / sqrt(BSIZE));
87 for (x = 1; x < BSIZE; x++) {
88 *column = line[x] * sqrt(2. / BSIZE);
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);
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];
109 static void idct_block(DCTdnoizContext *ctx, float *dst, int dst_linesize)
112 float *block = ctx->block;
113 float *tmp = ctx->tmp_block;
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);
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];
134 static int config_input(AVFilterLink *inlink)
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) },
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]];
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);
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);
169 s->weights = av_malloc(s->pr_height * linesize * sizeof(*s->weights));
171 return AVERROR(ENOMEM);
172 iweights = av_calloc(s->pr_height, linesize * sizeof(*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];
188 static av_cold int init(AVFilterContext *ctx)
190 DCTdnoizContext *s = ctx->priv;
193 int ret = av_expr_parse(&s->expr, s->expr_str, var_names,
194 NULL, NULL, NULL, NULL, 0, ctx);
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));
206 if (!s->dct || !s->idct || !s->tmp_block || !s->block)
207 return AVERROR(ENOMEM);
212 static int query_formats(AVFilterContext *ctx)
214 static const enum AVPixelFormat pix_fmts[] = {
215 AV_PIX_FMT_BGR24, AV_PIX_FMT_RGB24,
218 ff_set_common_formats(ctx, ff_make_format_list(pix_fmts));
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)
226 float *dstp_r = dst[0];
227 float *dstp_g = dst[1];
228 float *dstp_b = dst[2];
230 for (y = 0; y < h; y++) {
231 const uint8_t *srcp = src;
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];
240 dstp_r += dst_linesize;
241 dstp_g += dst_linesize;
242 dstp_b += dst_linesize;
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)
250 const float *src_r = src[0];
251 const float *src_g = src[1];
252 const float *src_b = src[2];
254 for (y = 0; y < h; y++) {
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]);
264 src_r += src_linesize;
265 src_g += src_linesize;
266 src_b += src_linesize;
270 static void filter_plane(AVFilterContext *ctx,
271 float *dst, int dst_linesize,
272 const float *src, int src_linesize,
276 DCTdnoizContext *s = ctx->priv;
278 const float *weights = s->weights;
281 memset(dst, 0, h * dst_linesize * sizeof(*dst));
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);
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);
296 for (by = 0; by < BSIZE; by++) {
297 for (bx = 0; bx < BSIZE; bx++) {
298 if (FFABS(*ftb) < s->th)
304 idct_block(s, dst + x, dst_linesize);
306 src += s->step * src_linesize;
307 dst += s->step * dst_linesize;
312 for (y = 0; y < h; y++) {
313 for (x = 0; x < w; x++)
314 dst[x] *= weights[x];
316 weights += dst_linesize;
320 static int filter_frame(AVFilterLink *inlink, AVFrame *in)
322 AVFilterContext *ctx = inlink->dst;
323 DCTdnoizContext *s = ctx->priv;
324 AVFilterLink *outlink = inlink->dst->outputs[0];
328 if (av_frame_is_writable(in)) {
333 out = ff_get_video_buffer(outlink, outlink->w, outlink->h);
336 return AVERROR(ENOMEM);
338 av_frame_copy_props(out, in);
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);
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);
360 uint8_t *dstp = dst + s->pr_width * 3;
361 const uint8_t *srcp = src + s->pr_width * 3;
363 for (y = 0; y < s->pr_height; y++) {
364 memcpy(dstp, srcp, hpad);
365 dstp += dst_linesize;
366 srcp += src_linesize;
370 uint8_t *dstp = dst + s->pr_height * dst_linesize;
371 const uint8_t *srcp = src + s->pr_height * src_linesize;
373 for (y = 0; y < vpad; y++) {
374 memcpy(dstp, srcp, inlink->w * 3);
375 dstp += dst_linesize;
376 srcp += src_linesize;
383 return ff_filter_frame(outlink, out);
386 static av_cold void uninit(AVFilterContext *ctx)
389 DCTdnoizContext *s = ctx->priv;
394 av_free(s->tmp_block);
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]);
401 av_expr_free(s->expr);
404 static const AVFilterPad dctdnoiz_inputs[] = {
407 .type = AVMEDIA_TYPE_VIDEO,
408 .filter_frame = filter_frame,
409 .config_props = config_input,
414 static const AVFilterPad dctdnoiz_outputs[] = {
417 .type = AVMEDIA_TYPE_VIDEO,
422 AVFilter ff_vf_dctdnoiz = {
424 .description = NULL_IF_CONFIG_SMALL("Denoise frames using 2D DCT."),
425 .priv_size = sizeof(DCTdnoizContext),
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,