libvo/filter_kernels.c

   1 /*
   2  * This file is part of mplayer2.
   3  *
   4  * Most code for computing the weights is taken from Anti-Grain Geometry (AGG)
   5  * (licensed under GPL 2 or later), with modifications.
   6  * Copyright (C) 2002-2006 Maxim Shemanarev
   7  * http://vector-agg.cvs.sourceforge.net/viewvc/vector-agg/agg-2.5/include/agg_image_filters.h?view=markup
   8  *
   9  * Also see glumpy (BSD licensed), contains the same code in Python:
  10  * http://code.google.com/p/glumpy/source/browse/glumpy/image/filter.py
  11  *
  12  * Also see: Paul Heckbert's "zoom"
  13  *
  14  * Also see XBMC: ConvolutionKernels.cpp etc.
  15  *
  16  * mplayer2 is free software; you can redistribute it and/or modify
  17  * it under the terms of the GNU General Public License as published by
  18  * the Free Software Foundation; either version 2 of the License, or
  19  * (at your option) any later version.
  20  *
  21  * mplayer2 is distributed in the hope that it will be useful,
  22  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  23  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  24  * GNU General Public License for more details.
  25  *
  26  * You should have received a copy of the GNU General Public License along
  27  * with mplayer2; if not, write to the Free Software Foundation, Inc.,
  28  * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
  29  */
  30
  31 #include <stddef.h>
  32 #include <string.h>
  33 #include <math.h>
  34 #include <assert.h>
  35
  36 #include "filter_kernels.h"
  37
  38 // NOTE: all filters are separable, symmetric, and are intended for use with
  39 //       a lookup table/texture.
  40
  41 const struct filter_kernel *mp_find_filter_kernel(const char *name)
  42 {
  43     for (const struct filter_kernel *k = mp_filter_kernels; k->name; k++) {
  44         if (strcmp(k->name, name) == 0)
  45             return k;
  46     }
  47     return NULL;
  48 }
  49
  50 // sizes = sorted list of available filter sizes, terminated with size 0
  51 // inv_scale = source_size / dest_size
  52 bool mp_init_filter(struct filter_kernel *filter, const int *sizes,
  53                     double inv_scale)
  54 {
  55     // only downscaling requires widening the filter
  56     filter->inv_scale = inv_scale >= 1.0 ? inv_scale : 1.0;
  57     double support = filter->radius * filter->inv_scale;
  58     int size = ceil(2.0 * support);
  59     // round up to smallest available size that's still large enough
  60     if (size < sizes[0])
  61         size = sizes[0];
  62     const int *cursize = sizes;
  63     while (size > *cursize && *cursize)
  64         cursize++;
  65     if (*cursize) {
  66         filter->size = *cursize;
  67         return true;
  68     } else {
  69         // The filter doesn't fit - instead of failing completely, use the
  70         // largest filter available. This is incorrect, but better than refusing
  71         // to do anything.
  72         filter->size = cursize[-1];
  73         filter->inv_scale = filter->size / 2.0 / filter->radius;
  74         return false;
  75     }
  76 }
  77
  78 // Calculate the 1D filtering kernel for N sample points.
  79 // N = number of samples, which is filter->size
  80 // The weights will be stored in out_w[0] to out_w[N - 1]
  81 // f = x0 - abs(x0), subpixel position in the range [0,1) or [0,1].
  82 void mp_compute_weights(struct filter_kernel *filter, double f, float *out_w)
  83 {
  84     assert(filter->size > 0);
  85     double sum = 0;
  86     for (int n = 0; n < filter->size; n++) {
  87         double x = f - (n - filter->size / 2 + 1);
  88         double w = filter->weight(filter, fabs(x) / filter->inv_scale);
  89         out_w[n] = w;
  90         sum += w;
  91     }
  92     //normalize
  93     for (int n = 0; n < filter->size; n++)
  94         out_w[n] /= sum;
  95 }
  96
  97 // Fill the given array with weights for the range [0.0, 1.0]. The array is
  98 // interpreted as rectangular array of count * filter->size items.
  99 void mp_compute_lut(struct filter_kernel *filter, int count, float *out_array)
 100 {
 101     for (int n = 0; n < count; n++) {
 102         mp_compute_weights(filter, n / (double)(count - 1),
 103                            out_array + filter->size * n);
 104     }
 105 }
 106
 107 typedef struct filter_kernel kernel;
 108
 109 static double bilinear(kernel *k, double x)
 110 {
 111     return 1.0 - x;
 112 }
 113
 114 static double hanning(kernel *k, double x)
 115 {
 116     return 0.5 + 0.5 * cos(M_PI * x);
 117 }
 118
 119 static double hamming(kernel *k, double x)
 120 {
 121     return 0.54 + 0.46 * cos(M_PI * x);
 122 }
 123
 124 static double hermite(kernel *k, double x)
 125 {
 126     return (2.0 * x - 3.0) * x * x + 1.0;
 127 }
 128
 129 static double quadric(kernel *k, double x)
 130 {
 131     // NOTE: glumpy uses 0.75, AGG uses 0.5
 132     if (x < 0.5)
 133         return 0.75 - x * x;
 134     if (x < 1.5)
 135         return 0.5 * (x - 1.5) * (x - 1.5);
 136     return 0;
 137 }
 138
 139 static double bc_pow3(double x)
 140 {
 141     return (x <= 0) ? 0 : x * x * x;
 142 }
 143
 144 static double bicubic(kernel *k, double x)
 145 {
 146     return (1.0/6.0) * (      bc_pow3(x + 2)
 147                         - 4 * bc_pow3(x + 1)
 148                         + 6 * bc_pow3(x)
 149                         - 4 * bc_pow3(x - 1));
 150 }
 151
 152 static double bessel_i0(double epsilon, double x)
 153 {
 154     double sum = 1;
 155     double y = x * x / 4;
 156     double t = y;
 157     for (int i = 2; t > epsilon; i++) {
 158         sum += t;
 159         t *= y / (i * i);
 160     }
 161     return sum;
 162 }
 163
 164 static double kaiser(kernel *k, double x)
 165 {
 166     double a = k->params[0];
 167     double b = k->params[1];
 168     double epsilon = 1e-12;
 169     double i0a = 1 / bessel_i0(epsilon, b);
 170     return bessel_i0(epsilon, a * sqrt(1 - x * x)) * i0a;
 171 }
 172
 173 static double catmull_rom(kernel *k, double x)
 174 {
 175     if (x < 1.0)
 176         return 0.5 * (2.0 + x * x * (-5.0 + x * 3.0));
 177     if (x < 2.0)
 178         return 0.5 * (4.0 + x * (-8.0 + x * (5.0 - x)));
 179     return 0;
 180 }
 181
 182 // Mitchell-Netravali
 183 static double mitchell(kernel *k, double x)
 184 {
 185     double b = k->params[0];
 186     double c = k->params[1];
 187     double
 188         p0 = (6.0 - 2.0 * b) / 6.0,
 189         p2 = (-18.0 + 12.0 * b + 6.0 * c) / 6.0,
 190         p3 = (12.0 - 9.0 * b - 6.0 * c) / 6.0,
 191         q0 = (8.0 * b + 24.0 * c) / 6.0,
 192         q1 = (-12.0 * b - 48.0 * c) / 6.0,
 193         q2 = (6.0 * b + 30.0 * c) / 6.0,
 194         q3 = (-b - 6.0 * c) / 6.0;
 195     if (x < 1.0)
 196         return p0 + x * x * (p2 + x * p3);
 197     if (x < 2.0)
 198         return q0 + x * (q1 + x * (q2 + x * q3));
 199     return 0;
 200 }
 201
 202 static double spline16(kernel *k, double x)
 203 {
 204     if (x < 1.0)
 205         return ((x - 9.0/5.0 ) * x - 1.0/5.0 ) * x + 1.0;
 206     return ((-1.0/3.0 * (x-1) + 4.0/5.0) * (x-1) - 7.0/15.0 ) * (x-1);
 207 }
 208
 209 static double spline36(kernel *k, double x)
 210 {
 211     if(x < 1.0)
 212         return ((13.0/11.0 * x - 453.0/209.0) * x - 3.0/209.0) * x + 1.0;
 213     if(x < 2.0)
 214         return ((-6.0/11.0 * (x - 1) + 270.0/209.0) * (x - 1) - 156.0/209.0)
 215                * (x - 1);
 216     return ((1.0/11.0 * (x - 2) - 45.0/209.0) * (x - 2) +  26.0/209.0)
 217            * (x - 2);
 218 }
 219
 220 static double gaussian(kernel *k, double x)
 221 {
 222     return exp(-2.0 * x * x) * sqrt(2.0 / M_PI);
 223 }
 224
 225 static double sinc(kernel *k, double x)
 226 {
 227     if (x == 0.0)
 228         return 1.0;
 229     double pix = M_PI * x;
 230     return sin(pix) / pix;
 231 }
 232
 233 static double lanczos(kernel *k, double x)
 234 {
 235     double radius = k->size / 2;
 236     if (x < -radius || x > radius)
 237         return 0;
 238     if (x == 0)
 239         return 1;
 240     double pix = M_PI * x;
 241     return radius * sin(pix) * sin(pix / radius) / (pix * pix);
 242 }
 243
 244 static double blackman(kernel *k, double x)
 245 {
 246     double radius = k->size / 2;
 247     if (x == 0.0)
 248         return 1.0;
 249     if (x > radius)
 250         return 0.0;
 251     x *= M_PI;
 252     double xr = x / radius;
 253     return (sin(x) / x) * (0.42 + 0.5 * cos(xr) + 0.08 * cos(2 * xr));
 254 }
 255
 256 const struct filter_kernel mp_filter_kernels[] = {
 257     {"bilinear_slow",  1,   bilinear},
 258     {"hanning",        1,   hanning},
 259     {"hamming",        1,   hamming},
 260     {"hermite",        1,   hermite},
 261     {"quadric",        1.5, quadric},
 262     {"bicubic",        2,   bicubic},
 263     {"kaiser",         1,   kaiser, .params = {6.33, 6.33} },
 264     {"catmull_rom",    2,   catmull_rom},
 265     {"mitchell",       2,   mitchell, .params = {1.0/3.0, 1.0/3.0} },
 266     {"spline16",       2,   spline16},
 267     {"spline36",       3,   spline36},
 268     {"gaussian",       2,   gaussian},
 269     {"sinc2",          2,   sinc},
 270     {"sinc3",          3,   sinc},
 271     {"sinc4",          4,   sinc},
 272     {"lanczos2",       2,   lanczos},
 273     {"lanczos3",       3,   lanczos},
 274     {"lanczos4",       4,   lanczos},
 275     {"blackman2",      2,   blackman},
 276     {"blackman3",      3,   blackman},
 277     {"blackman4",      4,   blackman},
 278     {0}
 279 };