Something went wrong. Try again.
🎥 Command line media player
Something went wrong. Try again.
123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404/* * 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: * - 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 * - Paul Heckbert's "zoom" * - XBMC: ConvolutionKernels.cpp etc. * * This file is part of mpv. * * This file can be distributed under the 3-clause license ("New BSD License"). * * You can alternatively redistribute the non-Glumpy parts of this file and/or * modify it under the terms of the GNU Lesser General Public * License as published by the Free Software Foundation; either * version 2.1 of the License, or (at your option) any later version. */
#include <stddef.h>#include <string.h>#include <math.h>#include <assert.h>
#include "filter_kernels.h"#include "common/common.h"
// NOTE: all filters are designed for discrete convolution
const struct filter_window *mp_find_filter_window(const char *name){ if (!name) return NULL; for (const struct filter_window *w = mp_filter_windows; w->name; w++) { if (strcmp(w->name, name) == 0) return w; } return NULL;}
const struct filter_kernel *mp_find_filter_kernel(const char *name){ if (!name) return NULL; for (const struct filter_kernel *k = mp_filter_kernels; k->f.name; k++) { if (strcmp(k->f.name, name) == 0) return k; } return NULL;}
// sizes = sorted list of available filter sizes, terminated with size 0// inv_scale = source_size / dest_sizebool mp_init_filter(struct filter_kernel *filter, const int *sizes, double inv_scale){ assert(filter->f.radius > 0); // Only downscaling requires widening the filter filter->filter_scale = MPMAX(1.0, inv_scale); double src_radius = filter->f.radius * filter->filter_scale; // Polar filters are dependent solely on the radius if (filter->polar) { filter->size = 1; // Not meaningful for EWA/polar scalers. // Safety precaution to avoid generating a gigantic shader if (src_radius > 16.0) { src_radius = 16.0; filter->filter_scale = src_radius / filter->f.radius; return false; } return true; } int size = ceil(2.0 * src_radius); // round up to smallest available size that's still large enough if (size < sizes[0]) size = sizes[0]; const int *cursize = sizes; while (size > *cursize && *cursize) cursize++; if (*cursize) { filter->size = *cursize; return true; } else { // The filter doesn't fit - instead of failing completely, use the // largest filter available. This is incorrect, but better than refusing // to do anything. filter->size = cursize[-1]; filter->filter_scale = (filter->size/2.0) / filter->f.radius; return false; }}
// Sample from a blurred and tapered windowstatic double sample_window(struct filter_window *kernel, double x){ if (!kernel->weight) return 1.0;
// All windows are symmetric, this makes life easier x = fabs(x);
// Stretch and taper the window size as needed x = kernel->blur > 0.0 ? x / kernel->blur : x; x = x <= kernel->taper ? 0.0 : (x - kernel->taper) / (1 - kernel->taper);
if (x < kernel->radius) return kernel->weight(kernel, x); return 0.0;}
// Evaluate a filter's kernel and window at a given absolute positionstatic double sample_filter(struct filter_kernel *filter, double x){ // The window is always stretched to the entire kernel double w = sample_window(&filter->w, x / filter->f.radius * filter->w.radius); double k = w * sample_window(&filter->f, x); return k < 0 ? (1 - filter->clamp) * k : k;}
// Calculate the 1D filtering kernel for N sample points.// N = number of samples, which is filter->size// The weights will be stored in out_w[0] to out_w[N - 1]// f = x0 - abs(x0), subpixel position in the range [0,1) or [0,1].static void mp_compute_weights(struct filter_kernel *filter, double f, float *out_w){ assert(filter->size > 0); double sum = 0; for (int n = 0; n < filter->size; n++) { double x = f - (n - filter->size / 2 + 1); double w = sample_filter(filter, x / filter->filter_scale); out_w[n] = w; sum += w; } // Normalize to preserve energy for (int n = 0; n < filter->size; n++) out_w[n] /= sum;}
// Fill the given array with weights for the range [0.0, 1.0]. The array is// interpreted as rectangular array of count * filter->size items, with a// stride of `stride` floats in between each array element. (For polar filters,// the `count` indicates the row size and filter->size/stride are ignored)//// There will be slight sampling error if these weights are used in a OpenGL// texture as LUT directly. The sampling point of a texel is located at its// center, so out_array[0] will end up at 0.5 / count instead of 0.0.// Correct lookup requires a linear coordinate mapping from [0.0, 1.0] to// [0.5 / count, 1.0 - 0.5 / count].void mp_compute_lut(struct filter_kernel *filter, int count, int stride, float *out_array){ if (filter->polar) { filter->radius_cutoff = 0.0; // Compute a 1D array indexed by radius for (int x = 0; x < count; x++) { double r = x * filter->f.radius / (count - 1); out_array[x] = sample_filter(filter, r);
if (fabs(out_array[x]) > filter->value_cutoff) filter->radius_cutoff = r; } } else { // Compute a 2D array indexed by subpixel position for (int n = 0; n < count; n++) { mp_compute_weights(filter, n / (double)(count - 1), out_array + stride * n); } }}
typedef struct filter_window params;
static double box(params *p, double x){ // This is mathematically 1.0 everywhere, the clipping is done implicitly // based on the radius. return 1.0;}
static double triangle(params *p, double x){ return fmax(0.0, 1.0 - fabs(x / p->radius));}
static double hanning(params *p, double x){ return 0.5 + 0.5 * cos(M_PI * x);}
static double hamming(params *p, double x){ return 0.54 + 0.46 * cos(M_PI * x);}
static double quadric(params *p, double x){ if (x < 0.75) { return 0.75 - 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) * ( POW3(x + 2) - 4 * POW3(x + 1) + 6 * POW3(x) - 4 * POW3(x - 1));}
static double bessel_i0(double x){ double s = 1.0; double y = x * x / 4.0; double t = y; int i = 2; while (t > 1e-12) { s += t; t *= y / (i * i); i += 1; } return s;}
static double kaiser(params *p, double x){ 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){ double a = p->params[0]; double a0 = (1-a)/2.0, a1 = 1/2.0, a2 = a/2.0; double pix = M_PI * x; return a0 + a1*cos(pix) + a2*cos(2 * pix);}
static double welch(params *p, double x){ return 1.0 - x*x;}
// Family of cubic B/C splinesstatic double cubic_bc(params *p, double x){ 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); } else if (x < 2.0) { return q0 + x * (q1 + x * (q2 + x * q3)); } return 0.0;}
static double spline16(params *p, double x){ if (x < 1.0) { return ((x - 9.0/5.0 ) * x - 1.0/5.0 ) * x + 1.0; } 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) { return ((13.0/11.0 * x - 453.0/209.0) * x - 3.0/209.0) * x + 1.0; } 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; } else if (x < 2.0) { return ((-24.0/41.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 exp(-2.0 * x * x / p->params[0]);}
static double sinc(params *p, double x){ if (fabs(x) < 1e-8) return 1.0; x *= M_PI; return sin(x) / x;}
static double jinc(params *p, double x){ if (fabs(x) < 1e-8) return 1.0; x *= M_PI; return 2.0 * j1(x) / x;}
static double sphinx(params *p, double x){ if (fabs(x) < 1e-8) return 1.0; x *= M_PI; return 3.0 * (sin(x) - x * cos(x)) / (x * x * x);}
const struct filter_window mp_filter_windows[] = { {"box", 1, box}, {"triangle", 1, triangle}, {"bartlett", 1, triangle}, {"hanning", 1, hanning}, {"tukey", 1, hanning, .taper = 0.5}, {"hamming", 1, hamming}, {"quadric", 1.5, quadric}, {"welch", 1, welch}, {"kaiser", 1, kaiser, .params = {6.33, NAN} }, {"blackman", 1, blackman, .params = {0.16, NAN} }, {"gaussian", 2, gaussian, .params = {1.00, NAN} }, {"sinc", 1, sinc}, {"jinc", 1.2196698912665045, jinc}, {"sphinx", 1.4302966531242027, sphinx}, {0}};
const struct filter_kernel mp_filter_kernels[] = { // Spline filters {{"spline16", 2, spline16}}, {{"spline36", 3, spline36}}, {{"spline64", 4, spline64}}, // Sinc filters {{"sinc", 2, sinc, .resizable = true}}, {{"lanczos", 3, sinc, .resizable = true}, .window = "sinc"}, {{"ginseng", 3, sinc, .resizable = true}, .window = "jinc"}, // Jinc filters {{"jinc", 3, jinc, .resizable = true}, .polar = true}, {{"ewa_lanczos", 3, jinc, .resizable = true}, .polar = true, .window = "jinc"}, {{"ewa_hanning", 3, jinc, .resizable = true}, .polar = true, .window = "hanning" }, {{"ewa_ginseng", 3, jinc, .resizable = true}, .polar = true, .window = "sinc"}, // Radius is based on the true jinc radius, slightly sharpened as per // calculations by Nicolas Robidoux. Source: Imagemagick's magick/resize.c {{"ewa_lanczossharp", 3.2383154841662362, jinc, .blur = 0.9812505644269356, .resizable = true}, .polar = true, .window = "jinc"}, // Similar to the above, but softened instead. This one makes hash patterns // disappear completely. Blur determined by trial and error. {{"ewa_lanczossoft", 3.2383154841662362, jinc, .blur = 1.015, .resizable = true}, .polar = true, .window = "jinc"}, // Very soft (blurred) hanning-windowed jinc; removes almost all aliasing. // Blur paramater picked to match orthogonal and diagonal contributions {{"haasnsoft", 3.2383154841662362, jinc, .blur = 1.11, .resizable = true}, .polar = true, .window = "hanning"}, // Cubic filters {{"bicubic", 2, bicubic}}, {{"bcspline", 2, cubic_bc, .params = {0.5, 0.5} }}, {{"catmull_rom", 2, cubic_bc, .params = {0.0, 0.5} }}, {{"mitchell", 2, cubic_bc, .params = {1.0/3.0, 1.0/3.0} }}, {{"robidoux", 2, cubic_bc, .params = {12 / (19 + 9 * M_SQRT2), 113 / (58 + 216 * M_SQRT2)} }}, {{"robidouxsharp", 2, cubic_bc, .params = {6 / (13 + 7 * M_SQRT2), 7 / (2 + 12 * M_SQRT2)} }}, {{"ewa_robidoux", 2, cubic_bc, .params = {12 / (19 + 9 * M_SQRT2), 113 / (58 + 216 * M_SQRT2)}}, .polar = true}, {{"ewa_robidouxsharp", 2,cubic_bc, .params = {6 / (13 + 7 * M_SQRT2), 7 / (2 + 12 * M_SQRT2)}}, .polar = true}, // Miscellaneous filters {{"box", 1, box, .resizable = true}}, {{"nearest", 0.5, box}}, {{"triangle", 1, triangle, .resizable = true}}, {{"gaussian", 2, gaussian, .params = {1.0, NAN}, .resizable = true}}, {{0}}};