avfilter: add vaguedenoiser filter
[ffmpeg.git] / libavfilter / vf_vaguedenoiser.c
1 /*
2  * Copyright (c) 2003 LeFunGus, lefungus@altern.org
3  *
4  * This file is part of FFmpeg
5  *
6  * FFmpeg is free software; you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation; either version 2 of the License, or
9  * (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
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License along
17  * with FFmpeg; if not, write to the Free Software Foundation, Inc.,
18  * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
19  */
20
21 #include <float.h>
22
23 #include "libavutil/imgutils.h"
24 #include "libavutil/attributes.h"
25 #include "libavutil/common.h"
26 #include "libavutil/pixdesc.h"
27 #include "libavutil/intreadwrite.h"
28 #include "libavutil/opt.h"
29
30 #include "avfilter.h"
31 #include "formats.h"
32 #include "internal.h"
33 #include "video.h"
34
35 typedef struct VagueDenoiserContext {
36     const AVClass *class;
37
38     float threshold;
39     float percent;
40     int method;
41     int nsteps;
42     int planes;
43
44     int depth;
45     int peak;
46     int nb_planes;
47     int planeheight[4];
48     int planewidth[4];
49
50     float *block;
51     float *in;
52     float *out;
53     float *tmp;
54
55     int hlowsize[32];
56     int hhighsize[32];
57     int vlowsize[32];
58     int vhighsize[32];
59 } VagueDenoiserContext;
60
61 #define OFFSET(x) offsetof(VagueDenoiserContext, x)
62 #define FLAGS AV_OPT_FLAG_VIDEO_PARAM | AV_OPT_FLAG_FILTERING_PARAM
63 static const AVOption vaguedenoiser_options[] = {
64     { "threshold", "set filtering strength",   OFFSET(threshold), AV_OPT_TYPE_FLOAT, {.dbl=2.},  0,DBL_MAX, FLAGS },
65     { "method",    "set filtering method",     OFFSET(method),    AV_OPT_TYPE_INT,   {.i64=2 },  0, 2,      FLAGS, "method" },
66         { "hard",   "hard thresholding",       0,                 AV_OPT_TYPE_CONST, {.i64=0},   0, 0,      FLAGS, "method" },
67         { "soft",   "soft thresholding",       0,                 AV_OPT_TYPE_CONST, {.i64=1},   0, 0,      FLAGS, "method" },
68         { "garrote", "garotte thresholding",   0,                 AV_OPT_TYPE_CONST, {.i64=2},   0, 0,      FLAGS, "method" },
69     { "nsteps",    "set number of steps",      OFFSET(nsteps),    AV_OPT_TYPE_INT,   {.i64=6 },  1, 32,     FLAGS },
70     { "percent", "set percent of full denoising", OFFSET(percent),AV_OPT_TYPE_FLOAT, {.dbl=85},  0,100,     FLAGS },
71     { "planes",    "set planes to filter",     OFFSET(planes),    AV_OPT_TYPE_INT,   {.i64=15 }, 0, 15,     FLAGS },
72     { NULL }
73 };
74
75 AVFILTER_DEFINE_CLASS(vaguedenoiser);
76
77 #define NPAD 10
78
79 static const float analysis_low[9] = {
80     0.037828455506995f, -0.023849465019380f, -0.110624404418423f, 0.377402855612654f,
81     0.852698679009403f, 0.377402855612654f, -0.110624404418423f, -0.023849465019380f, 0.037828455506995f
82 };
83
84 static const float analysis_high[7] = {
85     -0.064538882628938f, 0.040689417609558f, 0.418092273222212f, -0.788485616405664f,
86     0.418092273222212f, 0.040689417609558f, -0.064538882628938f
87 };
88
89 static const float synthesis_low[7] = {
90     -0.064538882628938f, -0.040689417609558f, 0.418092273222212f, 0.788485616405664f,
91     0.418092273222212f, -0.040689417609558f, -0.064538882628938f
92 };
93
94 static const float synthesis_high[9] = {
95     -0.037828455506995f, -0.023849465019380f, 0.110624404418423f, 0.377402855612654f,
96     -0.852698679009403f, 0.377402855612654f, 0.110624404418423f, -0.023849465019380f, -0.037828455506995f
97 };
98
99 static int query_formats(AVFilterContext *ctx)
100 {
101     static const enum AVPixelFormat pix_fmts[] = {
102         AV_PIX_FMT_GRAY8,
103         AV_PIX_FMT_GRAY16,
104         AV_PIX_FMT_YUV410P, AV_PIX_FMT_YUV411P,
105         AV_PIX_FMT_YUV420P, AV_PIX_FMT_YUV422P,
106         AV_PIX_FMT_YUV440P, AV_PIX_FMT_YUV444P,
107         AV_PIX_FMT_YUVJ420P, AV_PIX_FMT_YUVJ422P,
108         AV_PIX_FMT_YUVJ440P, AV_PIX_FMT_YUVJ444P,
109         AV_PIX_FMT_YUVJ411P,
110         AV_PIX_FMT_YUV420P9, AV_PIX_FMT_YUV422P9, AV_PIX_FMT_YUV444P9,
111         AV_PIX_FMT_YUV420P10, AV_PIX_FMT_YUV422P10, AV_PIX_FMT_YUV444P10,
112         AV_PIX_FMT_YUV440P10,
113         AV_PIX_FMT_YUV444P12, AV_PIX_FMT_YUV422P12, AV_PIX_FMT_YUV420P12,
114         AV_PIX_FMT_YUV440P12,
115         AV_PIX_FMT_YUV444P14, AV_PIX_FMT_YUV422P14, AV_PIX_FMT_YUV420P14,
116         AV_PIX_FMT_YUV420P16, AV_PIX_FMT_YUV422P16, AV_PIX_FMT_YUV444P16,
117         AV_PIX_FMT_GBRP, AV_PIX_FMT_GBRP9, AV_PIX_FMT_GBRP10,
118         AV_PIX_FMT_GBRP12, AV_PIX_FMT_GBRP14, AV_PIX_FMT_GBRP16,
119         AV_PIX_FMT_NONE
120     };
121     AVFilterFormats *fmts_list = ff_make_format_list(pix_fmts);
122     if (!fmts_list)
123         return AVERROR(ENOMEM);
124     return ff_set_common_formats(ctx, fmts_list);
125 }
126
127 static int config_input(AVFilterLink *inlink)
128 {
129     VagueDenoiserContext *s = inlink->dst->priv;
130     const AVPixFmtDescriptor *desc = av_pix_fmt_desc_get(inlink->format);
131     int nsteps_width, nsteps_height, nsteps_max;
132
133     s->depth = desc->comp[0].depth;
134     s->nb_planes = desc->nb_components;
135
136     s->planeheight[1] = s->planeheight[2] = AV_CEIL_RSHIFT(inlink->h, desc->log2_chroma_h);
137     s->planeheight[0] = s->planeheight[3] = inlink->h;
138     s->planewidth[1]  = s->planewidth[2]  = AV_CEIL_RSHIFT(inlink->w, desc->log2_chroma_w);
139     s->planewidth[0]  = s->planewidth[3]  = inlink->w;
140
141     s->block = av_malloc_array(inlink->w * inlink->h, sizeof(*s->block));
142     s->in    = av_malloc_array(32 + FFMAX(inlink->w, inlink->h), sizeof(*s->in));
143     s->out   = av_malloc_array(32 + FFMAX(inlink->w, inlink->h), sizeof(*s->out));
144     s->tmp   = av_malloc_array(32 + FFMAX(inlink->w, inlink->h), sizeof(*s->tmp));
145
146     if (!s->block || !s->in || !s->out || !s->tmp)
147         return AVERROR(ENOMEM);
148
149     s->threshold *= 1 << (s->depth - 8);
150     s->peak = (1 << s->depth) - 1;
151
152     nsteps_width  = ((s->planes & 2 || s->planes & 4) && s->nb_planes > 1) ? s->planewidth[1] : s->planewidth[0];
153     nsteps_height = ((s->planes & 2 || s->planes & 4) && s->nb_planes > 1) ? s->planeheight[1] : s->planeheight[0];
154
155     for (nsteps_max = 1; nsteps_max < 15; nsteps_max++) {
156         if (pow(2, nsteps_max) >= nsteps_width || pow(2, nsteps_max) >= nsteps_height)
157             break;
158     }
159
160     s->nsteps = FFMIN(s->nsteps, nsteps_max - 2);
161
162     return 0;
163 }
164
165 static inline void copy(const float *p1, float *p2, const int length)
166 {
167     memcpy(p2, p1, length * sizeof(float));
168 }
169
170 static inline void copyv(const float *p1, const int stride1, float *p2, const int length)
171 {
172     int i;
173
174     for (i = 0; i < length; i++) {
175         p2[i] = *p1;
176         p1 += stride1;
177     }
178 }
179
180 static inline void copyh(const float *p1, float *p2, const int stride2, const int length)
181 {
182     int i;
183
184     for (i = 0; i < length; i++) {
185         *p2 = p1[i];
186         p2 += stride2;
187     }
188 }
189
190 // Do symmetric extension of data using prescribed symmetries
191 // Original values are in output[npad] through output[npad+size-1]
192 // New values will be placed in output[0] through output[npad] and in output[npad+size] through output[2*npad+size-1] (note: end values may not be filled in)
193 // extension at left bdry is ... 3 2 1 0 | 0 1 2 3 ...
194 // same for right boundary
195 // if right_ext=1 then ... 3 2 1 0 | 1 2 3
196 static void symmetric_extension(float *output, const int size, const int left_ext, const int right_ext)
197 {
198     int first = NPAD;
199     int last = NPAD - 1 + size;
200     const int originalLast = last;
201     int i, nextend, idx;
202
203     if (left_ext == 2)
204         output[--first] = output[NPAD];
205     if (right_ext == 2)
206         output[++last] = output[originalLast];
207
208     // extend left end
209     nextend = first;
210     for (i = 0; i < nextend; i++)
211         output[--first] = output[NPAD + 1 + i];
212
213     idx = NPAD + NPAD - 1 + size;
214
215     // extend right end
216     nextend = idx - last;
217     for (i = 0; i < nextend; i++)
218         output[++last] = output[originalLast - 1 - i];
219 }
220
221 static void transform_step(float *input, float *output, const int size, const int low_size, VagueDenoiserContext *s)
222 {
223     int i;
224
225     symmetric_extension(input, size, 1, 1);
226
227     for (i = NPAD; i < NPAD + low_size; i++) {
228         const float a = input[2 * i - 14] * analysis_low[0];
229         const float b = input[2 * i - 13] * analysis_low[1];
230         const float c = input[2 * i - 12] * analysis_low[2];
231         const float d = input[2 * i - 11] * analysis_low[3];
232         const float e = input[2 * i - 10] * analysis_low[4];
233         const float f = input[2 * i -  9] * analysis_low[3];
234         const float g = input[2 * i -  8] * analysis_low[2];
235         const float h = input[2 * i -  7] * analysis_low[1];
236         const float k = input[2 * i -  6] * analysis_low[0];
237
238         output[i] = a + b + c + d + e + f + g + h + k;
239     }
240
241     for (i = NPAD; i < NPAD + low_size; i++) {
242         const float a = input[2 * i - 12] * analysis_high[0];
243         const float b = input[2 * i - 11] * analysis_high[1];
244         const float c = input[2 * i - 10] * analysis_high[2];
245         const float d = input[2 * i -  9] * analysis_high[3];
246         const float e = input[2 * i -  8] * analysis_high[2];
247         const float f = input[2 * i -  7] * analysis_high[1];
248         const float g = input[2 * i -  6] * analysis_high[0];
249
250         output[i + low_size] = a + b + c + d + e + f + g;
251     }
252 }
253
254 static void invert_step(const float *input, float *output, float *temp, const int size, VagueDenoiserContext *s)
255 {
256     const int low_size = (size + 1) >> 1;
257     const int high_size = size >> 1;
258     int left_ext = 1, right_ext, i;
259     int findex;
260
261     memcpy(temp + NPAD, input + NPAD, low_size * sizeof(float));
262
263     right_ext = (size % 2 == 0) ? 2 : 1;
264     symmetric_extension(temp, low_size, left_ext, right_ext);
265
266     memset(output, 0, (NPAD + NPAD + size) * sizeof(float));
267     findex = (size + 2) >> 1;
268
269     for (i = 9; i < findex + 11; i++) {
270         const float a = temp[i] * synthesis_low[0];
271         const float b = temp[i] * synthesis_low[1];
272         const float c = temp[i] * synthesis_low[2];
273         const float d = temp[i] * synthesis_low[3];
274
275         output[2 * i - 13] += a;
276         output[2 * i - 12] += b;
277         output[2 * i - 11] += c;
278         output[2 * i - 10] += d;
279         output[2 * i -  9] += c;
280         output[2 * i -  8] += b;
281         output[2 * i -  7] += a;
282     }
283
284     memcpy(temp + NPAD, input + NPAD + low_size, high_size * sizeof(float));
285
286     left_ext = 2;
287     right_ext = (size % 2 == 0) ? 1 : 2;
288     symmetric_extension(temp, high_size, left_ext, right_ext);
289
290     for (i = 8; i < findex + 11; i++) {
291         const float a = temp[i] * synthesis_high[0];
292         const float b = temp[i] * synthesis_high[1];
293         const float c = temp[i] * synthesis_high[2];
294         const float d = temp[i] * synthesis_high[3];
295         const float e = temp[i] * synthesis_high[4];
296
297         output[2 * i - 13] += a;
298         output[2 * i - 12] += b;
299         output[2 * i - 11] += c;
300         output[2 * i - 10] += d;
301         output[2 * i -  9] += e;
302         output[2 * i -  8] += d;
303         output[2 * i -  7] += c;
304         output[2 * i -  6] += b;
305         output[2 * i -  5] += a;
306     }
307 }
308
309 static void hard_thresholding(float *block, const int width, const int height,
310                               const int stride, const float threshold,
311                               const float percent)
312 {
313     const float frac = 1.f - percent * 0.01f;
314     int y, x;
315
316     for (y = 0; y < height; y++) {
317         for (x = 0; x < width; x++) {
318             if (FFABS(block[x]) <= threshold)
319                 block[x] *= frac;
320         }
321         block += stride;
322     }
323 }
324
325 static void soft_thresholding(float *block, const int width, const int height, const int stride,
326                               const float threshold, const float percent, const int nsteps)
327 {
328     const float frac = 1.f - percent * 0.01f;
329     const float shift = threshold * 0.01f * percent;
330     int w = width;
331     int h = height;
332     int y, x, l;
333
334     for (l = 0; l < nsteps; l++) {
335         w = (w + 1) >> 1;
336         h = (h + 1) >> 1;
337     }
338
339     for (y = 0; y < height; y++) {
340         const int x0 = (y < h) ? w : 0;
341         for (x = x0; x < width; x++) {
342             const float temp = FFABS(block[x]);
343             if (temp <= threshold)
344                 block[x] *= frac;
345             else
346                 block[x] = (block[x] < 0.f ? -1.f : (block[x] > 0.f ? 1.f : 0.f)) * (temp - shift);
347         }
348         block += stride;
349     }
350 }
351
352 static void qian_thresholding(float *block, const int width, const int height,
353                               const int stride, const float threshold,
354                               const float percent)
355 {
356     const float percent01 = percent * 0.01f;
357     const float tr2 = threshold * threshold * percent01;
358     const float frac = 1.f - percent01;
359     int y, x;
360
361     for (y = 0; y < height; y++) {
362         for (x = 0; x < width; x++) {
363             const float temp = FFABS(block[x]);
364             if (temp <= threshold) {
365                 block[x] *= frac;
366             } else {
367                 const float tp2 = temp * temp;
368                 block[x] *= (tp2 - tr2) / tp2;
369             }
370         }
371         block += stride;
372     }
373 }
374
375 static void filter(VagueDenoiserContext *s, AVFrame *in, AVFrame *out)
376 {
377     int p, y, x, i, j;
378
379     for (p = 0; p < s->nb_planes; p++) {
380         const int height = s->planeheight[p];
381         const int width = s->planewidth[p];
382         const uint8_t *srcp8 = in->data[p];
383         const uint16_t *srcp16 = (const uint16_t *)in->data[p];
384         uint8_t *dstp8 = out->data[p];
385         uint16_t *dstp16 = (uint16_t *)out->data[p];
386         float *output = s->block;
387         int h_low_size0 = width;
388         int v_low_size0 = height;
389         int nsteps_transform = s->nsteps;
390         int nsteps_invert = s->nsteps;
391         const float *input = s->block;
392
393         if (!((1 << p) & s->planes)) {
394             av_image_copy_plane(out->data[p], out->linesize[p], in->data[p], in->linesize[p],
395                                 s->planewidth[p], s->planeheight[p]);
396             continue;
397         }
398
399         if (s->depth <= 8) {
400             for (y = 0; y < height; y++) {
401                 for (x = 0; x < width; x++)
402                     output[x] = srcp8[x];
403                 srcp8 += in->linesize[p];
404                 output += width;
405             }
406         } else {
407             for (y = 0; y < height; y++) {
408                 for (x = 0; x < width; x++)
409                     output[x] = srcp16[x];
410                 srcp16 += in->linesize[p] / 2;
411                 output += width;
412             }
413         }
414
415         while (nsteps_transform--) {
416             int low_size = (h_low_size0 + 1) >> 1;
417             float *input = s->block;
418             for (j = 0; j < v_low_size0; j++) {
419                 copy(input, s->in + NPAD, h_low_size0);
420                 transform_step(s->in, s->out, h_low_size0, low_size, s);
421                 copy(s->out + NPAD, input, h_low_size0);
422                 input += width;
423             }
424
425             low_size = (v_low_size0 + 1) >> 1;
426             input = s->block;
427             for (j = 0; j < h_low_size0; j++) {
428                 copyv(input, width, s->in + NPAD, v_low_size0);
429                 transform_step(s->in, s->out, v_low_size0, low_size, s);
430                 copyh(s->out + NPAD, input, width, v_low_size0);
431                 input++;
432             }
433
434             h_low_size0 = (h_low_size0 + 1) >> 1;
435             v_low_size0 = (v_low_size0 + 1) >> 1;
436         }
437
438         if (s->method == 0)
439             hard_thresholding(s->block, width, height, width, s->threshold, s->percent);
440         else if (s->method == 1)
441             soft_thresholding(s->block, width, height, width, s->threshold, s->percent, s->nsteps);
442         else
443             qian_thresholding(s->block, width, height, width, s->threshold, s->percent);
444
445         s->hlowsize[0]  = (width + 1) >> 1;
446         s->hhighsize[0] = width >> 1;
447         s->vlowsize[0]  = (height + 1) >> 1;
448         s->vhighsize[0] = height >> 1;
449
450         for (i = 1; i < s->nsteps; i++) {
451             s->hlowsize[i]  = (s->hlowsize[i - 1] + 1) >> 1;
452             s->hhighsize[i] = s->hlowsize[i - 1] >> 1;
453             s->vlowsize[i]  = (s->vlowsize[i - 1] + 1) >> 1;
454             s->vhighsize[i] = s->vlowsize[i - 1] >> 1;
455         }
456
457         while (nsteps_invert--) {
458             const int idx = s->vlowsize[nsteps_invert] + s->vhighsize[nsteps_invert];
459             const int idx2 = s->hlowsize[nsteps_invert] + s->hhighsize[nsteps_invert];
460             float * idx3 = s->block;
461             for (i = 0; i < idx2; i++) {
462                 copyv(idx3, width, s->in + NPAD, idx);
463                 invert_step(s->in, s->out, s->tmp, idx, s);
464                 copyh(s->out + NPAD, idx3, width, idx);
465                 idx3++;
466             }
467
468             idx3 = s->block;
469             for (i = 0; i < idx; i++) {
470                 copy(idx3, s->in + NPAD, idx2);
471                 invert_step(s->in, s->out, s->tmp, idx2, s);
472                 copy(s->out + NPAD, idx3, idx2);
473                 idx3 += width;
474             }
475         }
476
477         if (s->depth <= 8) {
478             for (y = 0; y < height; y++) {
479                 for (x = 0; x < width; x++)
480                     dstp8[x] = av_clip_uint8(input[x] + 0.5f);
481                 input += width;
482                 dstp8 += out->linesize[p];
483             }
484         } else {
485             for (y = 0; y < height; y++) {
486                 for (x = 0; x < width; x++)
487                     dstp16[x] = av_clip(input[x] + 0.5f, 0, s->peak);
488                 input += width;
489                 dstp16 += out->linesize[p] / 2;
490             }
491         }
492     }
493 }
494
495 static int filter_frame(AVFilterLink *inlink, AVFrame *in)
496 {
497     AVFilterContext *ctx  = inlink->dst;
498     VagueDenoiserContext *s = ctx->priv;
499     AVFilterLink *outlink = ctx->outputs[0];
500     AVFrame *out;
501     int direct = av_frame_is_writable(in);
502
503     if (direct) {
504         out = in;
505     } else {
506         out = ff_get_video_buffer(outlink, outlink->w, outlink->h);
507         if (!out) {
508             av_frame_free(&in);
509             return AVERROR(ENOMEM);
510         }
511
512         av_frame_copy_props(out, in);
513     }
514
515     filter(s, in, out);
516
517     if (!direct)
518         av_frame_free(&in);
519
520     return ff_filter_frame(outlink, out);
521 }
522
523 static av_cold void uninit(AVFilterContext *ctx)
524 {
525     VagueDenoiserContext *s = ctx->priv;
526
527     av_freep(&s->block);
528     av_freep(&s->in);
529     av_freep(&s->out);
530     av_freep(&s->tmp);
531 }
532
533 static const AVFilterPad vaguedenoiser_inputs[] = {
534     {
535         .name         = "default",
536         .type         = AVMEDIA_TYPE_VIDEO,
537         .config_props = config_input,
538         .filter_frame = filter_frame,
539     },
540     { NULL }
541 };
542
543
544 static const AVFilterPad vaguedenoiser_outputs[] = {
545     {
546         .name = "default",
547         .type = AVMEDIA_TYPE_VIDEO
548     },
549     { NULL }
550 };
551
552 AVFilter ff_vf_vaguedenoiser = {
553     .name          = "vaguedenoiser",
554     .description   = NULL_IF_CONFIG_SMALL("Apply a Wavelet based Denoiser."),
555     .priv_size     = sizeof(VagueDenoiserContext),
556     .priv_class    = &vaguedenoiser_class,
557     .uninit        = uninit,
558     .query_formats = query_formats,
559     .inputs        = vaguedenoiser_inputs,
560     .outputs       = vaguedenoiser_outputs,
561     .flags         = AVFILTER_FLAG_SUPPORT_TIMELINE_GENERIC,
562 };