avfilter/dctdnoiz: rewrite [f/i]dct
[ffmpeg.git] / libavfilter / vf_dctdnoiz.c
1 /*
2  * Copyright (c) 2013-2014 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  *
24  * @see http://www.ipol.im/pub/art/2011/ys-dct/
25  *
26  * The DCT factorization used is based on "Fast and numerically stable
27  * algorithms for discrete cosine transforms" from Gerlind Plonkaa & Manfred
28  * Tasche (DOI: 10.1016/j.laa.2004.07.015).
29  */
30
31 #include "libavutil/eval.h"
32 #include "libavutil/opt.h"
33 #include "drawutils.h"
34 #include "internal.h"
35
36 #define NBITS 4
37 #define BSIZE (1<<(NBITS))
38
39 static const char *const var_names[] = { "c", NULL };
40 enum { VAR_C, VAR_VARS_NB };
41
42 typedef struct DCTdnoizContext {
43     const AVClass *class;
44
45     /* coefficient factor expression */
46     char *expr_str;
47     AVExpr *expr;
48     double var_values[VAR_VARS_NB];
49
50     int pr_width, pr_height;    // width and height to process
51     float sigma;                // used when no expression are st
52     float th;                   // threshold (3*sigma)
53     float color_dct[3][3];      // 3x3 DCT for color decorrelation
54     float *cbuf[2][3];          // two planar rgb color buffers
55     float *weights;             // dct coeff are cumulated with overlapping; these values are used for averaging
56     int p_linesize;             // line sizes for color and weights
57     int overlap;                // number of block overlapping pixels
58     int step;                   // block step increment (BSIZE - overlap)
59     void (*filter_freq_func)(struct DCTdnoizContext *s,
60                              const float *src, int src_linesize,
61                              float *dst, int dst_linesize);
62 } DCTdnoizContext;
63
64 #define OFFSET(x) offsetof(DCTdnoizContext, x)
65 #define FLAGS AV_OPT_FLAG_FILTERING_PARAM|AV_OPT_FLAG_VIDEO_PARAM
66 static const AVOption dctdnoiz_options[] = {
67     { "sigma",   "set noise sigma constant",               OFFSET(sigma),    AV_OPT_TYPE_FLOAT,  {.dbl=0},            0, 999,          .flags = FLAGS },
68     { "s",       "set noise sigma constant",               OFFSET(sigma),    AV_OPT_TYPE_FLOAT,  {.dbl=0},            0, 999,          .flags = FLAGS },
69     { "overlap", "set number of block overlapping pixels", OFFSET(overlap),  AV_OPT_TYPE_INT,    {.i64=(1<<NBITS)-1}, 0, (1<<NBITS)-1, .flags = FLAGS },
70     { "expr",    "set coefficient factor expression",      OFFSET(expr_str), AV_OPT_TYPE_STRING, {.str=NULL},                          .flags = FLAGS },
71     { "e",       "set coefficient factor expression",      OFFSET(expr_str), AV_OPT_TYPE_STRING, {.str=NULL},                          .flags = FLAGS },
72     { NULL }
73 };
74
75 AVFILTER_DEFINE_CLASS(dctdnoiz);
76
77 static void av_always_inline fdct16_1d(float *dst, const float *src,
78                                        int dst_stridea, int dst_strideb,
79                                        int src_stridea, int src_strideb)
80 {
81     int i;
82
83     for (i = 0; i < BSIZE; i++) {
84         const float x00 = src[ 0*src_stridea] + src[15*src_stridea];
85         const float x01 = src[ 1*src_stridea] + src[14*src_stridea];
86         const float x02 = src[ 2*src_stridea] + src[13*src_stridea];
87         const float x03 = src[ 3*src_stridea] + src[12*src_stridea];
88         const float x04 = src[ 4*src_stridea] + src[11*src_stridea];
89         const float x05 = src[ 5*src_stridea] + src[10*src_stridea];
90         const float x06 = src[ 6*src_stridea] + src[ 9*src_stridea];
91         const float x07 = src[ 7*src_stridea] + src[ 8*src_stridea];
92         const float x08 = src[ 0*src_stridea] - src[15*src_stridea];
93         const float x09 = src[ 1*src_stridea] - src[14*src_stridea];
94         const float x0a = src[ 2*src_stridea] - src[13*src_stridea];
95         const float x0b = src[ 3*src_stridea] - src[12*src_stridea];
96         const float x0c = src[ 4*src_stridea] - src[11*src_stridea];
97         const float x0d = src[ 5*src_stridea] - src[10*src_stridea];
98         const float x0e = src[ 6*src_stridea] - src[ 9*src_stridea];
99         const float x0f = src[ 7*src_stridea] - src[ 8*src_stridea];
100         const float x10 = x00 + x07;
101         const float x11 = x01 + x06;
102         const float x12 = x02 + x05;
103         const float x13 = x03 + x04;
104         const float x14 = x00 - x07;
105         const float x15 = x01 - x06;
106         const float x16 = x02 - x05;
107         const float x17 = x03 - x04;
108         const float x18 = x10 + x13;
109         const float x19 = x11 + x12;
110         const float x1a = x10 - x13;
111         const float x1b = x11 - x12;
112         const float x1c = 1.38703984532215*x14 + 0.275899379282943*x17;
113         const float x1d = 1.17587560241936*x15 + 0.785694958387102*x16;
114         const float x1e = -0.785694958387102*x15 + 1.17587560241936*x16;
115         const float x1f = 0.275899379282943*x14 - 1.38703984532215*x17;
116         const float x20 = 0.25 * (x1c - x1d);
117         const float x21 = 0.25 * (x1e - x1f);
118         const float x22 = 1.40740373752638*x08 + 0.138617169199091*x0f;
119         const float x23 = 1.35331800117435*x09 + 0.410524527522357*x0e;
120         const float x24 = 1.24722501298667*x0a + 0.666655658477747*x0d;
121         const float x25 = 1.09320186700176*x0b + 0.897167586342636*x0c;
122         const float x26 = -0.897167586342636*x0b + 1.09320186700176*x0c;
123         const float x27 = 0.666655658477747*x0a - 1.24722501298667*x0d;
124         const float x28 = -0.410524527522357*x09 + 1.35331800117435*x0e;
125         const float x29 = 0.138617169199091*x08 - 1.40740373752638*x0f;
126         const float x2a = x22 + x25;
127         const float x2b = x23 + x24;
128         const float x2c = x22 - x25;
129         const float x2d = x23 - x24;
130         const float x2e = 0.25 * (x2a - x2b);
131         const float x2f = 0.326640741219094*x2c + 0.135299025036549*x2d;
132         const float x30 = 0.135299025036549*x2c - 0.326640741219094*x2d;
133         const float x31 = x26 + x29;
134         const float x32 = x27 + x28;
135         const float x33 = x26 - x29;
136         const float x34 = x27 - x28;
137         const float x35 = 0.25 * (x31 - x32);
138         const float x36 = 0.326640741219094*x33 + 0.135299025036549*x34;
139         const float x37 = 0.135299025036549*x33 - 0.326640741219094*x34;
140         dst[ 0*dst_stridea] = 0.25 * (x18 + x19);
141         dst[ 1*dst_stridea] = 0.25 * (x2a + x2b);
142         dst[ 2*dst_stridea] = 0.25 * (x1c + x1d);
143         dst[ 3*dst_stridea] = 0.707106781186547 * (x2f - x37);
144         dst[ 4*dst_stridea] = 0.326640741219094*x1a + 0.135299025036549*x1b;
145         dst[ 5*dst_stridea] = 0.707106781186547 * (x2f + x37);
146         dst[ 6*dst_stridea] = 0.707106781186547 * (x20 - x21);
147         dst[ 7*dst_stridea] = 0.707106781186547 * (x2e + x35);
148         dst[ 8*dst_stridea] = 0.25 * (x18 - x19);
149         dst[ 9*dst_stridea] = 0.707106781186547 * (x2e - x35);
150         dst[10*dst_stridea] = 0.707106781186547 * (x20 + x21);
151         dst[11*dst_stridea] = 0.707106781186547 * (x30 - x36);
152         dst[12*dst_stridea] = 0.135299025036549*x1a - 0.326640741219094*x1b;
153         dst[13*dst_stridea] = 0.707106781186547 * (x30 + x36);
154         dst[14*dst_stridea] = 0.25 * (x1e + x1f);
155         dst[15*dst_stridea] = 0.25 * (x31 + x32);
156         dst += dst_strideb;
157         src += src_strideb;
158     }
159 }
160
161 static void av_always_inline idct16_1d(float *dst, const float *src,
162                                        int dst_stridea, int dst_strideb,
163                                        int src_stridea, int src_strideb,
164                                        int add)
165 {
166     int i;
167
168     for (i = 0; i < BSIZE; i++) {
169         const float x00 =  1.4142135623731  *src[ 0*src_stridea];
170         const float x01 =  1.40740373752638 *src[ 1*src_stridea] + 0.138617169199091*src[15*src_stridea];
171         const float x02 =  1.38703984532215 *src[ 2*src_stridea] + 0.275899379282943*src[14*src_stridea];
172         const float x03 =  1.35331800117435 *src[ 3*src_stridea] + 0.410524527522357*src[13*src_stridea];
173         const float x04 =  1.30656296487638 *src[ 4*src_stridea] + 0.541196100146197*src[12*src_stridea];
174         const float x05 =  1.24722501298667 *src[ 5*src_stridea] + 0.666655658477747*src[11*src_stridea];
175         const float x06 =  1.17587560241936 *src[ 6*src_stridea] + 0.785694958387102*src[10*src_stridea];
176         const float x07 =  1.09320186700176 *src[ 7*src_stridea] + 0.897167586342636*src[ 9*src_stridea];
177         const float x08 =  1.4142135623731  *src[ 8*src_stridea];
178         const float x09 = -0.897167586342636*src[ 7*src_stridea] + 1.09320186700176*src[ 9*src_stridea];
179         const float x0a =  0.785694958387102*src[ 6*src_stridea] - 1.17587560241936*src[10*src_stridea];
180         const float x0b = -0.666655658477747*src[ 5*src_stridea] + 1.24722501298667*src[11*src_stridea];
181         const float x0c =  0.541196100146197*src[ 4*src_stridea] - 1.30656296487638*src[12*src_stridea];
182         const float x0d = -0.410524527522357*src[ 3*src_stridea] + 1.35331800117435*src[13*src_stridea];
183         const float x0e =  0.275899379282943*src[ 2*src_stridea] - 1.38703984532215*src[14*src_stridea];
184         const float x0f = -0.138617169199091*src[ 1*src_stridea] + 1.40740373752638*src[15*src_stridea];
185         const float x12 = x00 + x08;
186         const float x13 = x01 + x07;
187         const float x14 = x02 + x06;
188         const float x15 = x03 + x05;
189         const float x16 = 1.4142135623731*x04;
190         const float x17 = x00 - x08;
191         const float x18 = x01 - x07;
192         const float x19 = x02 - x06;
193         const float x1a = x03 - x05;
194         const float x1d = x12 + x16;
195         const float x1e = x13 + x15;
196         const float x1f = 1.4142135623731*x14;
197         const float x20 = x12 - x16;
198         const float x21 = x13 - x15;
199         const float x22 = 0.25 * (x1d - x1f);
200         const float x23 = 0.25 * (x20 + x21);
201         const float x24 = 0.25 * (x20 - x21);
202         const float x25 = 1.4142135623731*x17;
203         const float x26 = 1.30656296487638*x18 + 0.541196100146197*x1a;
204         const float x27 = 1.4142135623731*x19;
205         const float x28 = -0.541196100146197*x18 + 1.30656296487638*x1a;
206         const float x29 = 0.176776695296637 * (x25 + x27) + 0.25*x26;
207         const float x2a = 0.25 * (x25 - x27);
208         const float x2b = 0.176776695296637 * (x25 + x27) - 0.25*x26;
209         const float x2c = 0.353553390593274*x28;
210         const float x1b = 0.707106781186547 * (x2a - x2c);
211         const float x1c = 0.707106781186547 * (x2a + x2c);
212         const float x2d = 1.4142135623731*x0c;
213         const float x2e = x0b + x0d;
214         const float x2f = x0a + x0e;
215         const float x30 = x09 + x0f;
216         const float x31 = x09 - x0f;
217         const float x32 = x0a - x0e;
218         const float x33 = x0b - x0d;
219         const float x37 = 1.4142135623731*x2d;
220         const float x38 = 1.30656296487638*x2e + 0.541196100146197*x30;
221         const float x39 = 1.4142135623731*x2f;
222         const float x3a = -0.541196100146197*x2e + 1.30656296487638*x30;
223         const float x3b = 0.176776695296637 * (x37 + x39) + 0.25*x38;
224         const float x3c = 0.25 * (x37 - x39);
225         const float x3d = 0.176776695296637 * (x37 + x39) - 0.25*x38;
226         const float x3e = 0.353553390593274*x3a;
227         const float x34 = 0.707106781186547 * (x3c - x3e);
228         const float x35 = 0.707106781186547 * (x3c + x3e);
229         const float x3f = 1.4142135623731*x32;
230         const float x40 = x31 + x33;
231         const float x41 = x31 - x33;
232         const float x42 = 0.25 * (x3f + x40);
233         const float x43 = 0.25 * (x3f - x40);
234         const float x44 = 0.353553390593274*x41;
235         const float x36 = -x43;
236         const float x10 = -x34;
237         const float x11 = -x3d;
238         dst[ 0*dst_stridea] = (add ? dst[ 0*dst_stridea] : 0) + 0.176776695296637 * (x1d + x1f) + 0.25*x1e;
239         dst[ 1*dst_stridea] = (add ? dst[ 1*dst_stridea] : 0) + 0.707106781186547 * (x29 - x11);
240         dst[ 2*dst_stridea] = (add ? dst[ 2*dst_stridea] : 0) + 0.707106781186547 * (x29 + x11);
241         dst[ 3*dst_stridea] = (add ? dst[ 3*dst_stridea] : 0) + 0.707106781186547 * (x23 + x36);
242         dst[ 4*dst_stridea] = (add ? dst[ 4*dst_stridea] : 0) + 0.707106781186547 * (x23 - x36);
243         dst[ 5*dst_stridea] = (add ? dst[ 5*dst_stridea] : 0) + 0.707106781186547 * (x1b - x35);
244         dst[ 6*dst_stridea] = (add ? dst[ 6*dst_stridea] : 0) + 0.707106781186547 * (x1b + x35);
245         dst[ 7*dst_stridea] = (add ? dst[ 7*dst_stridea] : 0) + 0.707106781186547 * (x22 + x44);
246         dst[ 8*dst_stridea] = (add ? dst[ 8*dst_stridea] : 0) + 0.707106781186547 * (x22 - x44);
247         dst[ 9*dst_stridea] = (add ? dst[ 9*dst_stridea] : 0) + 0.707106781186547 * (x1c - x10);
248         dst[10*dst_stridea] = (add ? dst[10*dst_stridea] : 0) + 0.707106781186547 * (x1c + x10);
249         dst[11*dst_stridea] = (add ? dst[11*dst_stridea] : 0) + 0.707106781186547 * (x24 + x42);
250         dst[12*dst_stridea] = (add ? dst[12*dst_stridea] : 0) + 0.707106781186547 * (x24 - x42);
251         dst[13*dst_stridea] = (add ? dst[13*dst_stridea] : 0) + 0.707106781186547 * (x2b - x3b);
252         dst[14*dst_stridea] = (add ? dst[14*dst_stridea] : 0) + 0.707106781186547 * (x2b + x3b);
253         dst[15*dst_stridea] = (add ? dst[15*dst_stridea] : 0) + 0.176776695296637 * (x1d + x1f) - 0.25*x1e;
254         dst += dst_strideb;
255         src += src_strideb;
256     }
257 }
258
259 static av_always_inline void filter_freq(const float *src, int src_linesize,
260                                          float *dst, int dst_linesize,
261                                          AVExpr *expr, double *var_values,
262                                          int sigma_th)
263 {
264     unsigned i;
265     DECLARE_ALIGNED(32, float, tmp_block1)[BSIZE * BSIZE];
266     DECLARE_ALIGNED(32, float, tmp_block2)[BSIZE * BSIZE];
267
268     /* forward DCT */
269     fdct16_1d(tmp_block1, src, 1, BSIZE, 1, src_linesize);
270     fdct16_1d(tmp_block2, tmp_block1, BSIZE, 1, BSIZE, 1);
271
272     for (i = 0; i < BSIZE*BSIZE; i++) {
273         float *b = &tmp_block2[i];
274         /* frequency filtering */
275         if (expr) {
276             var_values[VAR_C] = FFABS(*b);
277             *b *= av_expr_eval(expr, var_values, NULL);
278         } else {
279             if (FFABS(*b) < sigma_th)
280                 *b = 0;
281         }
282     }
283
284     /* inverse DCT */
285     idct16_1d(tmp_block1, tmp_block2, 1, BSIZE, 1, BSIZE, 0);
286     idct16_1d(dst, tmp_block1, dst_linesize, 1, BSIZE, 1, 1);
287 }
288
289 static void filter_freq_sigma(DCTdnoizContext *s,
290                               const float *src, int src_linesize,
291                               float *dst, int dst_linesize)
292 {
293     filter_freq(src, src_linesize, dst, dst_linesize, NULL, NULL, s->th);
294 }
295
296 static void filter_freq_expr(DCTdnoizContext *s,
297                              const float *src, int src_linesize,
298                              float *dst, int dst_linesize)
299 {
300     filter_freq(src, src_linesize, dst, dst_linesize, s->expr, s->var_values, 0);
301 }
302
303 static int config_input(AVFilterLink *inlink)
304 {
305     AVFilterContext *ctx = inlink->dst;
306     DCTdnoizContext *s = ctx->priv;
307     int i, x, y, bx, by, linesize, *iweights;
308     const float dct_3x3[3][3] = {
309         { 1./sqrt(3),  1./sqrt(3),  1./sqrt(3) },
310         { 1./sqrt(2),           0, -1./sqrt(2) },
311         { 1./sqrt(6), -2./sqrt(6),  1./sqrt(6) },
312     };
313     uint8_t rgba_map[4];
314
315     ff_fill_rgba_map(rgba_map, inlink->format);
316     for (y = 0; y < 3; y++)
317         for (x = 0; x < 3; x++)
318             s->color_dct[y][x] = dct_3x3[rgba_map[y]][rgba_map[x]];
319
320     s->pr_width  = inlink->w - (inlink->w - BSIZE) % s->step;
321     s->pr_height = inlink->h - (inlink->h - BSIZE) % s->step;
322     if (s->pr_width != inlink->w)
323         av_log(ctx, AV_LOG_WARNING, "The last %d horizontal pixels won't be denoised\n",
324                inlink->w - s->pr_width);
325     if (s->pr_height != inlink->h)
326         av_log(ctx, AV_LOG_WARNING, "The last %d vertical pixels won't be denoised\n",
327                inlink->h - s->pr_height);
328
329     s->p_linesize = linesize = FFALIGN(s->pr_width, 32);
330     for (i = 0; i < 2; i++) {
331         s->cbuf[i][0] = av_malloc(linesize * s->pr_height * sizeof(*s->cbuf[i][0]));
332         s->cbuf[i][1] = av_malloc(linesize * s->pr_height * sizeof(*s->cbuf[i][1]));
333         s->cbuf[i][2] = av_malloc(linesize * s->pr_height * sizeof(*s->cbuf[i][2]));
334         if (!s->cbuf[i][0] || !s->cbuf[i][1] || !s->cbuf[i][2])
335             return AVERROR(ENOMEM);
336     }
337
338     s->weights = av_malloc(s->pr_height * linesize * sizeof(*s->weights));
339     if (!s->weights)
340         return AVERROR(ENOMEM);
341     iweights = av_calloc(s->pr_height, linesize * sizeof(*iweights));
342     if (!iweights)
343         return AVERROR(ENOMEM);
344     for (y = 0; y < s->pr_height - BSIZE + 1; y += s->step)
345         for (x = 0; x < s->pr_width - BSIZE + 1; x += s->step)
346             for (by = 0; by < BSIZE; by++)
347                 for (bx = 0; bx < BSIZE; bx++)
348                     iweights[(y + by)*linesize + x + bx]++;
349     for (y = 0; y < s->pr_height; y++)
350         for (x = 0; x < s->pr_width; x++)
351             s->weights[y*linesize + x] = 1. / iweights[y*linesize + x];
352     av_free(iweights);
353
354     return 0;
355 }
356
357 static av_cold int init(AVFilterContext *ctx)
358 {
359     DCTdnoizContext *s = ctx->priv;
360
361     if (s->expr_str) {
362         int ret = av_expr_parse(&s->expr, s->expr_str, var_names,
363                                 NULL, NULL, NULL, NULL, 0, ctx);
364         if (ret < 0)
365             return ret;
366         s->filter_freq_func = filter_freq_expr;
367     } else {
368         s->filter_freq_func = filter_freq_sigma;
369     }
370
371     s->th   = s->sigma * 3.;
372     s->step = BSIZE - s->overlap;
373     return 0;
374 }
375
376 static int query_formats(AVFilterContext *ctx)
377 {
378     static const enum AVPixelFormat pix_fmts[] = {
379         AV_PIX_FMT_BGR24, AV_PIX_FMT_RGB24,
380         AV_PIX_FMT_NONE
381     };
382     ff_set_common_formats(ctx, ff_make_format_list(pix_fmts));
383     return 0;
384 }
385
386 static void color_decorrelation(float dct3ch[3][3], float **dst, int dst_linesize,
387                                 const uint8_t *src, int src_linesize, int w, int h)
388 {
389     int x, y;
390     float *dstp_r = dst[0];
391     float *dstp_g = dst[1];
392     float *dstp_b = dst[2];
393
394     for (y = 0; y < h; y++) {
395         const uint8_t *srcp = src;
396
397         for (x = 0; x < w; x++) {
398             dstp_r[x] = srcp[0] * dct3ch[0][0] + srcp[1] * dct3ch[0][1] + srcp[2] * dct3ch[0][2];
399             dstp_g[x] = srcp[0] * dct3ch[1][0] + srcp[1] * dct3ch[1][1] + srcp[2] * dct3ch[1][2];
400             dstp_b[x] = srcp[0] * dct3ch[2][0] + srcp[1] * dct3ch[2][1] + srcp[2] * dct3ch[2][2];
401             srcp += 3;
402         }
403         src += src_linesize;
404         dstp_r += dst_linesize;
405         dstp_g += dst_linesize;
406         dstp_b += dst_linesize;
407     }
408 }
409
410 static void color_correlation(float dct3ch[3][3], uint8_t *dst, int dst_linesize,
411                               float **src, int src_linesize, int w, int h)
412 {
413     int x, y;
414     const float *src_r = src[0];
415     const float *src_g = src[1];
416     const float *src_b = src[2];
417
418     for (y = 0; y < h; y++) {
419         uint8_t *dstp = dst;
420
421         for (x = 0; x < w; x++) {
422             dstp[0] = av_clip_uint8(src_r[x] * dct3ch[0][0] + src_g[x] * dct3ch[1][0] + src_b[x] * dct3ch[2][0]);
423             dstp[1] = av_clip_uint8(src_r[x] * dct3ch[0][1] + src_g[x] * dct3ch[1][1] + src_b[x] * dct3ch[2][1]);
424             dstp[2] = av_clip_uint8(src_r[x] * dct3ch[0][2] + src_g[x] * dct3ch[1][2] + src_b[x] * dct3ch[2][2]);
425             dstp += 3;
426         }
427         dst += dst_linesize;
428         src_r += src_linesize;
429         src_g += src_linesize;
430         src_b += src_linesize;
431     }
432 }
433
434 static void filter_plane(AVFilterContext *ctx,
435                          float *dst, int dst_linesize,
436                          const float *src, int src_linesize,
437                          int w, int h)
438 {
439     int x, y;
440     DCTdnoizContext *s = ctx->priv;
441     float *dst0 = dst;
442     const float *weights = s->weights;
443
444     // reset block sums
445     memset(dst, 0, h * dst_linesize * sizeof(*dst));
446
447     // block dct sums
448     for (y = 0; y < h - BSIZE + 1; y += s->step) {
449         for (x = 0; x < w - BSIZE + 1; x += s->step)
450             s->filter_freq_func(s, src + x, src_linesize,
451                                    dst + x, dst_linesize);
452         src += s->step * src_linesize;
453         dst += s->step * dst_linesize;
454     }
455
456     // average blocks
457     dst = dst0;
458     for (y = 0; y < h; y++) {
459         for (x = 0; x < w; x++)
460             dst[x] *= weights[x];
461         dst += dst_linesize;
462         weights += dst_linesize;
463     }
464 }
465
466 static int filter_frame(AVFilterLink *inlink, AVFrame *in)
467 {
468     AVFilterContext *ctx = inlink->dst;
469     DCTdnoizContext *s = ctx->priv;
470     AVFilterLink *outlink = inlink->dst->outputs[0];
471     int direct, plane;
472     AVFrame *out;
473
474     if (av_frame_is_writable(in)) {
475         direct = 1;
476         out = in;
477     } else {
478         direct = 0;
479         out = ff_get_video_buffer(outlink, outlink->w, outlink->h);
480         if (!out) {
481             av_frame_free(&in);
482             return AVERROR(ENOMEM);
483         }
484         av_frame_copy_props(out, in);
485     }
486
487     color_decorrelation(s->color_dct, s->cbuf[0], s->p_linesize,
488                         in->data[0], in->linesize[0], s->pr_width, s->pr_height);
489     for (plane = 0; plane < 3; plane++)
490         filter_plane(ctx, s->cbuf[1][plane], s->p_linesize,
491                           s->cbuf[0][plane], s->p_linesize,
492                           s->pr_width, s->pr_height);
493     color_correlation(s->color_dct, out->data[0], out->linesize[0],
494                       s->cbuf[1], s->p_linesize, s->pr_width, s->pr_height);
495
496     if (!direct) {
497         int y;
498         uint8_t *dst = out->data[0];
499         const uint8_t *src = in->data[0];
500         const int dst_linesize = out->linesize[0];
501         const int src_linesize = in->linesize[0];
502         const int hpad = (inlink->w - s->pr_width) * 3;
503         const int vpad = (inlink->h - s->pr_height);
504
505         if (hpad) {
506             uint8_t       *dstp = dst + s->pr_width * 3;
507             const uint8_t *srcp = src + s->pr_width * 3;
508
509             for (y = 0; y < s->pr_height; y++) {
510                 memcpy(dstp, srcp, hpad);
511                 dstp += dst_linesize;
512                 srcp += src_linesize;
513             }
514         }
515         if (vpad) {
516             uint8_t       *dstp = dst + s->pr_height * dst_linesize;
517             const uint8_t *srcp = src + s->pr_height * src_linesize;
518
519             for (y = 0; y < vpad; y++) {
520                 memcpy(dstp, srcp, inlink->w * 3);
521                 dstp += dst_linesize;
522                 srcp += src_linesize;
523             }
524         }
525
526         av_frame_free(&in);
527     }
528
529     return ff_filter_frame(outlink, out);
530 }
531
532 static av_cold void uninit(AVFilterContext *ctx)
533 {
534     int i;
535     DCTdnoizContext *s = ctx->priv;
536
537     av_free(s->weights);
538     for (i = 0; i < 2; i++) {
539         av_free(s->cbuf[i][0]);
540         av_free(s->cbuf[i][1]);
541         av_free(s->cbuf[i][2]);
542     }
543     av_expr_free(s->expr);
544 }
545
546 static const AVFilterPad dctdnoiz_inputs[] = {
547     {
548         .name         = "default",
549         .type         = AVMEDIA_TYPE_VIDEO,
550         .filter_frame = filter_frame,
551         .config_props = config_input,
552     },
553     { NULL }
554 };
555
556 static const AVFilterPad dctdnoiz_outputs[] = {
557     {
558         .name = "default",
559         .type = AVMEDIA_TYPE_VIDEO,
560     },
561     { NULL }
562 };
563
564 AVFilter ff_vf_dctdnoiz = {
565     .name          = "dctdnoiz",
566     .description   = NULL_IF_CONFIG_SMALL("Denoise frames using 2D DCT."),
567     .priv_size     = sizeof(DCTdnoizContext),
568     .init          = init,
569     .uninit        = uninit,
570     .query_formats = query_formats,
571     .inputs        = dctdnoiz_inputs,
572     .outputs       = dctdnoiz_outputs,
573     .priv_class    = &dctdnoiz_class,
574     .flags         = AVFILTER_FLAG_SUPPORT_TIMELINE_GENERIC,
575 };