avfilter/dctdnoiz: rewrite [f/i]dct

This removes the avcodec dependency and make the code almost twice as
fast. More to come.

The DCT factorization is based on "Fast and numerically stable
algorithms for discrete cosine transforms" from Gerlind Plonkaa &
Manfred Tasche (DOI: 10.1016/j.laa.2004.07.015).

configure

@@ -2530,8 +2530,6 @@ boxblur_filter_deps="gpl"
 bs2b_filter_deps="libbs2b"
 colormatrix_filter_deps="gpl"
 cropdetect_filter_deps="gpl"
-dctdnoiz_filter_deps="avcodec"
-dctdnoiz_filter_select="dct"
 delogo_filter_deps="gpl"
 deshake_filter_select="pixelutils"
 drawtext_filter_deps="libfreetype"

libavfilter/vf_dctdnoiz.c

 /*
- * Copyright (c) 2013 Clément Bœsch
+ * Copyright (c) 2013-2014 Clément Bœsch
  *
  * This file is part of FFmpeg.
  *
@@ -20,10 +20,14 @@
 
 /**
  * A simple, relatively efficient and extremely slow DCT image denoiser.
+ *
  * @see http://www.ipol.im/pub/art/2011/ys-dct/
+ *
+ * The DCT factorization used is based on "Fast and numerically stable
+ * algorithms for discrete cosine transforms" from Gerlind Plonkaa & Manfred
+ * Tasche (DOI: 10.1016/j.laa.2004.07.015).
  */
 
-#include "libavcodec/avfft.h"
 #include "libavutil/eval.h"
 #include "libavutil/opt.h"
 #include "drawutils.h"
@@ -35,7 +39,7 @@
 static const char *const var_names[] = { "c", NULL };
 enum { VAR_C, VAR_VARS_NB };
 
