FFmpeg coverage

Directory: ../../../ffmpeg/
File: src/libswscale/filters.c
Date: 2026-05-04 20:54:34
Exec Total Coverage
Lines: 120 214 56.1%
Functions: 8 24 33.3%
Branches: 55 94 58.5%
Line Branch Exec Source
1 /*
2 * Copyright (C) 2026 Niklas Haas
3 *
4 * This file is part of FFmpeg.
5 *
6 * FFmpeg is free software; you can redistribute it and/or
7 * modify it under the terms of the GNU Lesser General Public
8 * License as published by the Free Software Foundation; either
9 * version 2.1 of the License, or (at your option) any later version.
10 *
11 * FFmpeg is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 * Lesser General Public License for more details.
15 *
16 * You should have received a copy of the GNU Lesser General Public
17 * License along with FFmpeg; if not, write to the Free Software
18 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
19 */
20
21 #include <math.h>
22 #include <stdbool.h>
23
24 #include <libavutil/attributes.h>
25 #include <libavutil/avassert.h>
26 #include <libavutil/mem.h>
27
28 #include "filters.h"
29
30 #ifdef _WIN32
31 # define j1 _j1
32 #endif
33
34 /* Maximum (pre-stretching) radius (for tunable filters) */
35 #define RADIUS_MAX 10.0
36
37 /* Defined only on [0, radius]. */
38 typedef double (*SwsFilterKernel)(double x, const double *params);
39
40 typedef struct SwsFilterFunction {
41 char name[16];
42 double radius; /* negative means resizable */
43 SwsFilterKernel kernel;
44 SwsFilterKernel window; /* optional */
45 double params[SWS_NUM_SCALER_PARAMS]; /* default params */
46 } SwsFilterFunction;
47
48 static const SwsFilterFunction filter_functions[SWS_SCALE_NB];
49
50 984032 static double scaler_sample(const SwsFilterFunction *f, double x)
51 {
52 984032 x = fabs(x);
53
2/2
✓ Branch 0 taken 27328 times.
✓ Branch 1 taken 956704 times.
 984032 if (x > f->radius)
54 27328 return 0.0;
55
56 956704 double w = f->kernel(x, f->params);
57
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 956704 times.
 956704 if (f->window)
58 ✗ w *= f->window(x / f->radius, f->params);
59 956704 return w;
60 }
61
62 59136 static void compute_row(SwsFilterWeights *f, const SwsFilterFunction *fun,
63 double radius, double ratio_inv, double stretch_inv,
64 int dst_pos, double *tmp)
65 {
66 59136 int *out = &f->weights[dst_pos * f->filter_size];
67 59136 int *pos = &f->offsets[dst_pos];
68
69 /**
70 * Explanation of the 0.5 offsets: Normally, pixel samples are assumed
71 * to be representative of the center of their containing area; e.g. for
72 * a 2x2 image, the samples are located at {0.5, 1.5}^2. However, with
73 * integer indexing, we round sample positions down (0-based indexing).
74 * So the (0, 0) sample is actually located at (0.5, 0.5) and represents
75 * the entire square from (0,0) to (1,1). When normalizing between different
76 * image sizes, we therefore need to add/subtract off these 0.5 offsets.
77 */
78 59136 const double src_pos = (dst_pos + 0.5) * ratio_inv - 0.5;
79
2/2
✓ Branch 0 taken 34048 times.
✓ Branch 1 taken 25088 times.
 59136 if (f->filter_size == 1) {
80 34048 *pos = fmin(fmax(round(src_pos), 0.0), f->src_size - 1);
81 34048 *out = SWS_FILTER_SCALE;
82 34048 return;
83 }
84
85 /* First pixel that is actually within the filter envelope */
86 25088 const double start_pos = src_pos - radius;
87 25088 int64_t start_idx = ceil(start_pos);
88 25088 start_idx = FFMAX(start_idx, 0); /* edge clamping */
89 25088 start_idx = FFMIN(start_idx, f->src_size - f->filter_size);
90 25088 const double offset = start_idx - src_pos;
91 25088 *pos = start_idx;
92
93 /**
94 * Generate raw filter weights with maximum precision. Sum the positive
95 * and negative weights separately to avoid catastrophic cancellation. This
96 * summation order should already give the best precision because abs(w)
97 * is monotonically decreasing
98 */
99 25088 const double base = stretch_inv * offset;
100 25088 double wsum_pos = 0.0, wsum_neg = 0.0;
101
2/2
✓ Branch 0 taken 277760 times.
✓ Branch 1 taken 25088 times.
 302848 for (int i = 0; i < f->filter_size; i++) {
102 277760 tmp[i] = scaler_sample(fun, base + stretch_inv * i);
103
2/2
✓ Branch 0 taken 164304 times.
✓ Branch 1 taken 113456 times.
 277760 if (tmp[i] >= 0)
104 164304 wsum_pos += tmp[i];
105 else
106 113456 wsum_neg += tmp[i];
107 }
108
109 25088 const double wsum = wsum_pos + wsum_neg;
110
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 25088 times.
 25088 av_assert0(wsum > 0);
111
112 /* Generate correctly rounded filter weights with error diffusion */
113 25088 double error = 0.0;
114 25088 int sum_pos = 0, sum_neg = 0;
115
2/2
✓ Branch 0 taken 277760 times.
✓ Branch 1 taken 25088 times.
 302848 for (int i = 0; i < f->filter_size; i++) {
116
2/2
✓ Branch 0 taken 25088 times.
✓ Branch 1 taken 252672 times.
 277760 if (i == f->filter_size - 1) {
117 /* Ensure weights sum to exactly SWS_FILTER_SCALE */
118 25088 out[i] = SWS_FILTER_SCALE - sum_pos - sum_neg;
119 } else {
120 252672 const double w = tmp[i] / wsum + error;
121 252672 out[i] = round(w * SWS_FILTER_SCALE);
122 252672 error = w - (double) out[i] / SWS_FILTER_SCALE;
123 }
124
2/2
✓ Branch 0 taken 199920 times.
✓ Branch 1 taken 77840 times.
 277760 if (out[i] >= 0)
125 199920 sum_pos += out[i];
126 else
127 77840 sum_neg += out[i];
128 }
129
130
2/2
✓ Branch 0 taken 1512 times.
✓ Branch 1 taken 23576 times.
 25088 if (sum_pos > f->sum_positive)
131 1512 f->sum_positive = sum_pos;
132
2/2
✓ Branch 0 taken 1400 times.
✓ Branch 1 taken 23688 times.
 25088 if (sum_neg < f->sum_negative)
133 1400 f->sum_negative = sum_neg;
134 }
135
136 1344 static void sws_filter_free(AVRefStructOpaque opaque, void *obj)
137 {
138 1344 SwsFilterWeights *filter = obj;
139 1344 av_refstruct_unref(&filter->weights);
140 1344 av_refstruct_unref(&filter->offsets);
141 1344 }
142
143 1344 static bool validate_params(const SwsFilterFunction *fun, SwsScaler scaler)
144 {
145
1/4
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 1344 times.
 1344 switch (scaler) {
146 ✗ case SWS_SCALE_GAUSSIAN:
147 ✗ return fun->params[0] >= 0.0; /* sigma */
148 ✗ case SWS_SCALE_LANCZOS:
149 ✗ return fun->params[0] >= 1.0 && fun->params[0] <= RADIUS_MAX; /* radius */
150 ✗ case SWS_SCALE_BICUBIC:
151 ✗ return fun->params[0] < 3.0; /* B param (division by zero) */
152 1344 default:
153 1344 return true;
154 }
155 }
156
157 1344 static double filter_radius(const SwsFilterFunction *fun)
158 {
159 1344 const double bound = fun->radius;
160 1344 const double step = 1e-2;
161
162 1344 double radius = bound;
163 1344 double prev = 0.0, fprev = 1.0; /* f(0) is always 1.0 */
164 1344 double integral = 0.0;
165
2/2
✓ Branch 0 taken 706272 times.
✓ Branch 1 taken 1344 times.
 707616 for (double x = step; x < bound + step; x += step) {
166 706272 const double fx = scaler_sample(fun, x);
167 706272 integral += (fprev + fx) * step; /* trapezoidal rule (mirrored) */
168 706272 double cutoff = SWS_MAX_REDUCE_CUTOFF * integral;
169
8/8
✓ Branch 0 taken 363552 times.
✓ Branch 1 taken 342720 times.
✓ Branch 2 taken 359520 times.
✓ Branch 3 taken 4032 times.
✓ Branch 4 taken 328608 times.
✓ Branch 5 taken 373632 times.
✓ Branch 6 taken 3360 times.
✓ Branch 7 taken 325248 times.
 706272 if ((fprev > cutoff && fx <= cutoff) || (fprev < -cutoff && fx >= -cutoff)) {
170 /* estimate crossing with secant method; note that we have to
171 * bias by the cutoff to find the actual cutoff radius */
172
2/2
✓ Branch 0 taken 3360 times.
✓ Branch 1 taken 4032 times.
 7392 double estimate = fx + (fx > fprev ? cutoff : -cutoff);
173 7392 double root = x - estimate * (x - prev) / (fx - fprev);
174 7392 radius = fmin(root, bound);
175 }
176 706272 prev = x;
177 706272 fprev = fx;
178 }
179
180 1344 return radius;
181 }
182
183 1344 int ff_sws_filter_generate(void *log, const SwsFilterParams *params,
184 SwsFilterWeights **out)
185 {
186 1344 SwsScaler scaler = params->scaler;
187
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1344 times.
 1344 if (scaler >= SWS_SCALE_NB)
188 ✗ return AVERROR(EINVAL);
189
190
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1344 times.
 1344 if (scaler == SWS_SCALE_AUTO)
191 ✗ scaler = SWS_SCALE_BICUBIC;
192
193 1344 const double ratio = (double) params->dst_size / params->src_size;
194 1344 double stretch = 1.0;
195
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 1344 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 1344 if (ratio < 1.0 && scaler != SWS_SCALE_POINT) {
196 /* Widen filter for downscaling (anti-aliasing) */
197 ✗ stretch = 1.0 / ratio;
198 }
199
200
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1344 times.
 1344 if (scaler == SWS_SCALE_AREA) {
201 /**
202 * SWS_SCALE_AREA is a pseudo-filter that is equivalent to bilinear
203 * filtering for upscaling (since bilinear just evenly mixes samples
204 * according to the relative distance), and equivalent to (anti-aliased)
205 * point sampling for downscaling.
206 */
207 ✗ scaler = ratio >= 1.0 ? SWS_SCALE_BILINEAR : SWS_SCALE_POINT;
208 }
209
210 1344 SwsFilterFunction fun = filter_functions[scaler];
211
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1344 times.
 1344 if (!fun.kernel)
212 ✗ return AVERROR(EINVAL);
213
214
2/2
✓ Branch 0 taken 2688 times.
✓ Branch 1 taken 1344 times.
 4032 for (int i = 0; i < SWS_NUM_SCALER_PARAMS; i++) {
215
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2688 times.
 2688 if (params->scaler_params[i] != SWS_PARAM_DEFAULT)
216 ✗ fun.params[i] = params->scaler_params[i];
217 }
218
219
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1344 times.
 1344 if (!validate_params(&fun, scaler)) {
220 ✗ av_log(log, AV_LOG_ERROR, "Invalid parameters for scaler %s: {%f, %f}\n",
221 fun.name, fun.params[0], fun.params[1]);
222 ✗ return AVERROR(EINVAL);
223 }
224
225
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1344 times.
 1344 if (fun.radius < 0.0) /* tunable width kernels like lanczos */
226 ✗ fun.radius = fun.params[0];
227
228 1344 const double radius = filter_radius(&fun) * stretch;
229 1344 int filter_size = ceil(radius * 2.0);
230 1344 filter_size = FFMIN(filter_size, params->src_size);
231
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1344 times.
 1344 av_assert0(filter_size >= 1);
232
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1344 times.
 1344 if (filter_size > SWS_FILTER_SIZE_MAX)
233 ✗ return AVERROR(ENOTSUP);
234
235 SwsFilterWeights *filter;
236 1344 filter = av_refstruct_alloc_ext(sizeof(*filter), 0, NULL, sws_filter_free);
237
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1344 times.
 1344 if (!filter)
238 ✗ return AVERROR(ENOMEM);
239 1344 memcpy(filter->name, fun.name, sizeof(filter->name));
240 1344 filter->src_size = params->src_size;
241 1344 filter->dst_size = params->dst_size;
242 1344 filter->filter_size = filter_size;
243
2/2
✓ Branch 0 taken 784 times.
✓ Branch 1 taken 560 times.
 1344 if (filter->filter_size == 1)
244 784 filter->sum_positive = SWS_FILTER_SCALE;
245
246 1344 av_log(log, AV_LOG_DEBUG, "Generating %s filter with %d taps (radius = %f)\n",
247 1344 filter->name, filter->filter_size, radius);
248
249 1344 filter->num_weights = (size_t) params->dst_size * filter->filter_size;
250 1344 filter->weights = av_refstruct_allocz(filter->num_weights * sizeof(*filter->weights));
251
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1344 times.
 1344 if (!filter->weights) {
252 ✗ av_refstruct_unref(&filter);
253 ✗ return AVERROR(ENOMEM);
254 }
255
256 1344 filter->offsets = av_refstruct_allocz(params->dst_size * sizeof(*filter->offsets));
257
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1344 times.
 1344 if (!filter->offsets) {
258 ✗ av_refstruct_unref(&filter);
259 ✗ return AVERROR(ENOMEM);
260 }
261
262 1344 double *tmp = av_malloc(filter->filter_size * sizeof(*tmp));
263
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1344 times.
 1344 if (!tmp) {
264 ✗ av_refstruct_unref(&filter);
265 ✗ return AVERROR(ENOMEM);
266 }
267
268 1344 const double ratio_inv = 1.0 / ratio, stretch_inv = 1.0 / stretch;
269
2/2
✓ Branch 0 taken 59136 times.
✓ Branch 1 taken 1344 times.
 60480 for (int i = 0; i < params->dst_size; i++)
270 59136 compute_row(filter, &fun, radius, ratio_inv, stretch_inv, i, tmp);
271 1344 av_free(tmp);
272
273 1344 *out = filter;
274 1344 return 0;
275 }
276
277 /*
278 * Some of the filter code originally derives (via libplacebo/mpv) from Glumpy:
279 * # Copyright (c) 2009-2016 Nicolas P. Rougier. All rights reserved.
280 * # Distributed under the (new) BSD License.
281 * (https://github.com/glumpy/glumpy/blob/master/glumpy/library/build-spatial-filters.py)
282 *
283 * The math underlying each filter function was written from scratch, with
284 * some algorithms coming from a number of different sources, including:
285 * - https://en.wikipedia.org/wiki/Window_function
286 * - https://en.wikipedia.org/wiki/Jinc
287 * - http://vector-agg.cvs.sourceforge.net/viewvc/vector-agg/agg-2.5/include/agg_image_filters.h
288 * - Vapoursynth plugin fmtconv (WTFPL Licensed), which is based on
289 * dither plugin for avisynth from the same author:
290 * https://github.com/vapoursynth/fmtconv/tree/master/src/fmtc
291 * - Paul Heckbert's "zoom"
292 * - XBMC: ConvolutionKernels.cpp etc.
293 * - https://github.com/AviSynth/jinc-resize (only used to verify the math)
294 */
295
296 32928 av_unused static double box(double x, const double *params)
297 {
298 32928 return 1.0;
299 }
300
301 ✗ av_unused static double triangle(double x, const double *params)
302 {
303 ✗ return 1.0 - x;
304 }
305
306 ✗ av_unused static double cosine(double x, const double *params)
307 {
308 ✗ return cos(x);
309 }
310
311 ✗ av_unused static double hann(double x, const double *params)
312 {
313 ✗ return 0.5 + 0.5 * cos(M_PI * x);
314 }
315
316 ✗ av_unused static double hamming(double x, const double *params)
317 {
318 ✗ return 0.54 + 0.46 * cos(M_PI * x);
319 }
320
321 ✗ av_unused static double welch(double x, const double *params)
322 {
323 ✗ return 1.0 - x * x;
324 }
325
326 ✗ av_unused static double bessel_i0(double x)
327 {
328 ✗ double s = 1.0;
329 ✗ double y = x * x / 4.0;
330 ✗ double t = y;
331 ✗ int i = 2;
332 ✗ while (t > 1e-12) {
333 ✗ s += t;
334 ✗ t *= y / (i * i);
335 ✗ i += 1;
336 }
337 ✗ return s;
338 }
339
340 ✗ av_unused static double kaiser(double x, const double *params)
341 {
342 ✗ double alpha = fmax(params[0], 0.0);
343 ✗ double scale = bessel_i0(alpha);
344 ✗ return bessel_i0(alpha * sqrt(1.0 - x * x)) / scale;
345 }
346
347 ✗ av_unused static double blackman(double x, const double *params)
348 {
349 ✗ double a = params[0];
350 ✗ double a0 = (1 - a) / 2.0, a1 = 1 / 2.0, a2 = a / 2.0;
351 ✗ x *= M_PI;
352 ✗ return a0 + a1 * cos(x) + a2 * cos(2 * x);
353 }
354
355 ✗ av_unused static double bohman(double x, const double *params)
356 {
357 ✗ double pix = M_PI * x;
358 ✗ return (1.0 - x) * cos(pix) + sin(pix) / M_PI;
359 }
360
361 ✗ av_unused static double gaussian(double x, const double *params)
362 {
363 ✗ return exp(-params[0] * x * x);
364 }
365
366 ✗ av_unused static double quadratic(double x, const double *params)
367 {
368 ✗ if (x < 0.5) {
369 ✗ return 1.0 - 4.0/3.0 * (x * x);
370 } else {
371 ✗ return 2.0 / 3.0 * (x - 1.5) * (x - 1.5);
372 }
373 }
374
375 923776 av_unused static double sinc(double x, const double *params)
376 {
377
2/2
✓ Branch 0 taken 4480 times.
✓ Branch 1 taken 919296 times.
 923776 if (x < 1e-8)
378 4480 return 1.0;
379 919296 x *= M_PI;
380 919296 return sin(x) / x;
381 }
382
383 ✗ av_unused static double jinc(double x, const double *params)
384 {
385 ✗ if (x < 1e-8)
386 ✗ return 1.0;
387 ✗ x *= M_PI;
388 ✗ return 2.0 * j1(x) / x;
389 }
390
391 ✗ av_unused static double sphinx(double x, const double *params)
392 {
393 ✗ if (x < 1e-8)
394 ✗ return 1.0;
395 ✗ x *= M_PI;
396 ✗ return 3.0 * (sin(x) - x * cos(x)) / (x * x * x);
397 }
398
399 ✗ av_unused static double cubic(double x, const double *params)
400 {
401 ✗ const double b = params[0], c = params[1];
402 ✗ double p0 = 6.0 - 2.0 * b,
403 ✗ p2 = -18.0 + 12.0 * b + 6.0 * c,
404 ✗ p3 = 12.0 - 9.0 * b - 6.0 * c,
405 ✗ q0 = 8.0 * b + 24.0 * c,
406 ✗ q1 = -12.0 * b - 48.0 * c,
407 ✗ q2 = 6.0 * b + 30.0 * c,
408 ✗ q3 = -b - 6.0 * c;
409
410 ✗ if (x < 1.0) {
411 ✗ return (p0 + x * x * (p2 + x * p3)) / p0;
412 } else {
413 ✗ return (q0 + x * (q1 + x * (q2 + x * q3))) / p0;
414 }
415 }
416
417 ✗ static double spline_coeff(double a, double b, double c, double d, double x)
418 {
419 ✗ if (x <= 1.0) {
420 ✗ return ((d * x + c) * x + b) * x + a;
421 } else {
422 ✗ return spline_coeff(0.0,
423 ✗ b + 2.0 * c + 3.0 * d,
424 ✗ c + 3.0 * d,
425 ✗ -b - 3.0 * c - 6.0 * d,
426 x - 1.0);
427 }
428 }
429
430 ✗ av_unused static double spline(double x, const double *params)
431 {
432 ✗ const double p = -2.196152422706632;
433 ✗ return spline_coeff(1.0, 0.0, p, -p - 1.0, x);
434 }
435
436 static const SwsFilterFunction filter_functions[SWS_SCALE_NB] = {
437 [SWS_SCALE_BILINEAR] = { "bilinear", 1.0, triangle },
438 [SWS_SCALE_BICUBIC] = { "bicubic", 2.0, cubic, .params = { 0.0, 0.6 } },
439 [SWS_SCALE_POINT] = { "point", 0.5, box },
440 [SWS_SCALE_GAUSSIAN] = { "gaussian", 4.0, gaussian, .params = { 3.0 } },
441 [SWS_SCALE_SINC] = { "sinc", RADIUS_MAX, sinc },
442 [SWS_SCALE_LANCZOS] = { "lanczos", -1.0, sinc, sinc, .params = { 3.0 } },
443 [SWS_SCALE_SPLINE] = { "spline", RADIUS_MAX, spline },
444 /* SWS_SCALE_AREA is a pseudo-filter, see code above */
445 };
446