filter_kernels: replace AGG-based code

This commit replaces code based on AGG, taken from this source file:

http://vector-agg.cvs.sourceforge.net/viewvc/vector-agg/agg-2.5/include/agg_image_filters.h

The intention is that filter_kernels.c can be relicensed to LGPL or BSD.
Because the AGG author died, full replacement is the only way to achieve
it.

This affects only some filter functions. These are exclusively
mathematical functions for computing filter coefficients. (Other parts
in filter_kernel.c were originally written by me, with heavy additions
and refactoring done by other mpv contributors.) While the code is
mostly just well-known mathematical formulas written down in C form,
AGG copyright could perhaps be claimed anyway.

To remove the AGG code, I replaced it with the filter functions from:

https://github.com/glumpy/glumpy/blob/master/glumpy/library/build-spatial-filters.py

These functions conveniently compute exactly the same thing in mpv,
Glumpy, AGG (and about anything that will filter images using the same
mathematical principles).

First I ported the Python code in the file to C. Then I replaced all
functions in filter_kernels.c with this code that could be replaced.
Then I investigated whether the remaining functions were based on AGG
code and took appropriate action:

hanning(), hamming(), quadric(), bicubic(), kaiser(), blackman(),
spline16(), spline36(), gaussian(), sinc() were taken straight from
Glumpy.

For sinc(), re-add the "fabs(x) < 1e-8" check, which was added in commit
586dc557 for unknown reasons.

gaussian() loses its filter parameter for some reason. (Well, who cares,
not my problem.)

The really awkward thing is that the text for hanning() and hamming()
does not change. In theory these functions are now based on Glumpy code,
but it seems like this can be neither proven nor denied. (The same
happened in some other cases with at least a few lines of code.)

sphinx() was added in commit 586dc557, and looks suspiciously like
sinc() as well. Replace the first 3 lines of the body with the ported
function (of which 2 lines do not change; the first uses code only in
mpv, and the second is just "return 1.0;"). The 4th line is only similar
on an abstract level (and that because of the mathematical relation
between these functions). Although the original sinc() was probably used
as template for it, with the other lines replaced, I don't think you
could make the claim that it falls under AGG copyright.

jinc() was added in commit 26baf5b9, but the code for it might be based
on sinc(). Rewrite it based on the "new" sinc(). Some of the same
remarks as with sphinx() apply.

cubic_bc() was ported from Glumpy's Mitchell(). (As far as I'm aware,
with the default parameters it's called "the" Mitchell-Netravali filter,
but in mpv this function is used to generate a whole group of filters.)

spline64() was added in commit a8b67c66, and was probably derived from
spline36(). Re-derive it from the "new" spline36().

triangle() could be considered derived from the original bilinear().
This is this in the original commit:

    static double bilinear(kernel *k, double x)
    {
        return 1.0 - x;
    }

This _might_ be based on AGG's image_filter_bilinear:

    struct image_filter_bilinear
    {
        static double radius() { return 1.0; }
        static double calc_weight(double x)
        {
            return 1.0 - x;
        }
    };

Considering that the "framework" was written by me, and the only part
from AGG taken is "return 1.0 - x;", and this part is trivial and was
later thoroughly replaced, this is probably not under the AGG copyright.

I'm hoping this doesn't introduce regressions. But the main focus is not
being productive anyway, and I didn't rigorously check unintended
changes in functionality.

1 files changed, 69 insertions, 71 deletions

