git.sesse.net Git - ffmpeg/blob - 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 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/avassert.h"
  32 #include "libavutil/eval.h"
  33 #include "libavutil/opt.h"
  34 #include "internal.h"
  35
  36 static const char *const var_names[] = { "c", NULL };
  37 enum { VAR_C, VAR_VARS_NB };
  38
  39 typedef struct DCTdnoizContext {
  40     const AVClass *class;
  41
  42     /* coefficient factor expression */
  43     char *expr_str;
  44     AVExpr *expr;
  45     double var_values[VAR_VARS_NB];
  46
  47     int pr_width, pr_height;    // width and height to process
  48     float sigma;                // used when no expression are st
  49     float th;                   // threshold (3*sigma)
  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 (blocksize - overlap)
  55     int n;                      // 1<<n is the block size
  56     int bsize;                  // block size, 1<<n
  57     void (*filter_freq_func)(struct DCTdnoizContext *s,
  58                              const float *src, int src_linesize,
  59                              float *dst, int dst_linesize);
  60     void (*color_decorrelation)(float **dst, int dst_linesize,
  61                                 const uint8_t *src, int src_linesize,
  62                                 int w, int h);
  63     void (*color_correlation)(uint8_t *dst, int dst_linesize,
  64                               float **src, int src_linesize,
  65                               int w, int h);
  66 } DCTdnoizContext;
  67
  68 #define MIN_NBITS 3 /* blocksize = 1<<3 =  8 */
  69 #define MAX_NBITS 4 /* blocksize = 1<<4 = 16 */
  70 #define DEFAULT_NBITS 3
  71
  72 #define OFFSET(x) offsetof(DCTdnoizContext, x)
  73 #define FLAGS AV_OPT_FLAG_FILTERING_PARAM|AV_OPT_FLAG_VIDEO_PARAM
  74 static const AVOption dctdnoiz_options[] = {
  75     { "sigma",   "set noise sigma constant",               OFFSET(sigma),    AV_OPT_TYPE_FLOAT,  {.dbl=0},            0, 999,          .flags = FLAGS },
  76     { "s",       "set noise sigma constant",               OFFSET(sigma),    AV_OPT_TYPE_FLOAT,  {.dbl=0},            0, 999,          .flags = FLAGS },
  77     { "overlap", "set number of block overlapping pixels", OFFSET(overlap),  AV_OPT_TYPE_INT,    {.i64=-1}, -1, (1<<MAX_NBITS)-1, .flags = FLAGS },
  78     { "expr",    "set coefficient factor expression",      OFFSET(expr_str), AV_OPT_TYPE_STRING, {.str=NULL},                          .flags = FLAGS },
  79     { "e",       "set coefficient factor expression",      OFFSET(expr_str), AV_OPT_TYPE_STRING, {.str=NULL},                          .flags = FLAGS },
  80     { "n",       "set the block size, expressed in bits",  OFFSET(n),        AV_OPT_TYPE_INT,    {.i64=DEFAULT_NBITS}, MIN_NBITS, MAX_NBITS, .flags = FLAGS },
  81     { NULL }
  82 };
  83
  84 AVFILTER_DEFINE_CLASS(dctdnoiz);
  85
  86 static void av_always_inline fdct8_1d(float *dst, const float *src,
  87                                       int dst_stridea, int dst_strideb,
  88                                       int src_stridea, int src_strideb)
  89 {
  90     int i;
  91
  92     for (i = 0; i < 8; i++) {
  93         const float x00 = src[0*src_stridea] + src[7*src_stridea];
  94         const float x01 = src[1*src_stridea] + src[6*src_stridea];
  95         const float x02 = src[2*src_stridea] + src[5*src_stridea];
  96         const float x03 = src[3*src_stridea] + src[4*src_stridea];
  97         const float x04 = src[0*src_stridea] - src[7*src_stridea];
  98         const float x05 = src[1*src_stridea] - src[6*src_stridea];
  99         const float x06 = src[2*src_stridea] - src[5*src_stridea];
 100         const float x07 = src[3*src_stridea] - src[4*src_stridea];
 101         const float x08 = x00 + x03;
 102         const float x09 = x01 + x02;
 103         const float x0a = x00 - x03;
 104         const float x0b = x01 - x02;
 105         const float x0c = 1.38703984532215f*x04 + 0.275899379282943f*x07;
 106         const float x0d = 1.17587560241936f*x05 + 0.785694958387102f*x06;
 107         const float x0e = -0.785694958387102f*x05 + 1.17587560241936f*x06;
 108         const float x0f = 0.275899379282943f*x04 - 1.38703984532215f*x07;
 109         const float x10 = 0.353553390593274f * (x0c - x0d);
 110         const float x11 = 0.353553390593274f * (x0e - x0f);
 111         dst[0*dst_stridea] = 0.353553390593274f * (x08 + x09);
 112         dst[1*dst_stridea] = 0.353553390593274f * (x0c + x0d);
 113         dst[2*dst_stridea] = 0.461939766255643f*x0a + 0.191341716182545f*x0b;
 114         dst[3*dst_stridea] = 0.707106781186547f * (x10 - x11);
 115         dst[4*dst_stridea] = 0.353553390593274f * (x08 - x09);
 116         dst[5*dst_stridea] = 0.707106781186547f * (x10 + x11);
 117         dst[6*dst_stridea] = 0.191341716182545f*x0a - 0.461939766255643f*x0b;
 118         dst[7*dst_stridea] = 0.353553390593274f * (x0e + x0f);
 119         dst += dst_strideb;
 120         src += src_strideb;
 121     }
 122 }
 123
 124 static void av_always_inline idct8_1d(float *dst, const float *src,
 125                                       int dst_stridea, int dst_strideb,
 126                                       int src_stridea, int src_strideb,
 127                                       int add)
 128 {
 129     int i;
 130
 131     for (i = 0; i < 8; i++) {
 132         const float x00 =  1.4142135623731f  *src[0*src_stridea];
 133         const float x01 =  1.38703984532215f *src[1*src_stridea] + 0.275899379282943f*src[7*src_stridea];
 134         const float x02 =  1.30656296487638f *src[2*src_stridea] + 0.541196100146197f*src[6*src_stridea];
 135         const float x03 =  1.17587560241936f *src[3*src_stridea] + 0.785694958387102f*src[5*src_stridea];
 136         const float x04 =  1.4142135623731f  *src[4*src_stridea];
 137         const float x05 = -0.785694958387102f*src[3*src_stridea] + 1.17587560241936f*src[5*src_stridea];
 138         const float x06 =  0.541196100146197f*src[2*src_stridea] - 1.30656296487638f*src[6*src_stridea];
 139         const float x07 = -0.275899379282943f*src[1*src_stridea] + 1.38703984532215f*src[7*src_stridea];
 140         const float x09 = x00 + x04;
 141         const float x0a = x01 + x03;
 142         const float x0b = 1.4142135623731f*x02;
 143         const float x0c = x00 - x04;
 144         const float x0d = x01 - x03;
 145         const float x0e = 0.353553390593274f * (x09 - x0b);
 146         const float x0f = 0.353553390593274f * (x0c + x0d);
 147         const float x10 = 0.353553390593274f * (x0c - x0d);
 148         const float x11 = 1.4142135623731f*x06;
 149         const float x12 = x05 + x07;
 150         const float x13 = x05 - x07;
 151         const float x14 = 0.353553390593274f * (x11 + x12);
 152         const float x15 = 0.353553390593274f * (x11 - x12);
 153         const float x16 = 0.5f*x13;
 154         const float x08 = -x15;
 155         dst[0*dst_stridea] = (add ? dst[ 0*dst_stridea] : 0) + 0.25f * (x09 + x0b) + 0.353553390593274f*x0a;
 156         dst[1*dst_stridea] = (add ? dst[ 1*dst_stridea] : 0) + 0.707106781186547f * (x0f - x08);
 157         dst[2*dst_stridea] = (add ? dst[ 2*dst_stridea] : 0) + 0.707106781186547f * (x0f + x08);
 158         dst[3*dst_stridea] = (add ? dst[ 3*dst_stridea] : 0) + 0.707106781186547f * (x0e + x16);
 159         dst[4*dst_stridea] = (add ? dst[ 4*dst_stridea] : 0) + 0.707106781186547f * (x0e - x16);
 160         dst[5*dst_stridea] = (add ? dst[ 5*dst_stridea] : 0) + 0.707106781186547f * (x10 - x14);
 161         dst[6*dst_stridea] = (add ? dst[ 6*dst_stridea] : 0) + 0.707106781186547f * (x10 + x14);
 162         dst[7*dst_stridea] = (add ? dst[ 7*dst_stridea] : 0) + 0.25f * (x09 + x0b) - 0.353553390593274f*x0a;
 163         dst += dst_strideb;
 164         src += src_strideb;
 165     }
 166 }
 167
 168
 169 static void av_always_inline fdct16_1d(float *dst, const float *src,
 170                                        int dst_stridea, int dst_strideb,
 171                                        int src_stridea, int src_strideb)
 172 {
 173     int i;
 174
 175     for (i = 0; i < 16; i++) {
 176         const float x00 = src[ 0*src_stridea] + src[15*src_stridea];
 177         const float x01 = src[ 1*src_stridea] + src[14*src_stridea];
 178         const float x02 = src[ 2*src_stridea] + src[13*src_stridea];
 179         const float x03 = src[ 3*src_stridea] + src[12*src_stridea];
 180         const float x04 = src[ 4*src_stridea] + src[11*src_stridea];
 181         const float x05 = src[ 5*src_stridea] + src[10*src_stridea];
 182         const float x06 = src[ 6*src_stridea] + src[ 9*src_stridea];
 183         const float x07 = src[ 7*src_stridea] + src[ 8*src_stridea];
 184         const float x08 = src[ 0*src_stridea] - src[15*src_stridea];
 185         const float x09 = src[ 1*src_stridea] - src[14*src_stridea];
 186         const float x0a = src[ 2*src_stridea] - src[13*src_stridea];
 187         const float x0b = src[ 3*src_stridea] - src[12*src_stridea];
 188         const float x0c = src[ 4*src_stridea] - src[11*src_stridea];
 189         const float x0d = src[ 5*src_stridea] - src[10*src_stridea];
 190         const float x0e = src[ 6*src_stridea] - src[ 9*src_stridea];
 191         const float x0f = src[ 7*src_stridea] - src[ 8*src_stridea];
 192         const float x10 = x00 + x07;
 193         const float x11 = x01 + x06;
 194         const float x12 = x02 + x05;
 195         const float x13 = x03 + x04;
 196         const float x14 = x00 - x07;
 197         const float x15 = x01 - x06;
 198         const float x16 = x02 - x05;
 199         const float x17 = x03 - x04;
 200         const float x18 = x10 + x13;
 201         const float x19 = x11 + x12;
 202         const float x1a = x10 - x13;
 203         const float x1b = x11 - x12;
 204         const float x1c =   1.38703984532215f*x14 + 0.275899379282943f*x17;
 205         const float x1d =   1.17587560241936f*x15 + 0.785694958387102f*x16;
 206         const float x1e = -0.785694958387102f*x15 + 1.17587560241936f *x16;
 207         const float x1f =  0.275899379282943f*x14 - 1.38703984532215f *x17;
 208         const float x20 = 0.25f * (x1c - x1d);
 209         const float x21 = 0.25f * (x1e - x1f);
 210         const float x22 =  1.40740373752638f *x08 + 0.138617169199091f*x0f;
 211         const float x23 =  1.35331800117435f *x09 + 0.410524527522357f*x0e;
 212         const float x24 =  1.24722501298667f *x0a + 0.666655658477747f*x0d;
 213         const float x25 =  1.09320186700176f *x0b + 0.897167586342636f*x0c;
 214         const float x26 = -0.897167586342636f*x0b + 1.09320186700176f *x0c;
 215         const float x27 =  0.666655658477747f*x0a - 1.24722501298667f *x0d;
 216         const float x28 = -0.410524527522357f*x09 + 1.35331800117435f *x0e;
 217         const float x29 =  0.138617169199091f*x08 - 1.40740373752638f *x0f;
 218         const float x2a = x22 + x25;
 219         const float x2b = x23 + x24;
 220         const float x2c = x22 - x25;
 221         const float x2d = x23 - x24;
 222         const float x2e = 0.25f * (x2a - x2b);
 223         const float x2f = 0.326640741219094f*x2c + 0.135299025036549f*x2d;
 224         const float x30 = 0.135299025036549f*x2c - 0.326640741219094f*x2d;
 225         const float x31 = x26 + x29;
 226         const float x32 = x27 + x28;
 227         const float x33 = x26 - x29;
 228         const float x34 = x27 - x28;
 229         const float x35 = 0.25f * (x31 - x32);
 230         const float x36 = 0.326640741219094f*x33 + 0.135299025036549f*x34;
 231         const float x37 = 0.135299025036549f*x33 - 0.326640741219094f*x34;
 232         dst[ 0*dst_stridea] = 0.25f * (x18 + x19);
 233         dst[ 1*dst_stridea] = 0.25f * (x2a + x2b);
 234         dst[ 2*dst_stridea] = 0.25f * (x1c + x1d);
 235         dst[ 3*dst_stridea] = 0.707106781186547f * (x2f - x37);
 236         dst[ 4*dst_stridea] = 0.326640741219094f*x1a + 0.135299025036549f*x1b;
 237         dst[ 5*dst_stridea] = 0.707106781186547f * (x2f + x37);
 238         dst[ 6*dst_stridea] = 0.707106781186547f * (x20 - x21);
 239         dst[ 7*dst_stridea] = 0.707106781186547f * (x2e + x35);
 240         dst[ 8*dst_stridea] = 0.25f * (x18 - x19);
 241         dst[ 9*dst_stridea] = 0.707106781186547f * (x2e - x35);
 242         dst[10*dst_stridea] = 0.707106781186547f * (x20 + x21);
 243         dst[11*dst_stridea] = 0.707106781186547f * (x30 - x36);
 244         dst[12*dst_stridea] = 0.135299025036549f*x1a - 0.326640741219094f*x1b;
 245         dst[13*dst_stridea] = 0.707106781186547f * (x30 + x36);
 246         dst[14*dst_stridea] = 0.25f * (x1e + x1f);
 247         dst[15*dst_stridea] = 0.25f * (x31 + x32);
 248         dst += dst_strideb;
 249         src += src_strideb;
 250     }
 251 }
 252
 253 static void av_always_inline idct16_1d(float *dst, const float *src,
 254                                        int dst_stridea, int dst_strideb,
 255                                        int src_stridea, int src_strideb,
 256                                        int add)
 257 {
 258     int i;
 259
 260     for (i = 0; i < 16; i++) {
 261         const float x00 =  1.4142135623731f  *src[ 0*src_stridea];
 262         const float x01 =  1.40740373752638f *src[ 1*src_stridea] + 0.138617169199091f*src[15*src_stridea];
 263         const float x02 =  1.38703984532215f *src[ 2*src_stridea] + 0.275899379282943f*src[14*src_stridea];
 264         const float x03 =  1.35331800117435f *src[ 3*src_stridea] + 0.410524527522357f*src[13*src_stridea];
 265         const float x04 =  1.30656296487638f *src[ 4*src_stridea] + 0.541196100146197f*src[12*src_stridea];
 266         const float x05 =  1.24722501298667f *src[ 5*src_stridea] + 0.666655658477747f*src[11*src_stridea];
 267         const float x06 =  1.17587560241936f *src[ 6*src_stridea] + 0.785694958387102f*src[10*src_stridea];
 268         const float x07 =  1.09320186700176f *src[ 7*src_stridea] + 0.897167586342636f*src[ 9*src_stridea];
 269         const float x08 =  1.4142135623731f  *src[ 8*src_stridea];
 270         const float x09 = -0.897167586342636f*src[ 7*src_stridea] + 1.09320186700176f*src[ 9*src_stridea];
 271         const float x0a =  0.785694958387102f*src[ 6*src_stridea] - 1.17587560241936f*src[10*src_stridea];
 272         const float x0b = -0.666655658477747f*src[ 5*src_stridea] + 1.24722501298667f*src[11*src_stridea];
 273         const float x0c =  0.541196100146197f*src[ 4*src_stridea] - 1.30656296487638f*src[12*src_stridea];
 274         const float x0d = -0.410524527522357f*src[ 3*src_stridea] + 1.35331800117435f*src[13*src_stridea];
 275         const float x0e =  0.275899379282943f*src[ 2*src_stridea] - 1.38703984532215f*src[14*src_stridea];
 276         const float x0f = -0.138617169199091f*src[ 1*src_stridea] + 1.40740373752638f*src[15*src_stridea];
 277         const float x12 = x00 + x08;
 278         const float x13 = x01 + x07;
 279         const float x14 = x02 + x06;
 280         const float x15 = x03 + x05;
 281         const float x16 = 1.4142135623731f*x04;
 282         const float x17 = x00 - x08;
 283         const float x18 = x01 - x07;
 284         const float x19 = x02 - x06;
 285         const float x1a = x03 - x05;
 286         const float x1d = x12 + x16;
 287         const float x1e = x13 + x15;
 288         const float x1f = 1.4142135623731f*x14;
 289         const float x20 = x12 - x16;
 290         const float x21 = x13 - x15;
 291         const float x22 = 0.25f * (x1d - x1f);
 292         const float x23 = 0.25f * (x20 + x21);
 293         const float x24 = 0.25f * (x20 - x21);
 294         const float x25 = 1.4142135623731f*x17;
 295         const float x26 = 1.30656296487638f*x18 + 0.541196100146197f*x1a;
 296         const float x27 = 1.4142135623731f*x19;
 297         const float x28 = -0.541196100146197f*x18 + 1.30656296487638f*x1a;
 298         const float x29 = 0.176776695296637f * (x25 + x27) + 0.25f*x26;
 299         const float x2a = 0.25f * (x25 - x27);
 300         const float x2b = 0.176776695296637f * (x25 + x27) - 0.25f*x26;
 301         const float x2c = 0.353553390593274f*x28;
 302         const float x1b = 0.707106781186547f * (x2a - x2c);
 303         const float x1c = 0.707106781186547f * (x2a + x2c);
 304         const float x2d = 1.4142135623731f*x0c;
 305         const float x2e = x0b + x0d;
 306         const float x2f = x0a + x0e;
 307         const float x30 = x09 + x0f;
 308         const float x31 = x09 - x0f;
 309         const float x32 = x0a - x0e;
 310         const float x33 = x0b - x0d;
 311         const float x37 = 1.4142135623731f*x2d;
 312         const float x38 = 1.30656296487638f*x2e + 0.541196100146197f*x30;
 313         const float x39 = 1.4142135623731f*x2f;
 314         const float x3a = -0.541196100146197f*x2e + 1.30656296487638f*x30;
 315         const float x3b = 0.176776695296637f * (x37 + x39) + 0.25f*x38;
 316         const float x3c = 0.25f * (x37 - x39);
 317         const float x3d = 0.176776695296637f * (x37 + x39) - 0.25f*x38;
 318         const float x3e = 0.353553390593274f*x3a;
 319         const float x34 = 0.707106781186547f * (x3c - x3e);
 320         const float x35 = 0.707106781186547f * (x3c + x3e);
 321         const float x3f = 1.4142135623731f*x32;
 322         const float x40 = x31 + x33;
 323         const float x41 = x31 - x33;
 324         const float x42 = 0.25f * (x3f + x40);
 325         const float x43 = 0.25f * (x3f - x40);
 326         const float x44 = 0.353553390593274f*x41;
 327         const float x36 = -x43;
 328         const float x10 = -x34;
 329         const float x11 = -x3d;
 330         dst[ 0*dst_stridea] = (add ? dst[ 0*dst_stridea] : 0) + 0.176776695296637f * (x1d + x1f) + 0.25f*x1e;
 331         dst[ 1*dst_stridea] = (add ? dst[ 1*dst_stridea] : 0) + 0.707106781186547f * (x29 - x11);
 332         dst[ 2*dst_stridea] = (add ? dst[ 2*dst_stridea] : 0) + 0.707106781186547f * (x29 + x11);
 333         dst[ 3*dst_stridea] = (add ? dst[ 3*dst_stridea] : 0) + 0.707106781186547f * (x23 + x36);
 334         dst[ 4*dst_stridea] = (add ? dst[ 4*dst_stridea] : 0) + 0.707106781186547f * (x23 - x36);
 335         dst[ 5*dst_stridea] = (add ? dst[ 5*dst_stridea] : 0) + 0.707106781186547f * (x1b - x35);
 336         dst[ 6*dst_stridea] = (add ? dst[ 6*dst_stridea] : 0) + 0.707106781186547f * (x1b + x35);
 337         dst[ 7*dst_stridea] = (add ? dst[ 7*dst_stridea] : 0) + 0.707106781186547f * (x22 + x44);
 338         dst[ 8*dst_stridea] = (add ? dst[ 8*dst_stridea] : 0) + 0.707106781186547f * (x22 - x44);
 339         dst[ 9*dst_stridea] = (add ? dst[ 9*dst_stridea] : 0) + 0.707106781186547f * (x1c - x10);
 340         dst[10*dst_stridea] = (add ? dst[10*dst_stridea] : 0) + 0.707106781186547f * (x1c + x10);
 341         dst[11*dst_stridea] = (add ? dst[11*dst_stridea] : 0) + 0.707106781186547f * (x24 + x42);
 342         dst[12*dst_stridea] = (add ? dst[12*dst_stridea] : 0) + 0.707106781186547f * (x24 - x42);
 343         dst[13*dst_stridea] = (add ? dst[13*dst_stridea] : 0) + 0.707106781186547f * (x2b - x3b);
 344         dst[14*dst_stridea] = (add ? dst[14*dst_stridea] : 0) + 0.707106781186547f * (x2b + x3b);
 345         dst[15*dst_stridea] = (add ? dst[15*dst_stridea] : 0) + 0.176776695296637f * (x1d + x1f) - 0.25f*x1e;
 346         dst += dst_strideb;
 347         src += src_strideb;
 348     }
 349 }
 350
 351 #define DEF_FILTER_FREQ_FUNCS(bsize)                                                        \
 352 static av_always_inline void filter_freq_##bsize(const float *src, int src_linesize,        \
 353                                                  float *dst, int dst_linesize,              \
 354                                                  AVExpr *expr, double *var_values,          \
 355                                                  int sigma_th)                              \
 356 {                                                                                           \
 357     unsigned i;                                                                             \
 358     DECLARE_ALIGNED(32, float, tmp_block1)[bsize * bsize];                                  \
 359     DECLARE_ALIGNED(32, float, tmp_block2)[bsize * bsize];                                  \
 360                                                                                             \
 361     /* forward DCT */                                                                       \
 362     fdct##bsize##_1d(tmp_block1, src, 1, bsize, 1, src_linesize);                           \
 363     fdct##bsize##_1d(tmp_block2, tmp_block1, bsize, 1, bsize, 1);                           \
 364                                                                                             \
 365     for (i = 0; i < bsize*bsize; i++) {                                                     \
 366         float *b = &tmp_block2[i];                                                          \
 367         /* frequency filtering */                                                           \
 368         if (expr) {                                                                         \
 369             var_values[VAR_C] = FFABS(*b);                                                  \
 370             *b *= av_expr_eval(expr, var_values, NULL);                                     \
 371         } else {                                                                            \
 372             if (FFABS(*b) < sigma_th)                                                       \
 373                 *b = 0;                                                                     \
 374         }                                                                                   \
 375     }                                                                                       \
 376                                                                                             \
 377     /* inverse DCT */                                                                       \
 378     idct##bsize##_1d(tmp_block1, tmp_block2, 1, bsize, 1, bsize, 0);                        \
 379     idct##bsize##_1d(dst, tmp_block1, dst_linesize, 1, bsize, 1, 1);                        \
 380 }                                                                                           \
 381                                                                                             \
 382 static void filter_freq_sigma_##bsize(DCTdnoizContext *s,                                   \
 383                                       const float *src, int src_linesize,                   \
 384                                       float *dst, int dst_linesize)                         \
 385 {                                                                                           \
 386     filter_freq_##bsize(src, src_linesize, dst, dst_linesize, NULL, NULL, s->th);           \
 387 }                                                                                           \
 388                                                                                             \
 389 static void filter_freq_expr_##bsize(DCTdnoizContext *s,                                    \
 390                                      const float *src, int src_linesize,                    \
 391                                      float *dst, int dst_linesize)                          \
 392 {                                                                                           \
 393     filter_freq_##bsize(src, src_linesize, dst, dst_linesize, s->expr, s->var_values, 0);   \
 394 }
 395
 396 DEF_FILTER_FREQ_FUNCS(8)
 397 DEF_FILTER_FREQ_FUNCS(16)
 398
 399 #define DCT3X3_0_0  0.5773502691896258f /*  1/sqrt(3) */
 400 #define DCT3X3_0_1  0.5773502691896258f /*  1/sqrt(3) */
 401 #define DCT3X3_0_2  0.5773502691896258f /*  1/sqrt(3) */
 402 #define DCT3X3_1_0  0.7071067811865475f /*  1/sqrt(2) */
 403 #define DCT3X3_1_2 -0.7071067811865475f /* -1/sqrt(2) */
 404 #define DCT3X3_2_0  0.4082482904638631f /*  1/sqrt(6) */
 405 #define DCT3X3_2_1 -0.8164965809277261f /* -2/sqrt(6) */
 406 #define DCT3X3_2_2  0.4082482904638631f /*  1/sqrt(6) */
 407
 408 static av_always_inline void color_decorrelation(float **dst, int dst_linesize,
 409                                                  const uint8_t *src, int src_linesize,
 410                                                  int w, int h,
 411                                                  int r, int g, int b)
 412 {
 413     int x, y;
 414     float *dstp_r = dst[0];
 415     float *dstp_g = dst[1];
 416     float *dstp_b = dst[2];
 417
 418     for (y = 0; y < h; y++) {
 419         const uint8_t *srcp = src;
 420
 421         for (x = 0; x < w; x++) {
 422             dstp_r[x] = srcp[r] * DCT3X3_0_0 + srcp[g] * DCT3X3_0_1 + srcp[b] * DCT3X3_0_2;
 423             dstp_g[x] = srcp[r] * DCT3X3_1_0 +                        srcp[b] * DCT3X3_1_2;
 424             dstp_b[x] = srcp[r] * DCT3X3_2_0 + srcp[g] * DCT3X3_2_1 + srcp[b] * DCT3X3_2_2;
 425             srcp += 3;
 426         }
 427         src += src_linesize;
 428         dstp_r += dst_linesize;
 429         dstp_g += dst_linesize;
 430         dstp_b += dst_linesize;
 431     }
 432 }
 433
 434 static av_always_inline void color_correlation(uint8_t *dst, int dst_linesize,
 435                                                float **src, int src_linesize,
 436                                                int w, int h,
 437                                                int r, int g, int b)
 438 {
 439     int x, y;
 440     const float *src_r = src[0];
 441     const float *src_g = src[1];
 442     const float *src_b = src[2];
 443
 444     for (y = 0; y < h; y++) {
 445         uint8_t *dstp = dst;
 446
 447         for (x = 0; x < w; x++) {
 448             dstp[r] = av_clip_uint8(src_r[x] * DCT3X3_0_0 + src_g[x] * DCT3X3_1_0 + src_b[x] * DCT3X3_2_0);
 449             dstp[g] = av_clip_uint8(src_r[x] * DCT3X3_0_1 +                         src_b[x] * DCT3X3_2_1);
 450             dstp[b] = av_clip_uint8(src_r[x] * DCT3X3_0_2 + src_g[x] * DCT3X3_1_2 + src_b[x] * DCT3X3_2_2);
 451             dstp += 3;
 452         }
 453         dst += dst_linesize;
 454         src_r += src_linesize;
 455         src_g += src_linesize;
 456         src_b += src_linesize;
 457     }
 458 }
 459
 460 #define DECLARE_COLOR_FUNCS(name, r, g, b)                                          \
 461 static void color_decorrelation_##name(float **dst, int dst_linesize,               \
 462                                        const uint8_t *src, int src_linesize,        \
 463                                        int w, int h)                                \
 464 {                                                                                   \
 465     color_decorrelation(dst, dst_linesize, src, src_linesize, w, h, r, g, b);       \
 466 }                                                                                   \
 467                                                                                     \
 468 static void color_correlation_##name(uint8_t *dst, int dst_linesize,                \
 469                                      float **src, int src_linesize,                 \
 470                                      int w, int h)                                  \
 471 {                                                                                   \
 472     color_correlation(dst, dst_linesize, src, src_linesize, w, h, r, g, b);         \
 473 }
 474
 475 DECLARE_COLOR_FUNCS(rgb, 0, 1, 2)
 476 DECLARE_COLOR_FUNCS(bgr, 2, 1, 0)
 477
 478 static int config_input(AVFilterLink *inlink)
 479 {
 480     AVFilterContext *ctx = inlink->dst;
 481     DCTdnoizContext *s = ctx->priv;
 482     int i, x, y, bx, by, linesize, *iweights;
 483     const int bsize = 1 << s->n;
 484
 485     switch (inlink->format) {
 486     case AV_PIX_FMT_BGR24:
 487         s->color_decorrelation = color_decorrelation_bgr;
 488         s->color_correlation   = color_correlation_bgr;
 489         break;
 490     case AV_PIX_FMT_RGB24:
 491         s->color_decorrelation = color_decorrelation_rgb;
 492         s->color_correlation   = color_correlation_rgb;
 493         break;
 494     default:
 495         av_assert0(0);
 496     }
 497
 498     s->pr_width  = inlink->w - (inlink->w - bsize) % s->step;
 499     s->pr_height = inlink->h - (inlink->h - bsize) % s->step;
 500     if (s->pr_width != inlink->w)
 501         av_log(ctx, AV_LOG_WARNING, "The last %d horizontal pixels won't be denoised\n",
 502                inlink->w - s->pr_width);
 503     if (s->pr_height != inlink->h)
 504         av_log(ctx, AV_LOG_WARNING, "The last %d vertical pixels won't be denoised\n",
 505                inlink->h - s->pr_height);
 506
 507     s->p_linesize = linesize = FFALIGN(s->pr_width, 32);
 508     for (i = 0; i < 2; i++) {
 509         s->cbuf[i][0] = av_malloc(linesize * s->pr_height * sizeof(*s->cbuf[i][0]));
 510         s->cbuf[i][1] = av_malloc(linesize * s->pr_height * sizeof(*s->cbuf[i][1]));
 511         s->cbuf[i][2] = av_malloc(linesize * s->pr_height * sizeof(*s->cbuf[i][2]));
 512         if (!s->cbuf[i][0] || !s->cbuf[i][1] || !s->cbuf[i][2])
 513             return AVERROR(ENOMEM);
 514     }
 515
 516     s->weights = av_malloc(s->pr_height * linesize * sizeof(*s->weights));
 517     if (!s->weights)
 518         return AVERROR(ENOMEM);
 519     iweights = av_calloc(s->pr_height, linesize * sizeof(*iweights));
 520     if (!iweights)
 521         return AVERROR(ENOMEM);
 522     for (y = 0; y < s->pr_height - bsize + 1; y += s->step)
 523         for (x = 0; x < s->pr_width - bsize + 1; x += s->step)
 524             for (by = 0; by < bsize; by++)
 525                 for (bx = 0; bx < bsize; bx++)
 526                     iweights[(y + by)*linesize + x + bx]++;
 527     for (y = 0; y < s->pr_height; y++)
 528         for (x = 0; x < s->pr_width; x++)
 529             s->weights[y*linesize + x] = 1. / iweights[y*linesize + x];
 530     av_free(iweights);
 531
 532     return 0;
 533 }
 534
 535 static av_cold int init(AVFilterContext *ctx)
 536 {
 537     DCTdnoizContext *s = ctx->priv;
 538
 539     s->bsize = 1 << s->n;
 540     if (s->overlap == -1)
 541         s->overlap = s->bsize - 1;
 542
 543     if (s->overlap > s->bsize - 1) {
 544         av_log(s, AV_LOG_ERROR, "Overlap value can not except %d "
 545                "with a block size of %dx%d\n",
 546                s->bsize - 1, s->bsize, s->bsize);
 547         return AVERROR(EINVAL);
 548     }
 549
 550     if (s->expr_str) {
 551         int ret = av_expr_parse(&s->expr, s->expr_str, var_names,
 552                                 NULL, NULL, NULL, NULL, 0, ctx);
 553         if (ret < 0)
 554             return ret;
 555         switch (s->n) {
 556         case 3: s->filter_freq_func = filter_freq_expr_8;  break;
 557         case 4: s->filter_freq_func = filter_freq_expr_16; break;
 558         default: av_assert0(0);
 559         }
 560     } else {
 561         switch (s->n) {
 562         case 3: s->filter_freq_func = filter_freq_sigma_8;  break;
 563         case 4: s->filter_freq_func = filter_freq_sigma_16; break;
 564         default: av_assert0(0);
 565         }
 566     }
 567
 568     s->th   = s->sigma * 3.;
 569     s->step = s->bsize - s->overlap;
 570     return 0;
 571 }
 572
 573 static int query_formats(AVFilterContext *ctx)
 574 {
 575     static const enum AVPixelFormat pix_fmts[] = {
 576         AV_PIX_FMT_BGR24, AV_PIX_FMT_RGB24,
 577         AV_PIX_FMT_NONE
 578     };
 579     ff_set_common_formats(ctx, ff_make_format_list(pix_fmts));
 580     return 0;
 581 }
 582
 583 static void filter_plane(AVFilterContext *ctx,
 584                          float *dst, int dst_linesize,
 585                          const float *src, int src_linesize,
 586                          int w, int h)
 587 {
 588     int x, y;
 589     DCTdnoizContext *s = ctx->priv;
 590     float *dst0 = dst;
 591     const float *weights = s->weights;
 592
 593     // reset block sums
 594     memset(dst, 0, h * dst_linesize * sizeof(*dst));
 595
 596     // block dct sums
 597     for (y = 0; y < h - s->bsize + 1; y += s->step) {
 598         for (x = 0; x < w - s->bsize + 1; x += s->step)
 599             s->filter_freq_func(s, src + x, src_linesize,
 600                                    dst + x, dst_linesize);
 601         src += s->step * src_linesize;
 602         dst += s->step * dst_linesize;
 603     }
 604
 605     // average blocks
 606     dst = dst0;
 607     for (y = 0; y < h; y++) {
 608         for (x = 0; x < w; x++)
 609             dst[x] *= weights[x];
 610         dst += dst_linesize;
 611         weights += dst_linesize;
 612     }
 613 }
 614
 615 static int filter_frame(AVFilterLink *inlink, AVFrame *in)
 616 {
 617     AVFilterContext *ctx = inlink->dst;
 618     DCTdnoizContext *s = ctx->priv;
 619     AVFilterLink *outlink = inlink->dst->outputs[0];
 620     int direct, plane;
 621     AVFrame *out;
 622
 623     if (av_frame_is_writable(in)) {
 624         direct = 1;
 625         out = in;
 626     } else {
 627         direct = 0;
 628         out = ff_get_video_buffer(outlink, outlink->w, outlink->h);
 629         if (!out) {
 630             av_frame_free(&in);
 631             return AVERROR(ENOMEM);
 632         }
 633         av_frame_copy_props(out, in);
 634     }
 635
 636     s->color_decorrelation(s->cbuf[0], s->p_linesize,
 637                            in->data[0], in->linesize[0],
 638                            s->pr_width, s->pr_height);
 639     for (plane = 0; plane < 3; plane++)
 640         filter_plane(ctx, s->cbuf[1][plane], s->p_linesize,
 641                           s->cbuf[0][plane], s->p_linesize,
 642                           s->pr_width, s->pr_height);
 643     s->color_correlation(out->data[0], out->linesize[0],
 644                          s->cbuf[1], s->p_linesize,
 645                          s->pr_width, s->pr_height);
 646
 647     if (!direct) {
 648         int y;
 649         uint8_t *dst = out->data[0];
 650         const uint8_t *src = in->data[0];
 651         const int dst_linesize = out->linesize[0];
 652         const int src_linesize = in->linesize[0];
 653         const int hpad = (inlink->w - s->pr_width) * 3;
 654         const int vpad = (inlink->h - s->pr_height);
 655
 656         if (hpad) {
 657             uint8_t       *dstp = dst + s->pr_width * 3;
 658             const uint8_t *srcp = src + s->pr_width * 3;
 659
 660             for (y = 0; y < s->pr_height; y++) {
 661                 memcpy(dstp, srcp, hpad);
 662                 dstp += dst_linesize;
 663                 srcp += src_linesize;
 664             }
 665         }
 666         if (vpad) {
 667             uint8_t       *dstp = dst + s->pr_height * dst_linesize;
 668             const uint8_t *srcp = src + s->pr_height * src_linesize;
 669
 670             for (y = 0; y < vpad; y++) {
 671                 memcpy(dstp, srcp, inlink->w * 3);
 672                 dstp += dst_linesize;
 673                 srcp += src_linesize;
 674             }
 675         }
 676
 677         av_frame_free(&in);
 678     }
 679
 680     return ff_filter_frame(outlink, out);
 681 }
 682
 683 static av_cold void uninit(AVFilterContext *ctx)
 684 {
 685     int i;
 686     DCTdnoizContext *s = ctx->priv;
 687
 688     av_free(s->weights);
 689     for (i = 0; i < 2; i++) {
 690         av_free(s->cbuf[i][0]);
 691         av_free(s->cbuf[i][1]);
 692         av_free(s->cbuf[i][2]);
 693     }
 694     av_expr_free(s->expr);
 695 }
 696
 697 static const AVFilterPad dctdnoiz_inputs[] = {
 698     {
 699         .name         = "default",
 700         .type         = AVMEDIA_TYPE_VIDEO,
 701         .filter_frame = filter_frame,
 702         .config_props = config_input,
 703     },
 704     { NULL }
 705 };
 706
 707 static const AVFilterPad dctdnoiz_outputs[] = {
 708     {
 709         .name = "default",
 710         .type = AVMEDIA_TYPE_VIDEO,
 711     },
 712     { NULL }
 713 };
 714
 715 AVFilter ff_vf_dctdnoiz = {
 716     .name          = "dctdnoiz",
 717     .description   = NULL_IF_CONFIG_SMALL("Denoise frames using 2D DCT."),
 718     .priv_size     = sizeof(DCTdnoizContext),
 719     .init          = init,
 720     .uninit        = uninit,
 721     .query_formats = query_formats,
 722     .inputs        = dctdnoiz_inputs,
 723     .outputs       = dctdnoiz_outputs,
 724     .priv_class    = &dctdnoiz_class,
 725     .flags         = AVFILTER_FLAG_SUPPORT_TIMELINE_GENERIC,
 726 };