-typedef struct {
+typedef struct DCTdnoizContext {
     const AVClass *class;
 
     /* coefficient factor expression */
@@ -52,8 +56,9 @@ typedef struct {
     int p_linesize;             // line sizes for color and weights
     int overlap;                // number of block overlapping pixels
     int step;                   // block step increment (BSIZE - overlap)
-    DCTContext *dct, *idct;     // DCT and inverse DCT contexts
-    float *block, *tmp_block;   // two BSIZE x BSIZE block buffers
+    void (*filter_freq_func)(struct DCTdnoizContext *s,
+                             const float *src, int src_linesize,
+                             float *dst, int dst_linesize);
 } DCTdnoizContext;
 
 #define OFFSET(x) offsetof(DCTdnoizContext, x)
@@ -69,66 +74,230 @@ static const AVOption dctdnoiz_options[] = {
 
 AVFILTER_DEFINE_CLASS(dctdnoiz);
 
-static float *dct_block(DCTdnoizContext *ctx, const float *src, int src_linesize)
+static void av_always_inline fdct16_1d(float *dst, const float *src,
+                                       int dst_stridea, int dst_strideb,
+                                       int src_stridea, int src_strideb)
 {
-    int x, y;
-    float *column;
-
-    for (y = 0; y < BSIZE; y++) {
-        float *line = ctx->block;
+    int i;
 
-        memcpy(line, src, BSIZE * sizeof(*line));
-        src += src_linesize;
-        av_dct_calc(ctx->dct, line);
-
-        column = ctx->tmp_block + y;
-        column[0] = line[0] * (1. / sqrt(BSIZE));
-        column += BSIZE;
-        for (x = 1; x < BSIZE; x++) {
-            *column = line[x] * sqrt(2. / BSIZE);
-            column += BSIZE;
-        }
+    for (i = 0; i < BSIZE; i++) {
+        const float x00 = src[ 0*src_stridea] + src[15*src_stridea];
+        const float x01 = src[ 1*src_stridea] + src[14*src_stridea];
+        const float x02 = src[ 2*src_stridea] + src[13*src_stridea];
+        const float x03 = src[ 3*src_stridea] + src[12*src_stridea];
+        const float x04 = src[ 4*src_stridea] + src[11*src_stridea];
+        const float x05 = src[ 5*src_stridea] + src[10*src_stridea];
+        const float x06 = src[ 6*src_stridea] + src[ 9*src_stridea];
+        const float x07 = src[ 7*src_stridea] + src[ 8*src_stridea];
+        const float x08 = src[ 0*src_stridea] - src[15*src_stridea];
+        const float x09 = src[ 1*src_stridea] - src[14*src_stridea];
+        const float x0a = src[ 2*src_stridea] - src[13*src_stridea];
+        const float x0b = src[ 3*src_stridea] - src[12*src_stridea];
+        const float x0c = src[ 4*src_stridea] - src[11*src_stridea];
+        const float x0d = src[ 5*src_stridea] - src[10*src_stridea];
+        const float x0e = src[ 6*src_stridea] - src[ 9*src_stridea];
+        const float x0f = src[ 7*src_stridea] - src[ 8*src_stridea];
+        const float x10 = x00 + x07;
+        const float x11 = x01 + x06;
+        const float x12 = x02 + x05;
+        const float x13 = x03 + x04;
+        const float x14 = x00 - x07;
+        const float x15 = x01 - x06;
+        const float x16 = x02 - x05;
+        const float x17 = x03 - x04;
+        const float x18 = x10 + x13;
+        const float x19 = x11 + x12;
+        const float x1a = x10 - x13;
+        const float x1b = x11 - x12;
+        const float x1c = 1.38703984532215*x14 + 0.275899379282943*x17;
+        const float x1d = 1.17587560241936*x15 + 0.785694958387102*x16;
+        const float x1e = -0.785694958387102*x15 + 1.17587560241936*x16;
+        const float x1f = 0.275899379282943*x14 - 1.38703984532215*x17;
+        const float x20 = 0.25 * (x1c - x1d);
+        const float x21 = 0.25 * (x1e - x1f);
+        const float x22 = 1.40740373752638*x08 + 0.138617169199091*x0f;
+        const float x23 = 1.35331800117435*x09 + 0.410524527522357*x0e;
+        const float x24 = 1.24722501298667*x0a + 0.666655658477747*x0d;
+        const float x25 = 1.09320186700176*x0b + 0.897167586342636*x0c;
+        const float x26 = -0.897167586342636*x0b + 1.09320186700176*x0c;
+        const float x27 = 0.666655658477747*x0a - 1.24722501298667*x0d;
+        const float x28 = -0.410524527522357*x09 + 1.35331800117435*x0e;
+        const float x29 = 0.138617169199091*x08 - 1.40740373752638*x0f;
+        const float x2a = x22 + x25;
+        const float x2b = x23 + x24;
+        const float x2c = x22 - x25;
+        const float x2d = x23 - x24;
+        const float x2e = 0.25 * (x2a - x2b);
+        const float x2f = 0.326640741219094*x2c + 0.135299025036549*x2d;
+        const float x30 = 0.135299025036549*x2c - 0.326640741219094*x2d;
+        const float x31 = x26 + x29;
+        const float x32 = x27 + x28;
+        const float x33 = x26 - x29;
+        const float x34 = x27 - x28;
+        const float x35 = 0.25 * (x31 - x32);
+        const float x36 = 0.326640741219094*x33 + 0.135299025036549*x34;
+        const float x37 = 0.135299025036549*x33 - 0.326640741219094*x34;
+        dst[ 0*dst_stridea] = 0.25 * (x18 + x19);
+        dst[ 1*dst_stridea] = 0.25 * (x2a + x2b);
+        dst[ 2*dst_stridea] = 0.25 * (x1c + x1d);
+        dst[ 3*dst_stridea] = 0.707106781186547 * (x2f - x37);
+        dst[ 4*dst_stridea] = 0.326640741219094*x1a + 0.135299025036549*x1b;
+        dst[ 5*dst_stridea] = 0.707106781186547 * (x2f + x37);
+        dst[ 6*dst_stridea] = 0.707106781186547 * (x20 - x21);
+        dst[ 7*dst_stridea] = 0.707106781186547 * (x2e + x35);
+        dst[ 8*dst_stridea] = 0.25 * (x18 - x19);
+        dst[ 9*dst_stridea] = 0.707106781186547 * (x2e - x35);
+        dst[10*dst_stridea] = 0.707106781186547 * (x20 + x21);
+        dst[11*dst_stridea] = 0.707106781186547 * (x30 - x36);
+        dst[12*dst_stridea] = 0.135299025036549*x1a - 0.326640741219094*x1b;
+        dst[13*dst_stridea] = 0.707106781186547 * (x30 + x36);
+        dst[14*dst_stridea] = 0.25 * (x1e + x1f);
+        dst[15*dst_stridea] = 0.25 * (x31 + x32);
+        dst += dst_strideb;
+        src += src_strideb;
     }
+}
+
+static void av_always_inline idct16_1d(float *dst, const float *src,
+                                       int dst_stridea, int dst_strideb,
+                                       int src_stridea, int src_strideb,
+                                       int add)
+{
+    int i;
 
-    column = ctx->tmp_block;
-    for (x = 0; x < BSIZE; x++) {
-        av_dct_calc(ctx->dct, column);
-        column[0] *= 1. / sqrt(BSIZE);
-        for (y = 1; y < BSIZE; y++)
-            column[y] *= sqrt(2. / BSIZE);
-        column += BSIZE;
+    for (i = 0; i < BSIZE; i++) {
+        const float x00 =  1.4142135623731  *src[ 0*src_stridea];
+        const float x01 =  1.40740373752638 *src[ 1*src_stridea] + 0.138617169199091*src[15*src_stridea];
+        const float x02 =  1.38703984532215 *src[ 2*src_stridea] + 0.275899379282943*src[14*src_stridea];
+        const float x03 =  1.35331800117435 *src[ 3*src_stridea] + 0.410524527522357*src[13*src_stridea];
+        const float x04 =  1.30656296487638 *src[ 4*src_stridea] + 0.541196100146197*src[12*src_stridea];
+        const float x05 =  1.24722501298667 *src[ 5*src_stridea] + 0.666655658477747*src[11*src_stridea];
+        const float x06 =  1.17587560241936 *src[ 6*src_stridea] + 0.785694958387102*src[10*src_stridea];
+        const float x07 =  1.09320186700176 *src[ 7*src_stridea] + 0.897167586342636*src[ 9*src_stridea];
+        const float x08 =  1.4142135623731  *src[ 8*src_stridea];
+        const float x09 = -0.897167586342636*src[ 7*src_stridea] + 1.09320186700176*src[ 9*src_stridea];
+        const float x0a =  0.785694958387102*src[ 6*src_stridea] - 1.17587560241936*src[10*src_stridea];
+        const float x0b = -0.666655658477747*src[ 5*src_stridea] + 1.24722501298667*src[11*src_stridea];
+        const float x0c =  0.541196100146197*src[ 4*src_stridea] - 1.30656296487638*src[12*src_stridea];
+        const float x0d = -0.410524527522357*src[ 3*src_stridea] + 1.35331800117435*src[13*src_stridea];
+        const float x0e =  0.275899379282943*src[ 2*src_stridea] - 1.38703984532215*src[14*src_stridea];
+        const float x0f = -0.138617169199091*src[ 1*src_stridea] + 1.40740373752638*src[15*src_stridea];
+        const float x12 = x00 + x08;
+        const float x13 = x01 + x07;
+        const float x14 = x02 + x06;
+        const float x15 = x03 + x05;
+        const float x16 = 1.4142135623731*x04;
+        const float x17 = x00 - x08;
+        const float x18 = x01 - x07;
+        const float x19 = x02 - x06;
+        const float x1a = x03 - x05;
+        const float x1d = x12 + x16;
+        const float x1e = x13 + x15;
+        const float x1f = 1.4142135623731*x14;
+        const float x20 = x12 - x16;
+        const float x21 = x13 - x15;
+        const float x22 = 0.25 * (x1d - x1f);
+        const float x23 = 0.25 * (x20 + x21);
+        const float x24 = 0.25 * (x20 - x21);
+        const float x25 = 1.4142135623731*x17;
+        const float x26 = 1.30656296487638*x18 + 0.541196100146197*x1a;
+        const float x27 = 1.4142135623731*x19;
+        const float x28 = -0.541196100146197*x18 + 1.30656296487638*x1a;
+        const float x29 = 0.176776695296637 * (x25 + x27) + 0.25*x26;
+        const float x2a = 0.25 * (x25 - x27);
+        const float x2b = 0.176776695296637 * (x25 + x27) - 0.25*x26;
+        const float x2c = 0.353553390593274*x28;
+        const float x1b = 0.707106781186547 * (x2a - x2c);
+        const float x1c = 0.707106781186547 * (x2a + x2c);
+        const float x2d = 1.4142135623731*x0c;
+        const float x2e = x0b + x0d;
+        const float x2f = x0a + x0e;
+        const float x30 = x09 + x0f;
+        const float x31 = x09 - x0f;
+        const float x32 = x0a - x0e;
+        const float x33 = x0b - x0d;
+        const float x37 = 1.4142135623731*x2d;
+        const float x38 = 1.30656296487638*x2e + 0.541196100146197*x30;
+        const float x39 = 1.4142135623731*x2f;
+        const float x3a = -0.541196100146197*x2e + 1.30656296487638*x30;
+        const float x3b = 0.176776695296637 * (x37 + x39) + 0.25*x38;
+        const float x3c = 0.25 * (x37 - x39);
+        const float x3d = 0.176776695296637 * (x37 + x39) - 0.25*x38;
+        const float x3e = 0.353553390593274*x3a;
+        const float x34 = 0.707106781186547 * (x3c - x3e);
+        const float x35 = 0.707106781186547 * (x3c + x3e);
+        const float x3f = 1.4142135623731*x32;
+        const float x40 = x31 + x33;
+        const float x41 = x31 - x33;
+        const float x42 = 0.25 * (x3f + x40);
+        const float x43 = 0.25 * (x3f - x40);
+        const float x44 = 0.353553390593274*x41;
+        const float x36 = -x43;
+        const float x10 = -x34;
+        const float x11 = -x3d;
+        dst[ 0*dst_stridea] = (add ? dst[ 0*dst_stridea] : 0) + 0.176776695296637 * (x1d + x1f) + 0.25*x1e;
+        dst[ 1*dst_stridea] = (add ? dst[ 1*dst_stridea] : 0) + 0.707106781186547 * (x29 - x11);
+        dst[ 2*dst_stridea] = (add ? dst[ 2*dst_stridea] : 0) + 0.707106781186547 * (x29 + x11);
+        dst[ 3*dst_stridea] = (add ? dst[ 3*dst_stridea] : 0) + 0.707106781186547 * (x23 + x36);
+        dst[ 4*dst_stridea] = (add ? dst[ 4*dst_stridea] : 0) + 0.707106781186547 * (x23 - x36);
+        dst[ 5*dst_stridea] = (add ? dst[ 5*dst_stridea] : 0) + 0.707106781186547 * (x1b - x35);
+        dst[ 6*dst_stridea] = (add ? dst[ 6*dst_stridea] : 0) + 0.707106781186547 * (x1b + x35);
+        dst[ 7*dst_stridea] = (add ? dst[ 7*dst_stridea] : 0) + 0.707106781186547 * (x22 + x44);
+        dst[ 8*dst_stridea] = (add ? dst[ 8*dst_stridea] : 0) + 0.707106781186547 * (x22 - x44);
+        dst[ 9*dst_stridea] = (add ? dst[ 9*dst_stridea] : 0) + 0.707106781186547 * (x1c - x10);
+        dst[10*dst_stridea] = (add ? dst[10*dst_stridea] : 0) + 0.707106781186547 * (x1c + x10);
+        dst[11*dst_stridea] = (add ? dst[11*dst_stridea] : 0) + 0.707106781186547 * (x24 + x42);
+        dst[12*dst_stridea] = (add ? dst[12*dst_stridea] : 0) + 0.707106781186547 * (x24 - x42);
+        dst[13*dst_stridea] = (add ? dst[13*dst_stridea] : 0) + 0.707106781186547 * (x2b - x3b);
+        dst[14*dst_stridea] = (add ? dst[14*dst_stridea] : 0) + 0.707106781186547 * (x2b + x3b);
+        dst[15*dst_stridea] = (add ? dst[15*dst_stridea] : 0) + 0.176776695296637 * (x1d + x1f) - 0.25*x1e;
+        dst += dst_strideb;
+        src += src_strideb;
     }
+}
 
-    for (y = 0; y < BSIZE; y++)
-        for (x = 0; x < BSIZE; x++)
-            ctx->block[y*BSIZE + x] = ctx->tmp_block[x*BSIZE + y];
+static av_always_inline void filter_freq(const float *src, int src_linesize,
+                                         float *dst, int dst_linesize,
+                                         AVExpr *expr, double *var_values,
+                                         int sigma_th)
+{
+    unsigned i;
+    DECLARE_ALIGNED(32, float, tmp_block1)[BSIZE * BSIZE];
+    DECLARE_ALIGNED(32, float, tmp_block2)[BSIZE * BSIZE];
+
+    /* forward DCT */
+    fdct16_1d(tmp_block1, src, 1, BSIZE, 1, src_linesize);
+    fdct16_1d(tmp_block2, tmp_block1, BSIZE, 1, BSIZE, 1);
+
+    for (i = 0; i < BSIZE*BSIZE; i++) {
+        float *b = &tmp_block2[i];
+        /* frequency filtering */
+        if (expr) {
+            var_values[VAR_C] = FFABS(*b);
+            *b *= av_expr_eval(expr, var_values, NULL);
+        } else {
+            if (FFABS(*b) < sigma_th)
+                *b = 0;
+        }
+    }
 
-    return ctx->block;
+    /* inverse DCT */
+    idct16_1d(tmp_block1, tmp_block2, 1, BSIZE, 1, BSIZE, 0);
+    idct16_1d(dst, tmp_block1, dst_linesize, 1, BSIZE, 1, 1);
 }
 
-static void idct_block(DCTdnoizContext *ctx, float *dst, int dst_linesize)
+static void filter_freq_sigma(DCTdnoizContext *s,
+                              const float *src, int src_linesize,
+                              float *dst, int dst_linesize)
 {
-    int x, y;
-    float *block = ctx->block;
-    float *tmp = ctx->tmp_block;
-
-    for (y = 0; y < BSIZE; y++) {
-        block[0] *= sqrt(BSIZE);
-        for (x = 1; x < BSIZE; x++)
-            block[x] *= 1./sqrt(2. / BSIZE);
-        av_dct_calc(ctx->idct, block);
-        block += BSIZE;
-    }
+    filter_freq(src, src_linesize, dst, dst_linesize, NULL, NULL, s->th);
+}
 
-    block = ctx->block;
-    for (y = 0; y < BSIZE; y++) {
-        tmp[0] = block[y] * sqrt(BSIZE);
-        for (x = 1; x < BSIZE; x++)
-            tmp[x] = block[x*BSIZE + y] * (1./sqrt(2. / BSIZE));
-        av_dct_calc(ctx->idct, tmp);
-        for (x = 0; x < BSIZE; x++)
-            dst[x*dst_linesize + y] += tmp[x];
-    }
+static void filter_freq_expr(DCTdnoizContext *s,
+                             const float *src, int src_linesize,
+                             float *dst, int dst_linesize)
+{
+    filter_freq(src, src_linesize, dst, dst_linesize, s->expr, s->var_values, 0);
 }
 
 static int config_input(AVFilterLink *inlink)
@@ -194,18 +363,13 @@ static av_cold int init(AVFilterContext *ctx)
                                 NULL, NULL, NULL, NULL, 0, ctx);
         if (ret < 0)
             return ret;
+        s->filter_freq_func = filter_freq_expr;
+    } else {
+        s->filter_freq_func = filter_freq_sigma;
     }
 
     s->th   = s->sigma * 3.;
     s->step = BSIZE - s->overlap;
-    s->dct  = av_dct_init(NBITS, DCT_II);
-    s->idct = av_dct_init(NBITS, DCT_III);
-    s->block     = av_malloc(BSIZE * BSIZE * sizeof(*s->block));
-    s->tmp_block = av_malloc(BSIZE * BSIZE * sizeof(*s->tmp_block));
-
-    if (!s->dct || !s->idct || !s->tmp_block || !s->block)
-        return AVERROR(ENOMEM);
-
     return 0;
 }
 
@@ -272,7 +436,7 @@ static void filter_plane(AVFilterContext *ctx,
                          const float *src, int src_linesize,
                          int w, int h)
 {
-    int x, y, bx, by;
+    int x, y;
     DCTdnoizContext *s = ctx->priv;
     float *dst0 = dst;
     const float *weights = s->weights;
@@ -282,27 +446,9 @@ static void filter_plane(AVFilterContext *ctx,
 
     // block dct sums
     for (y = 0; y < h - BSIZE + 1; y += s->step) {
-        for (x = 0; x < w - BSIZE + 1; x += s->step) {
-            float *ftb = dct_block(s, src + x, src_linesize);
-
-            if (s->expr) {
-                for (by = 0; by < BSIZE; by++) {
-                    for (bx = 0; bx < BSIZE; bx++) {
-                        s->var_values[VAR_C] = FFABS(*ftb);
-                        *ftb++ *= av_expr_eval(s->expr, s->var_values, s);
-                    }
-                }
-            } else {
-                for (by = 0; by < BSIZE; by++) {
-                    for (bx = 0; bx < BSIZE; bx++) {
-                        if (FFABS(*ftb) < s->th)
-                            *ftb = 0;
-                        ftb++;
-                    }
-                }
-            }
-            idct_block(s, dst + x, dst_linesize);
-        }
+        for (x = 0; x < w - BSIZE + 1; x += s->step)
+            s->filter_freq_func(s, src + x, src_linesize,
+                                   dst + x, dst_linesize);
         src += s->step * src_linesize;
         dst += s->step * dst_linesize;
     }
@@ -388,10 +534,6 @@ static av_cold void uninit(AVFilterContext *ctx)
     int i;
     DCTdnoizContext *s = ctx->priv;
 
-    av_dct_end(s->dct);
-    av_dct_end(s->idct);
-    av_free(s->block);
-    av_free(s->tmp_block);
     av_free(s->weights);
     for (i = 0; i < 2; i++) {
         av_free(s->cbuf[i][0]);