diff --git a/video/out/filter_kernels.c b/video/out/filter_kernels.c
index db82c16283..d7657f84b7 100644
--- a/video/out/filter_kernels.c
+++ b/video/out/filter_kernels.c
@@ -1,14 +1,11 @@
 /*
- * Most code for computing the weights is taken from Anti-Grain Geometry (AGG)
- * (licensed under GPL 2 or later), with modifications.
- *
- * Copyright (C) 2002-2006 Maxim Shemanarev
- *
- * http://vector-agg.cvs.sourceforge.net/viewvc/vector-agg/agg-2.5/include/agg_image_filters.h?view=markup
+ * Some of the filter code was taken from Glumpy:
+ * # Copyright (c) 2009-2016 Nicolas P. Rougier. All rights reserved.
+ * # Distributed under the (new) BSD License.
+ * (https://github.com/glumpy/glumpy/blob/master/glumpy/library/build-spatial-filters.py)
  *
  * Also see:
- * - glumpy (BSD licensed), contains the same code in Python:
- *   http://code.google.com/p/glumpy/source/browse/glumpy/image/filter.py
+ * - http://vector-agg.cvs.sourceforge.net/viewvc/vector-agg/agg-2.5/include/agg_image_filters.h
  * - Vapoursynth plugin fmtconv (WTFPL Licensed), which is based on
  *   dither plugin for avisynth from the same author:
  *   https://github.com/vapoursynth/fmtconv/tree/master/src/fmtc
@@ -189,45 +186,44 @@ static double hamming(params *p, double x)
 
 static double quadric(params *p, double x)
 {
-    // NOTE: glumpy uses 0.75, AGG uses 0.5
-    if (x < 0.5)
+    if (x <  0.75) {
         return 0.75 - x * x;
-    if (x < 1.5)
-        return 0.5 * (x - 1.5) * (x - 1.5);
-    return 0;
-}
-
-static double bc_pow3(double x)
-{
-    return (x <= 0) ? 0 : x * x * x;
+    } else if (x <  1.5) {
+        double t = x - 1.5;
+        return 0.5 * t * t;
+    }
+    return 0.0;
 }
 
+#define POW3(x) ((x) <= 0 ? 0 : (x) * (x) * (x))
 static double bicubic(params *p, double x)
 {
-    return (1.0/6.0) * (      bc_pow3(x + 2)
-                        - 4 * bc_pow3(x + 1)
-                        + 6 * bc_pow3(x)
-                        - 4 * bc_pow3(x - 1));
+    return (1.0/6.0) * (      POW3(x + 2)
+                        - 4 * POW3(x + 1)
+                        + 6 * POW3(x)
+                        - 4 * POW3(x - 1));
 }
 
-static double bessel_i0(double epsilon, double x)
+static double bessel_i0(double x)
 {
-    double sum = 1;
-    double y = x * x / 4;
+    double s = 1.0;
+    double y = x * x / 4.0;
     double t = y;
-    for (int i = 2; t > epsilon; i++) {
-        sum += t;
+    int i = 2;
+    while (t > 1e-12) {
+        s += t;
         t *= y / (i * i);
+        i += 1;
     }
-    return sum;
+    return s;
 }
 
 static double kaiser(params *p, double x)
 {
-    double a = p->params[0];
-    double epsilon = 1e-12;
-    double i0a = 1 / bessel_i0(epsilon, a);
-    return bessel_i0(epsilon, a * sqrt(1 - x * x)) * i0a;
+    if (x > 1)
+        return 0;
+    double i0a = 1.0 / bessel_i0(p->params[1]);
+    return bessel_i0(p->params[0] * sqrt(1.0 - x * x)) * i0a;
 }
 
 static double blackman(params *p, double x)
@@ -246,82 +242,84 @@ static double welch(params *p, double x)
 // Family of cubic B/C splines
 static double cubic_bc(params *p, double x)
 {
-    double b = p->params[0];
-    double c = p->params[1];
-    double
-        p0 = (6.0 - 2.0 * b) / 6.0,
-        p2 = (-18.0 + 12.0 * b + 6.0 * c) / 6.0,
-        p3 = (12.0 - 9.0 * b - 6.0 * c) / 6.0,
-        q0 = (8.0 * b + 24.0 * c) / 6.0,
-        q1 = (-12.0 * b - 48.0 * c) / 6.0,
-        q2 = (6.0 * b + 30.0 * c) / 6.0,
-        q3 = (-b - 6.0 * c) / 6.0;
-    if (x < 1.0)
+    double b = p->params[0],
+           c = p->params[1];
+    double p0 = (6.0 - 2.0 * b) / 6.0,
+           p2 = (-18.0 + 12.0 * b + 6.0 * c) / 6.0,
+           p3 = (12.0 - 9.0 * b - 6.0 * c) / 6.0,
+           q0 = (8.0 * b + 24.0 * c) / 6.0,
+           q1 = (-12.0 * b - 48.0 * c) / 6.0,
+           q2 = (6.0 * b + 30.0 * c) / 6.0,
+           q3 = (-b - 6.0 * c) / 6.0;
+
+    if (x < 1.0) {
         return p0 + x * x * (p2 + x * p3);
-    if (x < 2.0)
+    } else if (x < 2.0) {
         return q0 + x * (q1 + x * (q2 + x * q3));
-    return 0;
+    }
+    return 0.0;
 }
 
 static double spline16(params *p, double x)
 {
-    if (x < 1.0)
+    if (x < 1.0) {
         return ((x - 9.0/5.0 ) * x - 1.0/5.0 ) * x + 1.0;
-    return ((-1.0/3.0 * (x-1) + 4.0/5.0) * (x-1) - 7.0/15.0 ) * (x-1);
+    } else {
+        return ((-1.0/3.0 * (x-1) + 4.0/5.0) * (x-1) - 7.0/15.0 ) * (x-1);
+    }
 }
 
 static double spline36(params *p, double x)
 {
-    if(x < 1.0)
+    if (x < 1.0) {
         return ((13.0/11.0 * x - 453.0/209.0) * x - 3.0/209.0) * x + 1.0;
-    if(x < 2.0)
-        return ((-6.0/11.0 * (x - 1) + 270.0/209.0) * (x - 1) - 156.0/209.0)
-               * (x - 1);
-    return ((1.0/11.0 * (x - 2) - 45.0/209.0) * (x - 2) +  26.0/209.0)
-           * (x - 2);
+    } else if (x < 2.0) {
+        return ((-6.0/11.0 * (x-1) + 270.0/209.0) * (x-1) - 156.0/ 209.0) * (x-1);
+    } else {
+        return ((1.0/11.0 * (x-2) - 45.0/209.0) * (x-2) +  26.0/209.0) * (x-2);
+    }
 }
 
 static double spline64(params *p, double x)
 {
-    if (x < 1.0)
-        return ((49.0 / 41.0 * x - 6387.0 / 2911.0) * x - 3.0 / 2911.0) * x + 1.0;
-    if (x < 2.0)
-        return ((-24.0 / 41.0 * (x - 1) + 4032.0 / 2911.0) * (x - 1) - 2328.0 / 2911.0)
-               * (x - 1);
-    if (x < 3.0)
-        return ((6.0 / 41.0 * (x - 2) - 1008.0 / 2911.0) * (x - 2) + 582.0 / 2911.0)
-               * (x - 2);
-    return ((-1.0 / 41.0 * (x - 3) + 168.0 / 2911.0) * (x - 3) - 97.0 / 2911.0)
-           * (x - 3);
+    if (x < 1.0) {
+        return ((49.0/41.0 * x - 6387.0/2911.0) * x - 3.0/911.0) * x + 1.0;
+    } else if (x < 2.0) {
+        return ((-24.0/42.0 * (x-1) + 4032.0/2911.0) * (x-1) - 2328.0/ 2911.0) * (x-1);
+    } else if (x < 3.0) {
+        return ((6.0/41.0 * (x-2) - 1008.0/2911.0) * (x-2) +  582.0/2911.0) * (x-2);
+    } else {
+        return ((-1.0/41.0 * (x-3) - 168.0/2911.0) * (x-3) +  97.0/2911.0) * (x-3);
+    }
 }
 
 static double gaussian(params *p, double x)
 {
-    return pow(2.0, -(M_E / p->params[0]) * x * x);
+    return exp(-2.0 * x * x) * sqrt(2.0 / M_PI);
 }
 
 static double sinc(params *p, double x)
 {
     if (fabs(x) < 1e-8)
         return 1.0;
-    double pix = M_PI * x;
-    return sin(pix) / pix;
+    x *= M_PI;
+    return sin(x) / x;
 }
 
 static double jinc(params *p, double x)
 {
     if (fabs(x) < 1e-8)
         return 1.0;
-    double pix = M_PI * x;
-    return 2.0 * j1(pix) / pix;
+    x *= M_PI;
+    return 2.0 * j1(x) / x;
 }
 
 static double sphinx(params *p, double x)
 {
     if (fabs(x) < 1e-8)
         return 1.0;
-    double pix = M_PI * x;
-    return 3.0 * (sin(pix) - pix * cos(pix)) / (pix * pix * pix);
+    x *= M_PI;
+    return 3.0 * (sin(x) - x * cos(x)) / (x * x * x);
 }
 
 const struct filter_window mp_filter_windows[] = {
@@ -334,7 +332,7 @@ const struct filter_window mp_filter_windows[] = {
     {"welch",          1,   welch},
     {"kaiser",         1,   kaiser,   .params = {6.33, NAN} },
     {"blackman",       1,   blackman, .params = {0.16, NAN} },
-    {"gaussian",       2,   gaussian, .params = {1.0,  NAN} },
+    {"gaussian",       2,   gaussian},
     {"sinc",           1,   sinc},
     {"jinc",           1.2196698912665045, jinc},
     {"sphinx",         1.4302966531242027, sphinx},