1 /*
2 * Some of the filter code was taken from Glumpy:
3 * # Copyright (c) 2009-2016 Nicolas P. Rougier. All rights reserved.
4 * # Distributed under the (new) BSD License.
5 * (https://github.com/glumpy/glumpy/blob/master/glumpy/library/build-spatial-filters.py)
6 *
7 * Also see:
8 * - http://vector-agg.cvs.sourceforge.net/viewvc/vector-agg/agg-2.5/include/agg_image_filters.h
9 * - Vapoursynth plugin fmtconv (WTFPL Licensed), which is based on
10 * dither plugin for avisynth from the same author:
11 * https://github.com/vapoursynth/fmtconv/tree/master/src/fmtc
12 * - Paul Heckbert's "zoom"
13 * - XBMC: ConvolutionKernels.cpp etc.
14 *
15 * This file is part of mpv.
16 *
17 * This file can be distributed under the 3-clause license ("New BSD License").
18 *
19 * You can alternatively redistribute the non-Glumpy parts of this file and/or
20 * modify it under the terms of the GNU Lesser General Public
21 * License as published by the Free Software Foundation; either
22 * version 2.1 of the License, or (at your option) any later version.
23 */
24
25 #include <stddef.h>
26 #include <string.h>
27 #include <math.h>
28 #include <assert.h>
29
30 #include "filter_kernels.h"
31 #include "common/common.h"
32
33 // NOTE: all filters are designed for discrete convolution
34
mp_find_filter_window(const char * name)35 const struct filter_window *mp_find_filter_window(const char *name)
36 {
37 if (!name)
38 return NULL;
39 for (const struct filter_window *w = mp_filter_windows; w->name; w++) {
40 if (strcmp(w->name, name) == 0)
41 return w;
42 }
43 return NULL;
44 }
45
mp_find_filter_kernel(const char * name)46 const struct filter_kernel *mp_find_filter_kernel(const char *name)
47 {
48 if (!name)
49 return NULL;
50 for (const struct filter_kernel *k = mp_filter_kernels; k->f.name; k++) {
51 if (strcmp(k->f.name, name) == 0)
52 return k;
53 }
54 return NULL;
55 }
56
57 // sizes = sorted list of available filter sizes, terminated with size 0
58 // inv_scale = source_size / dest_size
mp_init_filter(struct filter_kernel * filter,const int * sizes,double inv_scale)59 bool mp_init_filter(struct filter_kernel *filter, const int *sizes,
60 double inv_scale)
61 {
62 assert(filter->f.radius > 0);
63 // Only downscaling requires widening the filter
64 filter->filter_scale = MPMAX(1.0, inv_scale);
65 double src_radius = filter->f.radius * filter->filter_scale;
66 // Polar filters are dependent solely on the radius
67 if (filter->polar) {
68 filter->size = 1; // Not meaningful for EWA/polar scalers.
69 // Safety precaution to avoid generating a gigantic shader
70 if (src_radius > 16.0) {
71 src_radius = 16.0;
72 filter->filter_scale = src_radius / filter->f.radius;
73 return false;
74 }
75 return true;
76 }
77 int size = ceil(2.0 * src_radius);
78 // round up to smallest available size that's still large enough
79 if (size < sizes[0])
80 size = sizes[0];
81 const int *cursize = sizes;
82 while (size > *cursize && *cursize)
83 cursize++;
84 if (*cursize) {
85 filter->size = *cursize;
86 return true;
87 } else {
88 // The filter doesn't fit - instead of failing completely, use the
89 // largest filter available. This is incorrect, but better than refusing
90 // to do anything.
91 filter->size = cursize[-1];
92 filter->filter_scale = (filter->size/2.0) / filter->f.radius;
93 return false;
94 }
95 }
96
97 // Sample from a blurred and tapered window
sample_window(struct filter_window * kernel,double x)98 static double sample_window(struct filter_window *kernel, double x)
99 {
100 if (!kernel->weight)
101 return 1.0;
102
103 // All windows are symmetric, this makes life easier
104 x = fabs(x);
105
106 // Stretch and taper the window size as needed
107 x = kernel->blur > 0.0 ? x / kernel->blur : x;
108 x = x <= kernel->taper ? 0.0 : (x - kernel->taper) / (1 - kernel->taper);
109
110 if (x < kernel->radius)
111 return kernel->weight(kernel, x);
112 return 0.0;
113 }
114
115 // Evaluate a filter's kernel and window at a given absolute position
sample_filter(struct filter_kernel * filter,double x)116 static double sample_filter(struct filter_kernel *filter, double x)
117 {
118 // The window is always stretched to the entire kernel
119 double w = sample_window(&filter->w, x / filter->f.radius * filter->w.radius);
120 double k = w * sample_window(&filter->f, x);
121 return k < 0 ? (1 - filter->clamp) * k : k;
122 }
123
124 // Calculate the 1D filtering kernel for N sample points.
125 // N = number of samples, which is filter->size
126 // The weights will be stored in out_w[0] to out_w[N - 1]
127 // f = x0 - abs(x0), subpixel position in the range [0,1) or [0,1].
mp_compute_weights(struct filter_kernel * filter,double f,float * out_w)128 static void mp_compute_weights(struct filter_kernel *filter, double f,
129 float *out_w)
130 {
131 assert(filter->size > 0);
132 double sum = 0;
133 for (int n = 0; n < filter->size; n++) {
134 double x = f - (n - filter->size / 2 + 1);
135 double w = sample_filter(filter, x / filter->filter_scale);
136 out_w[n] = w;
137 sum += w;
138 }
139 // Normalize to preserve energy
140 for (int n = 0; n < filter->size; n++)
141 out_w[n] /= sum;
142 }
143
144 // Fill the given array with weights for the range [0.0, 1.0]. The array is
145 // interpreted as rectangular array of count * filter->size items, with a
146 // stride of `stride` floats in between each array element. (For polar filters,
147 // the `count` indicates the row size and filter->size/stride are ignored)
148 //
149 // There will be slight sampling error if these weights are used in a OpenGL
150 // texture as LUT directly. The sampling point of a texel is located at its
151 // center, so out_array[0] will end up at 0.5 / count instead of 0.0.
152 // Correct lookup requires a linear coordinate mapping from [0.0, 1.0] to
153 // [0.5 / count, 1.0 - 0.5 / count].
mp_compute_lut(struct filter_kernel * filter,int count,int stride,float * out_array)154 void mp_compute_lut(struct filter_kernel *filter, int count, int stride,
155 float *out_array)
156 {
157 if (filter->polar) {
158 filter->radius_cutoff = 0.0;
159 // Compute a 1D array indexed by radius
160 for (int x = 0; x < count; x++) {
161 double r = x * filter->f.radius / (count - 1);
162 out_array[x] = sample_filter(filter, r);
163
164 if (fabs(out_array[x]) > filter->value_cutoff)
165 filter->radius_cutoff = r;
166 }
167 } else {
168 // Compute a 2D array indexed by subpixel position
169 for (int n = 0; n < count; n++) {
170 mp_compute_weights(filter, n / (double)(count - 1),
171 out_array + stride * n);
172 }
173 }
174 }
175
176 typedef struct filter_window params;
177
box(params * p,double x)178 static double box(params *p, double x)
179 {
180 // This is mathematically 1.0 everywhere, the clipping is done implicitly
181 // based on the radius.
182 return 1.0;
183 }
184
triangle(params * p,double x)185 static double triangle(params *p, double x)
186 {
187 return fmax(0.0, 1.0 - fabs(x / p->radius));
188 }
189
hanning(params * p,double x)190 static double hanning(params *p, double x)
191 {
192 return 0.5 + 0.5 * cos(M_PI * x);
193 }
194
hamming(params * p,double x)195 static double hamming(params *p, double x)
196 {
197 return 0.54 + 0.46 * cos(M_PI * x);
198 }
199
quadric(params * p,double x)200 static double quadric(params *p, double x)
201 {
202 if (x < 0.5) {
203 return 0.75 - x * x;
204 } else if (x < 1.5) {
205 double t = x - 1.5;
206 return 0.5 * t * t;
207 }
208 return 0.0;
209 }
210
211 #define POW3(x) ((x) <= 0 ? 0 : (x) * (x) * (x))
bicubic(params * p,double x)212 static double bicubic(params *p, double x)
213 {
214 return (1.0/6.0) * ( POW3(x + 2)
215 - 4 * POW3(x + 1)
216 + 6 * POW3(x)
217 - 4 * POW3(x - 1));
218 }
219
bessel_i0(double x)220 static double bessel_i0(double x)
221 {
222 double s = 1.0;
223 double y = x * x / 4.0;
224 double t = y;
225 int i = 2;
226 while (t > 1e-12) {
227 s += t;
228 t *= y / (i * i);
229 i += 1;
230 }
231 return s;
232 }
233
kaiser(params * p,double x)234 static double kaiser(params *p, double x)
235 {
236 if (x > 1)
237 return 0;
238 double i0a = 1.0 / bessel_i0(p->params[1]);
239 return bessel_i0(p->params[0] * sqrt(1.0 - x * x)) * i0a;
240 }
241
blackman(params * p,double x)242 static double blackman(params *p, double x)
243 {
244 double a = p->params[0];
245 double a0 = (1-a)/2.0, a1 = 1/2.0, a2 = a/2.0;
246 double pix = M_PI * x;
247 return a0 + a1*cos(pix) + a2*cos(2 * pix);
248 }
249
welch(params * p,double x)250 static double welch(params *p, double x)
251 {
252 return 1.0 - x*x;
253 }
254
255 // Family of cubic B/C splines
cubic_bc(params * p,double x)256 static double cubic_bc(params *p, double x)
257 {
258 double b = p->params[0],
259 c = p->params[1];
260 double p0 = (6.0 - 2.0 * b) / 6.0,
261 p2 = (-18.0 + 12.0 * b + 6.0 * c) / 6.0,
262 p3 = (12.0 - 9.0 * b - 6.0 * c) / 6.0,
263 q0 = (8.0 * b + 24.0 * c) / 6.0,
264 q1 = (-12.0 * b - 48.0 * c) / 6.0,
265 q2 = (6.0 * b + 30.0 * c) / 6.0,
266 q3 = (-b - 6.0 * c) / 6.0;
267
268 if (x < 1.0) {
269 return p0 + x * x * (p2 + x * p3);
270 } else if (x < 2.0) {
271 return q0 + x * (q1 + x * (q2 + x * q3));
272 }
273 return 0.0;
274 }
275
spline16(params * p,double x)276 static double spline16(params *p, double x)
277 {
278 if (x < 1.0) {
279 return ((x - 9.0/5.0 ) * x - 1.0/5.0 ) * x + 1.0;
280 } else {
281 return ((-1.0/3.0 * (x-1) + 4.0/5.0) * (x-1) - 7.0/15.0 ) * (x-1);
282 }
283 }
284
spline36(params * p,double x)285 static double spline36(params *p, double x)
286 {
287 if (x < 1.0) {
288 return ((13.0/11.0 * x - 453.0/209.0) * x - 3.0/209.0) * x + 1.0;
289 } else if (x < 2.0) {
290 return ((-6.0/11.0 * (x-1) + 270.0/209.0) * (x-1) - 156.0/ 209.0) * (x-1);
291 } else {
292 return ((1.0/11.0 * (x-2) - 45.0/209.0) * (x-2) + 26.0/209.0) * (x-2);
293 }
294 }
295
spline64(params * p,double x)296 static double spline64(params *p, double x)
297 {
298 if (x < 1.0) {
299 return ((49.0/41.0 * x - 6387.0/2911.0) * x - 3.0/2911.0) * x + 1.0;
300 } else if (x < 2.0) {
301 return ((-24.0/41.0 * (x-1) + 4032.0/2911.0) * (x-1) - 2328.0/2911.0) * (x-1);
302 } else if (x < 3.0) {
303 return ((6.0/41.0 * (x-2) - 1008.0/2911.0) * (x-2) + 582.0/2911.0) * (x-2);
304 } else {
305 return ((-1.0/41.0 * (x-3) + 168.0/2911.0) * (x-3) - 97.0/2911.0) * (x-3);
306 }
307 }
308
gaussian(params * p,double x)309 static double gaussian(params *p, double x)
310 {
311 return exp(-2.0 * x * x / p->params[0]);
312 }
313
sinc(params * p,double x)314 static double sinc(params *p, double x)
315 {
316 if (fabs(x) < 1e-8)
317 return 1.0;
318 x *= M_PI;
319 return sin(x) / x;
320 }
321
jinc(params * p,double x)322 static double jinc(params *p, double x)
323 {
324 if (fabs(x) < 1e-8)
325 return 1.0;
326 x *= M_PI;
327 return 2.0 * j1(x) / x;
328 }
329
sphinx(params * p,double x)330 static double sphinx(params *p, double x)
331 {
332 if (fabs(x) < 1e-8)
333 return 1.0;
334 x *= M_PI;
335 return 3.0 * (sin(x) - x * cos(x)) / (x * x * x);
336 }
337
338 const struct filter_window mp_filter_windows[] = {
339 {"box", 1, box},
340 {"triangle", 1, triangle},
341 {"bartlett", 1, triangle},
342 {"hanning", 1, hanning},
343 {"tukey", 1, hanning, .taper = 0.5},
344 {"hamming", 1, hamming},
345 {"quadric", 1.5, quadric},
346 {"welch", 1, welch},
347 {"kaiser", 1, kaiser, .params = {6.33, NAN} },
348 {"blackman", 1, blackman, .params = {0.16, NAN} },
349 {"gaussian", 2, gaussian, .params = {1.00, NAN} },
350 {"sinc", 1, sinc},
351 {"jinc", 1.2196698912665045, jinc},
352 {"sphinx", 1.4302966531242027, sphinx},
353 {0}
354 };
355
356 const struct filter_kernel mp_filter_kernels[] = {
357 // Spline filters
358 {{"spline16", 2, spline16}},
359 {{"spline36", 3, spline36}},
360 {{"spline64", 4, spline64}},
361 // Sinc filters
362 {{"sinc", 2, sinc, .resizable = true}},
363 {{"lanczos", 3, sinc, .resizable = true}, .window = "sinc"},
364 {{"ginseng", 3, sinc, .resizable = true}, .window = "jinc"},
365 // Jinc filters
366 {{"jinc", 3, jinc, .resizable = true}, .polar = true},
367 {{"ewa_lanczos", 3, jinc, .resizable = true}, .polar = true, .window = "jinc"},
368 {{"ewa_hanning", 3, jinc, .resizable = true}, .polar = true, .window = "hanning" },
369 {{"ewa_ginseng", 3, jinc, .resizable = true}, .polar = true, .window = "sinc"},
370 // Radius is based on the true jinc radius, slightly sharpened as per
371 // calculations by Nicolas Robidoux. Source: Imagemagick's magick/resize.c
372 {{"ewa_lanczossharp", 3.2383154841662362, jinc, .blur = 0.9812505644269356,
373 .resizable = true}, .polar = true, .window = "jinc"},
374 // Similar to the above, but softened instead. This one makes hash patterns
375 // disappear completely. Blur determined by trial and error.
376 {{"ewa_lanczossoft", 3.2383154841662362, jinc, .blur = 1.015,
377 .resizable = true}, .polar = true, .window = "jinc"},
378 // Very soft (blurred) hanning-windowed jinc; removes almost all aliasing.
379 // Blur paramater picked to match orthogonal and diagonal contributions
380 {{"haasnsoft", 3.2383154841662362, jinc, .blur = 1.11, .resizable = true},
381 .polar = true, .window = "hanning"},
382 // Cubic filters
383 {{"bicubic", 2, bicubic}},
384 {{"bcspline", 2, cubic_bc, .params = {0.5, 0.5} }},
385 {{"catmull_rom", 2, cubic_bc, .params = {0.0, 0.5} }},
386 {{"mitchell", 2, cubic_bc, .params = {1.0/3.0, 1.0/3.0} }},
387 {{"robidoux", 2, cubic_bc, .params = {12 / (19 + 9 * M_SQRT2),
388 113 / (58 + 216 * M_SQRT2)} }},
389 {{"robidouxsharp", 2, cubic_bc, .params = {6 / (13 + 7 * M_SQRT2),
390 7 / (2 + 12 * M_SQRT2)} }},
391 {{"ewa_robidoux", 2, cubic_bc, .params = {12 / (19 + 9 * M_SQRT2),
392 113 / (58 + 216 * M_SQRT2)}},
393 .polar = true},
394 {{"ewa_robidouxsharp", 2,cubic_bc, .params = {6 / (13 + 7 * M_SQRT2),
395 7 / (2 + 12 * M_SQRT2)}},
396 .polar = true},
397 // Miscellaneous filters
398 {{"box", 1, box, .resizable = true}},
399 {{"nearest", 0.5, box}},
400 {{"triangle", 1, triangle, .resizable = true}},
401 {{"gaussian", 2, gaussian, .params = {1.0, NAN}, .resizable = true}},
402 {{0}}
403 };
404