FFmpeg coverage

Directory: ../../../ffmpeg/
File: src/libswscale/filters.c
Date: 2026-10-29 23:19:09
Exec Total Coverage
Lines: 154 229 67.2%
Functions: 10 24 41.7%
Branches: 69 99 69.7%
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 179005272 static double scaler_sample(const SwsFilterFunction *f, double x)
51 {
52 179005272 x = fabs(x);
53
2/2
✓ Branch 0 taken 3731059 times.
✓ Branch 1 taken 175274213 times.
 179005272 if (x > f->radius)
54 3731059 return 0.0;
55
56 175274213 double w = f->kernel(x, f->params);
57
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 175274213 times.
 175274213 if (f->window)
58 ✗ w *= f->window(x / f->radius, f->params);
59 175274213 return w;
60 }
61
62 52298976 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 52298976 int *out = &f->weights[dst_pos * f->filter_size];
67 52298976 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 52298976 const double src_pos = (dst_pos + 0.5) * ratio_inv - 0.5 + f->offset;
79
2/2
✓ Branch 0 taken 2612160 times.
✓ Branch 1 taken 49686816 times.
 52298976 if (f->filter_size == 1) {
80 2612160 *pos = fmin(fmax(round(src_pos), 0.0), f->src_size - 1);
81 2612160 *out = SWS_FILTER_SCALE;
82 2612160 return;
83 }
84
85 /* First pixel that is actually within the filter envelope */
86 49686816 const double start_pos = src_pos - radius;
87 49686816 int64_t start_idx = ceil(start_pos);
88 49686816 start_idx = FFMAX(start_idx, 0); /* edge clamping */
89 49686816 start_idx = FFMIN(start_idx, f->src_size - f->filter_size);
90 49686816 const double offset = start_idx - src_pos;
91 49686816 *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 49686816 const double base = stretch_inv * offset;
100 49686816 double wsum_pos = 0.0, wsum_neg = 0.0;
101
2/2
✓ Branch 0 taken 137989056 times.
✓ Branch 1 taken 49686816 times.
 187675872 for (int i = 0; i < f->filter_size; i++) {
102 137989056 tmp[i] = scaler_sample(fun, base + stretch_inv * i);
103
2/2
✓ Branch 0 taken 106112652 times.
✓ Branch 1 taken 31876404 times.
 137989056 if (tmp[i] >= 0)
104 106112652 wsum_pos += tmp[i];
105 else
106 31876404 wsum_neg += tmp[i];
107 }
108
109 49686816 const double wsum = wsum_pos + wsum_neg;
110
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 49686816 times.
 49686816 av_assert0(wsum > 0);
111
112 /* Generate correctly rounded filter weights with error diffusion */
113 49686816 double error = 0.0;
114 49686816 int sum_pos = 0, sum_neg = 0;
115
2/2
✓ Branch 0 taken 137989056 times.
✓ Branch 1 taken 49686816 times.
 187675872 for (int i = 0; i < f->filter_size; i++) {
116
2/2
✓ Branch 0 taken 49686816 times.
✓ Branch 1 taken 88302240 times.
 137989056 if (i == f->filter_size - 1) {
117 /* Ensure weights sum to exactly SWS_FILTER_SCALE */
118 49686816 out[i] = SWS_FILTER_SCALE - sum_pos - sum_neg;
119 } else {
120 88302240 const double w = tmp[i] / wsum + error;
121 88302240 out[i] = round(w * SWS_FILTER_SCALE);
122 88302240 error = w - (double) out[i] / SWS_FILTER_SCALE;
123 }
124
2/2
✓ Branch 0 taken 106219500 times.
✓ Branch 1 taken 31769556 times.
 137989056 if (out[i] >= 0)
125 106219500 sum_pos += out[i];
126 else
127 31769556 sum_neg += out[i];
128 }
129
130
2/2
✓ Branch 0 taken 1496648 times.
✓ Branch 1 taken 48190168 times.
 49686816 if (sum_pos > f->sum_positive)
131 1496648 f->sum_positive = sum_pos;
132
2/2
✓ Branch 0 taken 518949 times.
✓ Branch 1 taken 49167867 times.
 49686816 if (sum_neg < f->sum_negative)
