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