vo_opengl: add ewa_lanczos upscaler (aka jinc)

This is the polar (elliptic weighted average) version of lanczos.
This introduces a general new form of polar filters.

1 files changed, 71 insertions, 0 deletions

diff --git a/video/out/filter_kernels.c b/video/out/filter_kernels.c
index 26f62ff567..4e5ca471bb 100644
--- a/video/out/filter_kernels.c
+++ b/video/out/filter_kernels.c
@@ -58,6 +58,11 @@ bool mp_init_filter(struct filter_kernel *filter, const int *sizes,
 {
     if (filter->radius < 0)
         filter->radius = 2.0;
+    // polar filters can be of any radius, and nothing special is needed
+    if (filter->polar) {
+        filter->size = filter->radius;
+        return true;
+    }
     // only downscaling requires widening the filter
     filter->inv_scale = inv_scale >= 1.0 ? inv_scale : 1.0;
     double support = filter->radius * filter->inv_scale;
@@ -111,6 +116,18 @@ void mp_compute_lut(struct filter_kernel *filter, int count, float *out_array)
     }
 }
 
+// Fill the given array with weights for the range [0, R], where R is the
+// radius of hte filter. The array is interpreted as a one-dimensional array
+// of count items.
+void mp_compute_lut_polar(struct filter_kernel *filter, int count, float *out_array)
+{
+    assert(filter->radius > 0);
+    for (int x = 0; x < count; x++) {
+        double r = x * filter->radius / (count - 1);
+        out_array[x] = r <= filter->radius ? filter->weight(filter, r) : 0;
+    }
+}
+
 typedef struct filter_kernel kernel;
 
 static double nearest(kernel *k, double x)
@@ -261,6 +278,14 @@ static double sinc(kernel *k, double x)
     return sin(pix) / pix;
 }
 
+static double jinc(kernel *k, double x)
+{
+    if (x == 0.0)
+        return 1.0;
+    double pix = M_PI * x;
+    return 2.0 * j1(pix) / pix;
+}
+
 static double lanczos(kernel *k, double x)
 {
     double radius = k->size / 2;
@@ -272,6 +297,48 @@ static double lanczos(kernel *k, double x)
     return radius * sin(pix) * sin(pix / radius) / (pix * pix);
 }
 
+static double ewa_lanczos(kernel *k, double x)
+{
+    double radius = k->radius;
+    assert(radius >= 1.0);
+
+    // This is already three orders of magnitude slower than anything you could
+    // possibly hope to play back in realtime and results in tons of ringing
+    // artifacts, so I doubt anybody will complain.
+    if (radius > 16)
+        radius = 16;
+
+    if (fabs(x) < 1e-8)
+        return 1.0;
+    if (fabs(x) >= radius)
+        return 0.0;
+
+    // Precomputed zeros of the jinc() function, needed to adjust the
+    // window size. Computing this at runtime is nontrivial.
+    // Copied from: https://github.com/AviSynth/jinc-resize/blob/master/JincResize/JincFilter.cpp#L171
+    static double jinc_zeros[16] = {
+        1.2196698912665045,
+        2.2331305943815286,
+        3.2383154841662362,
+        4.2410628637960699,
+        5.2427643768701817,
+        6.2439216898644877,
+        7.2447598687199570,
+        8.2453949139520427,
+        9.2458926849494673,
+        10.246293348754916,
+        11.246622794877883,
+        12.246898461138105,
+        13.247132522181061,
+        14.247333735806849,
+        15.247508563037300,
+        16.247661874700962
+    };
+
+    double window = jinc_zeros[0] / jinc_zeros[(int)radius - 1];
+    return jinc(k, x) * jinc(k, x*window);
+}
+
 static double blackman(kernel *k, double x)
 {
     double radius = k->size / 2;
@@ -303,6 +370,10 @@ const struct filter_kernel mp_filter_kernels[] = {
     {"sinc3",          3,   sinc},
     {"sinc4",          4,   sinc},
     {"sinc",           -1,  sinc},
+    {"ewa_lanczos2",   2,   ewa_lanczos, .polar = true},
+    {"ewa_lanczos3",   3,   ewa_lanczos, .polar = true},
+    {"ewa_lanczos4",   4,   ewa_lanczos, .polar = true},
+    {"ewa_lanczos",    -1,  ewa_lanczos, .polar = true},
     {"lanczos2",       2,   lanczos},
     {"lanczos3",       3,   lanczos},
     {"lanczos4",       4,   lanczos},