133 518949 f->sum_negative = sum_neg;
134 }
135
136 1202597 static void sws_filter_free(AVRefStructOpaque opaque, void *obj)
137 {
138 1202597 SwsFilterWeights *filter = obj;
139 1202597 av_refstruct_unref(&filter->weights);
140 1202597 av_refstruct_unref(&filter->offsets);
141 1202597 }
142
143 1202597 static bool validate_params(const SwsFilterFunction *fun, SwsScaler scaler)
144 {
145
2/4
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✓ Branch 2 taken 194991 times.
✓ Branch 3 taken 1007606 times.
 1202597 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 194991 case SWS_SCALE_BICUBIC:
151 194991 return fun->params[0] < 3.0; /* B param (division by zero) */
152 1007606 default:
153 1007606 return true;
154 }
155 }
156
157 /**
158 * Numerically estimate the last intersection between the function value
159 * and the cutoff domain [-SWS_MAX_REDUCE_CUTOFF, SWS_MAX_REDUCE_CUTOFF].
160 */
161 197007 static double est_filter_radius(const SwsFilterFunction *fun)
162 {
163 197007 const double bound = fun->radius;
164 197007 const double step = 1e-2;
165
166 197007 double radius = bound;
167 197007 double prev = 0.0, fprev = 1.0; /* f(0) is always 1.0 */
168 197007 double integral = 0.0;
169
2/2
✓ Branch 0 taken 41016216 times.
✓ Branch 1 taken 197007 times.
 41213223 for (double x = step; x < bound + step; x += step) {
170 41016216 const double fx = scaler_sample(fun, x);
171 41016216 integral += (fprev + fx) * step; /* trapezoidal rule (mirrored) */
172 41016216 double cutoff = SWS_MAX_REDUCE_CUTOFF * integral;
173
8/8
✓ Branch 0 taken 20488956 times.
✓ Branch 1 taken 20527260 times.
✓ Branch 2 taken 20283885 times.
✓ Branch 3 taken 205071 times.
✓ Branch 4 taken 19314978 times.
✓ Branch 5 taken 21496167 times.
✓ Branch 6 taken 205071 times.
✓ Branch 7 taken 19109907 times.
 41016216 if ((fprev > cutoff && fx <= cutoff) || (fprev < -cutoff && fx >= -cutoff)) {
174 /* estimate crossing with secant method; note that we have to
175 * bias by the cutoff to find the actual cutoff radius */
176
2/2
✓ Branch 0 taken 205071 times.
✓ Branch 1 taken 205071 times.
 410142 double estimate = fx + (fx > fprev ? cutoff : -cutoff);
177 410142 double root = x - estimate * (x - prev) / (fx - fprev);
178 410142 radius = fmin(root, bound);
179 }
180 41016216 prev = x;
181 41016216 fprev = fx;
182 }
183
184 197007 return radius;
185 }
186
187 1202597 int ff_sws_filter_generate(void *log, const SwsFilterParams *params,
188 SwsFilterWeights **out)
189 {
190 1202597 SwsScaler scaler = params->scaler;
191
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1202597 times.
 1202597 if (scaler >= SWS_SCALE_NB)
192 ✗ return AVERROR(EINVAL);
193
194
2/2
✓ Branch 0 taken 73992 times.
✓ Branch 1 taken 1128605 times.
 1202597 if (scaler == SWS_SCALE_AUTO)
195 73992 scaler = SWS_SCALE_BICUBIC;
196
197 1202597 double virtual_size = params->virtual_size;
198
1/2
✓ Branch 0 taken 1202597 times.
✗ Branch 1 not taken.
 1202597 if (!virtual_size)
199 1202597 virtual_size = params->dst_size;
200
201 1202597 const double ratio = virtual_size / params->src_size;
202 1202597 double stretch = 1.0;
203
4/4
✓ Branch 0 taken 107997 times.
✓ Branch 1 taken 1094600 times.
✓ Branch 2 taken 94794 times.
✓ Branch 3 taken 13203 times.
 1202597 if (ratio < 1.0 && scaler != SWS_SCALE_POINT) {
204 /* Widen filter for downscaling (anti-aliasing) */
205 94794 stretch = 1.0 / ratio;
206 }
207
208
2/2
✓ Branch 0 taken 25398 times.
✓ Branch 1 taken 1177199 times.
 1202597 if (scaler == SWS_SCALE_AREA) {
209 /**
210 * SWS_SCALE_AREA is a pseudo-filter that is equivalent to bilinear
211 * filtering for upscaling (since bilinear just evenly mixes samples
212 * according to the relative distance), and equivalent to (anti-aliased)
213 * point sampling for downscaling.
214 */
215
2/2
✓ Branch 0 taken 12792 times.
✓ Branch 1 taken 12606 times.
 25398 scaler = ratio >= 1.0 ? SWS_SCALE_BILINEAR : SWS_SCALE_POINT;
216 }
217
218 1202597 SwsFilterFunction fun = filter_functions[scaler];
219
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1202597 times.
 1202597 if (!fun.kernel)
220 ✗ return AVERROR(EINVAL);
221
222
2/2
✓ Branch 0 taken 2405194 times.
✓ Branch 1 taken 1202597 times.
 3607791 for (int i = 0; i < SWS_NUM_SCALER_PARAMS; i++) {
223
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2405194 times.
 2405194 if (params->scaler_params[i] != SWS_PARAM_DEFAULT)
224 ✗ fun.params[i] = params->scaler_params[i];
225 }
226
227
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1202597 times.
 1202597 if (!validate_params(&fun, scaler)) {
228 ✗ av_log(log, AV_LOG_ERROR, "Invalid parameters for scaler %s: {%f, %f}\n",
229 fun.name, fun.params[0], fun.params[1]);
230 ✗ return AVERROR(EINVAL);
231 }
232
233
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1202597 times.
 1202597 if (fun.radius < 0.0) /* tunable width kernels like lanczos */
234 ✗ fun.radius = fun.params[0];
235
236 double radius;
237
3/3
✓ Branch 0 taken 40833 times.
✓ Branch 1 taken 964757 times.
✓ Branch 2 taken 197007 times.
 1202597 switch (scaler) {
238 40833 case SWS_SCALE_POINT:
239 40833 radius = 0.5;
240 40833 break;
241 964757 case SWS_SCALE_BILINEAR:
242 964757 radius = 1.0 - SWS_MAX_REDUCE_CUTOFF;
243 964757 break;
244 197007 default:
245 /* Numerically estimate radius of nontrivial or parametric kernels */
246 197007 radius = est_filter_radius(&fun);
247 197007 break;
248 }
249 1202597 radius *= stretch;
250
251 1202597 int filter_size = ceil(radius * 2.0);
252 1202597 filter_size = FFMIN(filter_size, params->src_size);
253
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1202597 times.
 1202597 av_assert0(filter_size >= 1);
254
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1202597 times.
 1202597 if (filter_size > SWS_FILTER_SIZE_MAX)
255 ✗ return AVERROR(ENOTSUP);
256
257 SwsFilterWeights *filter;
258 1202597 filter = av_refstruct_alloc_ext(sizeof(*filter), 0, NULL, sws_filter_free);
259
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1202597 times.
 1202597 if (!filter)
260 ✗ return AVERROR(ENOMEM);
261 1202597 memcpy(filter->name, fun.name, sizeof(filter->name));
262 1202597 filter->src_size = params->src_size;
263 1202597 filter->dst_size = params->dst_size;
264 1202597 filter->virtual_size = virtual_size;
265 1202597 filter->offset = params->offset;
266 1202597 filter->filter_size = filter_size;
267
2/2
✓ Branch 0 taken 28563 times.
✓ Branch 1 taken 1174034 times.
 1202597 if (filter->filter_size == 1)
268 28563 filter->sum_positive = SWS_FILTER_SCALE;
269
270 1202597 av_log(log, AV_LOG_DEBUG, "Generating %s filter with %d taps (radius = %f)\n",
271 1202597 filter->name, filter->filter_size, radius);
272
273 1202597 filter->num_weights = (size_t) params->dst_size * filter->filter_size;
274 1202597 filter->weights = av_refstruct_allocz(filter->num_weights * sizeof(*filter->weights));
275
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1202597 times.
 1202597 if (!filter->weights) {
276 ✗ av_refstruct_unref(&filter);
277 ✗ return AVERROR(ENOMEM);
278 }
279
280 1202597 filter->offsets = av_refstruct_allocz(params->dst_size * sizeof(*filter->offsets));
281
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1202597 times.
 1202597 if (!filter->offsets) {
282 ✗ av_refstruct_unref(&filter);
283 ✗ return AVERROR(ENOMEM);
284 }
285
286 1202597 double *tmp = av_malloc(filter->filter_size * sizeof(*tmp));
287
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1202597 times.
 1202597 if (!tmp) {
288 ✗ av_refstruct_unref(&filter);
289 ✗ return AVERROR(ENOMEM);
290 }
291
292 1202597 const double ratio_inv = 1.0 / ratio, stretch_inv = 1.0 / stretch;
293
2/2
✓ Branch 0 taken 52298976 times.
✓ Branch 1 taken 1202597 times.
 53501573 for (int i = 0; i < params->dst_size; i++)
294 52298976 compute_row(filter, &fun, radius, ratio_inv, stretch_inv, i, tmp);
295 1202597 av_free(tmp);
296
297 1202597 *out = filter;
298 1202597 return 0;
299 }
300
301 /*
302 * Some of the filter code originally derives (via libplacebo/mpv) from Glumpy:
303 * # Copyright (c) 2009-2016 Nicolas P. Rougier. All rights reserved.
304 * # Distributed under the (new) BSD License.
305 * (https://github.com/glumpy/glumpy/blob/master/glumpy/library/build-spatial-filters.py)
306 *
307 * The math underlying each filter function was written from scratch, with
308 * some algorithms coming from a number of different sources, including:
309 * - https://en.wikipedia.org/wiki/Window_function
310 * - https://en.wikipedia.org/wiki/Jinc
311 * - http://vector-agg.cvs.sourceforge.net/viewvc/vector-agg/agg-2.5/include/agg_image_filters.h
312 * - Vapoursynth plugin fmtconv (WTFPL Licensed), which is based on
313 * dither plugin for avisynth from the same author:
314 * https://github.com/vapoursynth/fmtconv/tree/master/src/fmtc
315 * - Paul Heckbert's "zoom"
316 * - XBMC: ConvolutionKernels.cpp etc.
317 * - https://github.com/AviSynth/jinc-resize (only used to verify the math)
318 */
319
320 1613568 av_unused static double box(double x, const double *params)
321 {
322 1613568 return 1.0;
323 }
324
325 68483606 av_unused static double triangle(double x, const double *params)
326 {
327 68483606 return 1.0 - x;
328 }
329
330 ✗ av_unused static double cosine(double x, const double *params)
331 {
332 ✗ return cos(x);
333 }
334
335 ✗ av_unused static double hann(double x, const double *params)
336 {
337 ✗ return 0.5 + 0.5 * cos(M_PI * x);
338 }
339
340 ✗ av_unused static double hamming(double x, const double *params)
341 {
342 ✗ return 0.54 + 0.46 * cos(M_PI * x);
343 }
344
345 ✗ av_unused static double welch(double x, const double *params)
346 {
347 ✗ return 1.0 - x * x;
348 }
349
350 ✗ av_unused static double bessel_i0(double x)
351 {
352 ✗ double s = 1.0;
353 ✗ double y = x * x / 4.0;
354 ✗ double t = y;
355 ✗ int i = 2;
356 ✗ while (t > 1e-12) {
357 ✗ s += t;
358 ✗ t *= y / (i * i);
359 ✗ i += 1;
360 }
361 ✗ return s;
362 }
363
364 ✗ av_unused static double kaiser(double x, const double *params)
365 {
366 ✗ double alpha = fmax(params[0], 0.0);
367 ✗ double scale = bessel_i0(alpha);
368 ✗ return bessel_i0(alpha * sqrt(1.0 - x * x)) / scale;
369 }
370
371 ✗ av_unused static double blackman(double x, const double *params)
372 {
373 ✗ double a = params[0];
374 ✗ double a0 = (1 - a) / 2.0, a1 = 1 / 2.0, a2 = a / 2.0;
375 ✗ x *= M_PI;
376 ✗ return a0 + a1 * cos(x) + a2 * cos(2 * x);
377 }
378
379 ✗ av_unused static double bohman(double x, const double *params)
380 {
381 ✗ double pix = M_PI * x;
382 ✗ return (1.0 - x) * cos(pix) + sin(pix) / M_PI;
383 }
384
385 ✗ av_unused static double gaussian(double x, const double *params)
386 {
387 ✗ return exp(-params[0] * x * x);
388 }
389
390 ✗ av_unused static double quadratic(double x, const double *params)
391 {
392 ✗ if (x < 0.5) {
393 ✗ return 1.0 - 4.0/3.0 * (x * x);
394 } else {
395 ✗ return 2.0 / 3.0 * (x - 1.5) * (x - 1.5);
396 }
397 }
398
399 2771328 av_unused static double sinc(double x, const double *params)
400 {
401
2/2
✓ Branch 0 taken 13440 times.
✓ Branch 1 taken 2757888 times.
 2771328 if (x < 1e-8)
402 13440 return 1.0;
403 2757888 x *= M_PI;
404 2757888 return sin(x) / x;
405 }
406
407 ✗ av_unused static double jinc(double x, const double *params)
408 {
409 ✗ if (x < 1e-8)
410 ✗ return 1.0;
411 ✗ x *= M_PI;
412 ✗ return 2.0 * j1(x) / x;
413 }
414
415 ✗ av_unused static double sphinx(double x, const double *params)
416 {
417 ✗ if (x < 1e-8)
418 ✗ return 1.0;
419 ✗ x *= M_PI;
420 ✗ return 3.0 * (sin(x) - x * cos(x)) / (x * x * x);
421 }
422
423 102405711 av_unused static double cubic(double x, const double *params)
424 {
425 102405711 const double b = params[0], c = params[1];
426 102405711 double p0 = 6.0 - 2.0 * b,
427 102405711 p2 = -18.0 + 12.0 * b + 6.0 * c,
428 102405711 p3 = 12.0 - 9.0 * b - 6.0 * c,
429 102405711 q0 = 8.0 * b + 24.0 * c,
430 102405711 q1 = -12.0 * b - 48.0 * c,
431 102405711 q2 = 6.0 * b + 30.0 * c,
432 102405711 q3 = -b - 6.0 * c;
433
434
2/2
✓ Branch 0 taken 51370575 times.
✓ Branch 1 taken 51035136 times.
 102405711 if (x < 1.0) {
435 51370575 return (p0 + x * x * (p2 + x * p3)) / p0;
436 } else {
437 51035136 return (q0 + x * (q1 + x * (q2 + x * q3))) / p0;
438 }
439 }
440
441 ✗ static double spline_coeff(double a, double b, double c, double d, double x)
442 {
443 ✗ if (x <= 1.0) {
444 ✗ return ((d * x + c) * x + b) * x + a;
445 } else {
446 ✗ return spline_coeff(0.0,
447 ✗ b + 2.0 * c + 3.0 * d,
448 ✗ c + 3.0 * d,
449 ✗ -b - 3.0 * c - 6.0 * d,
450 x - 1.0);
451 }
452 }
453
454 ✗ av_unused static double spline(double x, const double *params)
455 {
456 ✗ const double p = -2.196152422706632;
457 ✗ return spline_coeff(1.0, 0.0, p, -p - 1.0, x);
458 }
459
460 static const SwsFilterFunction filter_functions[SWS_SCALE_NB] = {
461 [SWS_SCALE_BILINEAR] = { "bilinear", 1.0, triangle },
462 [SWS_SCALE_BICUBIC] = { "bicubic", 2.0, cubic, .params = { 0.0, 0.6 } },
463 [SWS_SCALE_POINT] = { "point", 0.5, box },
464 [SWS_SCALE_GAUSSIAN] = { "gaussian", 4.0, gaussian, .params = { 3.0 } },
465 [SWS_SCALE_SINC] = { "sinc", RADIUS_MAX, sinc },
466 [SWS_SCALE_LANCZOS] = { "lanczos", -1.0, sinc, sinc, .params = { 3.0 } },
467 [SWS_SCALE_SPLINE] = { "spline", RADIUS_MAX, spline },
468 /* SWS_SCALE_AREA is a pseudo-filter, see code above */
469 };
470