filters/resample_functions.cpp
| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | // Avisynth v2.5. Copyright 2002 Ben Rudiak-Gould et al. | ||
| 2 | // http://avisynth.nl | ||
| 3 | |||
| 4 | // This program is free software; you can redistribute it and/or modify | ||
| 5 | // it under the terms of the GNU General Public License as published by | ||
| 6 | // the Free Software Foundation; either version 2 of the License, or | ||
| 7 | // (at your option) any later version. | ||
| 8 | // | ||
| 9 | // This program is distributed in the hope that it will be useful, | ||
| 10 | // but WITHOUT ANY WARRANTY; without even the implied warranty of | ||
| 11 | // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | ||
| 12 | // GNU General Public License for more details. | ||
| 13 | // | ||
| 14 | // You should have received a copy of the GNU General Public License | ||
| 15 | // along with this program; if not, write to the Free Software | ||
| 16 | // Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA, or visit | ||
| 17 | // http://www.gnu.org/copyleft/gpl.html . | ||
| 18 | // | ||
| 19 | // Linking Avisynth statically or dynamically with other modules is making a | ||
| 20 | // combined work based on Avisynth. Thus, the terms and conditions of the GNU | ||
| 21 | // General Public License cover the whole combination. | ||
| 22 | // | ||
| 23 | // As a special exception, the copyright holders of Avisynth give you | ||
| 24 | // permission to link Avisynth with independent modules that communicate with | ||
| 25 | // Avisynth solely through the interfaces defined in avisynth.h, regardless of the license | ||
| 26 | // terms of these independent modules, and to copy and distribute the | ||
| 27 | // resulting combined work under terms of your choice, provided that | ||
| 28 | // every copy of the combined work is accompanied by a complete copy of | ||
| 29 | // the source code of Avisynth (the version of Avisynth used to produce the | ||
| 30 | // combined work), being distributed under the terms of the GNU General | ||
| 31 | // Public License plus this exception. An independent module is a module | ||
| 32 | // which is not derived from or based on Avisynth, such as 3rd-party filters, | ||
| 33 | // import and export plugins, or graphical user interfaces. | ||
| 34 | |||
| 35 | #include "resample_functions.h" | ||
| 36 | #include <cmath> | ||
| 37 | #include <vector> | ||
| 38 | #include <algorithm> | ||
| 39 | #include <avs/minmax.h> | ||
| 40 | #include <avs/alignment.h> | ||
| 41 | |||
| 42 | |||
| 43 | /******************************************* | ||
| 44 | *************************************** | ||
| 45 | ** Helper classes for resample.cpp ** | ||
| 46 | *************************************** | ||
| 47 | *******************************************/ | ||
| 48 | |||
| 49 | /*************************** | ||
| 50 | ***** Point filter ***** | ||
| 51 | **************************/ | ||
| 52 | |||
| 53 | 333 | double PointFilter::f(double x) { | |
| 54 | AVS_UNUSED(x); | ||
| 55 | 333 | return 1.0; | |
| 56 | } | ||
| 57 | |||
| 58 | |||
| 59 | /*************************** | ||
| 60 | ***** Triangle filter ***** | ||
| 61 | **************************/ | ||
| 62 | |||
| 63 | 8759 | double TriangleFilter::f(double x) { | |
| 64 | 8759 | x = fabs(x); | |
| 65 |
2/2✓ Branch 2 → 3 taken 7781 times.
✓ Branch 2 → 4 taken 978 times.
|
8759 | return (x<1.0) ? 1.0-x : 0.0; |
| 66 | } | ||
| 67 | |||
| 68 | |||
| 69 | |||
| 70 | |||
| 71 | |||
| 72 | /********************************* | ||
| 73 | *** Mitchell-Netravali filter *** | ||
| 74 | *********************************/ | ||
| 75 | |||
| 76 | 6 | MitchellNetravaliFilter::MitchellNetravaliFilter (double b, double c) { | |
| 77 | 6 | p0 = ( 6. - 2.*b ) / 6.; | |
| 78 | 6 | p2 = ( -18. + 12.*b + 6.*c ) / 6.; | |
| 79 | 6 | p3 = ( 12. - 9.*b - 6.*c ) / 6.; | |
| 80 | 6 | q0 = ( 8.*b + 24.*c ) / 6.; | |
| 81 | 6 | q1 = ( - 12.*b - 48.*c ) / 6.; | |
| 82 | 6 | q2 = ( 6.*b + 30.*c ) / 6.; | |
| 83 | 6 | q3 = ( - b - 6.*c ) / 6.; | |
| 84 | 6 | } | |
| 85 | |||
| 86 | 608 | double MitchellNetravaliFilter::f (double x) { | |
| 87 | 608 | x = fabs(x); | |
| 88 |
4/4✓ Branch 2 → 3 taken 262 times.
✓ Branch 2 → 4 taken 346 times.
✓ Branch 4 → 5 taken 296 times.
✓ Branch 4 → 6 taken 50 times.
|
608 | return (x<1) ? (p0+x*x*(p2+x*p3)) : (x<2) ? (q0+x*(q1+x*(q2+x*q3))) : 0.0; |
| 89 | } | ||
| 90 | |||
| 91 | |||
| 92 | /*********************** | ||
| 93 | *** Lanczos3 filter *** | ||
| 94 | ***********************/ | ||
| 95 |
1/2✗ Branch 3 → 4 not taken.
✓ Branch 3 → 5 taken 60 times.
|
60 | LanczosFilter::LanczosFilter(int _taps) { |
| 96 | 60 | taps = (double)clamp(_taps, 1, 100); | |
| 97 | 60 | } | |
| 98 | |||
| 99 | 193828 | double LanczosFilter::sinc(double value) { | |
| 100 |
2/2✓ Branch 2 → 3 taken 193752 times.
✓ Branch 2 → 4 taken 76 times.
|
193828 | if (value > 0.000001) { |
| 101 | 193752 | value *= M_PI; | |
| 102 | 193752 | return sin(value) / value; | |
| 103 | } else { | ||
| 104 | 76 | return 1.0; | |
| 105 | } | ||
| 106 | } | ||
| 107 | |||
| 108 | 96956 | double LanczosFilter::f(double value) { | |
| 109 | 96956 | value = fabs(value); | |
| 110 | |||
| 111 |
2/2✓ Branch 2 → 3 taken 96914 times.
✓ Branch 2 → 6 taken 42 times.
|
96956 | if (value < taps) { |
| 112 | 96914 | return (sinc(value) * sinc(value / taps)); | |
| 113 | } else { | ||
| 114 | 42 | return 0.0; | |
| 115 | } | ||
| 116 | } | ||
| 117 | |||
| 118 | |||
| 119 | /*********************** | ||
| 120 | *** Blackman filter *** | ||
| 121 | ***********************/ | ||
| 122 | ✗ | BlackmanFilter::BlackmanFilter(int _taps) { | |
| 123 | ✗ | taps = (double)clamp(_taps, 1, 100); | |
| 124 | ✗ | rtaps = 1.0/taps; | |
| 125 | ✗ | } | |
| 126 | |||
| 127 | ✗ | double BlackmanFilter::f(double value) { | |
| 128 | ✗ | value = fabs(value); | |
| 129 | |||
| 130 | ✗ | if (value < taps) { | |
| 131 | ✗ | if (value > 0.000001) { | |
| 132 | ✗ | value *= M_PI; | |
| 133 | ✗ | return (sin(value) / value) * (0.42 + 0.5*cos(value*rtaps) + 0.08*cos(2*value*rtaps)); | |
| 134 | } else { | ||
| 135 | ✗ | return 1.0; | |
| 136 | } | ||
| 137 | } else { | ||
| 138 | ✗ | return 0.0; | |
| 139 | } | ||
| 140 | } | ||
| 141 | |||
| 142 | |||
| 143 | /*********************** | ||
| 144 | *** Spline16 filter *** | ||
| 145 | ***********************/ | ||
| 146 | |||
| 147 | ✗ | double Spline16Filter::f(double value) { | |
| 148 | ✗ | value = fabs(value); | |
| 149 | |||
| 150 | ✗ | if (value < 1.0) { | |
| 151 | ✗ | return ( ( value - 9.0/5.0 ) * value - 1.0/5.0 ) * value + 1.0; | |
| 152 | ✗ | } else if (value < 2.0) { | |
| 153 | ✗ | return ( ( -1.0/3.0 * (value-1.0) + 4.0/5.0 ) * (value-1.0) - 7.0/15.0 ) * (value-1.0); | |
| 154 | } | ||
| 155 | ✗ | return 0.0; | |
| 156 | } | ||
| 157 | |||
| 158 | /*********************** | ||
| 159 | *** Spline36 filter *** | ||
| 160 | ***********************/ | ||
| 161 | |||
| 162 | 294 | double Spline36Filter::f(double value) { | |
| 163 | 294 | value = fabs(value); | |
| 164 | |||
| 165 |
2/2✓ Branch 2 → 3 taken 92 times.
✓ Branch 2 → 4 taken 202 times.
|
294 | if (value < 1.0) { |
| 166 | 92 | return ( ( 13.0/11.0 * (value ) - 453.0/ 209.0 ) * (value ) - 3.0/ 209.0 ) *(value ) + 1.0; | |
| 167 |
2/2✓ Branch 4 → 5 taken 96 times.
✓ Branch 4 → 6 taken 106 times.
|
202 | } else if (value < 2.0) { |
| 168 | 96 | return ( ( -6.0/11.0 * (value-1.0) + 270.0/ 209.0 ) * (value-1.0) - 156.0/ 209.0 ) *(value-1.0); | |
| 169 |
2/2✓ Branch 6 → 7 taken 95 times.
✓ Branch 6 → 8 taken 11 times.
|
106 | } else if (value < 3.0) { |
| 170 | 95 | return ( ( 1.0/11.0 * (value-2.0) - 45.0/ 209.0 ) * (value-2.0) + 26.0/ 209.0 ) *(value-2.0); | |
| 171 | } | ||
| 172 | 11 | return 0.0; | |
| 173 | } | ||
| 174 | |||
| 175 | /*********************** | ||
| 176 | *** Spline64 filter *** | ||
| 177 | ***********************/ | ||
| 178 | |||
| 179 | ✗ | double Spline64Filter::f(double value) { | |
| 180 | ✗ | value = fabs(value); | |
| 181 | |||
| 182 | ✗ | if (value < 1.0) { | |
| 183 | ✗ | return (( 49.0/41.0 * (value ) - 6387.0/2911.0) * (value ) - 3.0/2911.0) * (value ) + 1.0; | |
| 184 | ✗ | } else if (value < 2.0) { | |
| 185 | ✗ | return ((-24.0/41.0 * (value-1.0) + 4032.0/2911.0) * (value-1.0) - 2328.0/2911.0) * (value-1.0); | |
| 186 | ✗ | } else if (value < 3.0) { | |
| 187 | ✗ | return (( 6.0/41.0 * (value-2.0) - 1008.0/2911.0) * (value-2.0) + 582.0/2911.0) * (value-2.0); | |
| 188 | ✗ | } else if (value < 4.0) { | |
| 189 | ✗ | return ((- 1.0/41.0 * (value-3.0) + 168.0/2911.0) * (value-3.0) - 97.0/2911.0) * (value-3.0); | |
| 190 | } | ||
| 191 | ✗ | return 0.0; | |
| 192 | } | ||
| 193 | |||
| 194 | /*********************** | ||
| 195 | *** Gaussian filter *** | ||
| 196 | ***********************/ | ||
| 197 | |||
| 198 | /* Solve taps from p*value*value < 9 as pow(2.0, -9.0) == 1.0/512.0 i.e 0.5 bit | ||
| 199 | value*value < 9/p p = param*0.1; | ||
| 200 | value*value < 90/param | ||
| 201 | value*value < 90/{0.1, 22.5, 30.0, 100.0} | ||
| 202 | value*value < {900, 4.0, 3.0, 0.9} | ||
| 203 | value < {30, 2.0, 1.73, 0.949} */ | ||
| 204 | |||
| 205 | ✗ | GaussianFilter::GaussianFilter(double p, double _b, double _s) { | |
| 206 | ✗ | param = clamp(p, 0.01, 100.0); | |
| 207 | ✗ | b = clamp(_b, 1.5, 3.5); | |
| 208 | ✗ | s = _s; | |
| 209 | ✗ | if (_s == 0) // auto-support signal | |
| 210 | { | ||
| 211 | // get support from b and param for 0.01 of resudual kernel value | ||
| 212 | // equatiion is s = sqrt(-ln(0.01)/(param*ln(b)) | ||
| 213 | // where ln(0.01) is about -4.6 and -ln(0.01) is 4.6 | ||
| 214 | ✗ | s = sqrt(4.6 / ((param * 0.1) * log(b))); | |
| 215 | } | ||
| 216 | ✗ | s = clamp(s, 0.1, 150.0); | |
| 217 | ✗ | } | |
| 218 | |||
| 219 | ✗ | double GaussianFilter::f(double value) { | |
| 220 | ✗ | double p = param * 0.1; | |
| 221 | ✗ | return pow(b, -p * value * value); // <3.7.4: b was fixed at 2.0 | |
| 222 | } | ||
| 223 | |||
| 224 | /*********************** | ||
| 225 | *** Sinc filter *** | ||
| 226 | ***********************/ | ||
| 227 |
1/2✗ Branch 3 → 4 not taken.
✓ Branch 3 → 5 taken 3 times.
|
3 | SincFilter::SincFilter(int _taps) { |
| 228 | 3 | taps = (double)clamp(_taps, 1, 150); | |
| 229 | 3 | } | |
| 230 | |||
| 231 | 524 | double SincFilter::f(double value) { | |
| 232 | 524 | value = fabs(value); | |
| 233 | |||
| 234 |
2/2✓ Branch 2 → 3 taken 506 times.
✓ Branch 2 → 4 taken 18 times.
|
524 | if (value > 0.000001) { |
| 235 | 506 | value *= M_PI; | |
| 236 | 506 | return sin(value)/value; | |
| 237 | } else { | ||
| 238 | 18 | return 1.0; | |
| 239 | } | ||
| 240 | } | ||
| 241 | |||
| 242 | |||
| 243 | /********************** | ||
| 244 | *** SinPower filter *** | ||
| 245 | ***********************/ | ||
| 246 | |||
| 247 | ✗ | SinPowerFilter::SinPowerFilter(double p) { | |
| 248 | ✗ | param = clamp(p, 1.0, 10.0); | |
| 249 | ✗ | } | |
| 250 | |||
| 251 | ✗ | double SinPowerFilter::f(double value) { | |
| 252 | ✗ | value = fabs(value); | |
| 253 | ✗ | value *= M_PI / param; | |
| 254 | |||
| 255 | ✗ | if (value < (M_PI / 2)) return pow(cos(value), 1.8); | |
| 256 | else | ||
| 257 | { | ||
| 258 | ✗ | if (value < M_PI) return -(cos(value) * cos(value)) / (0.9 * value); | |
| 259 | ✗ | else return 0; | |
| 260 | } | ||
| 261 | } | ||
| 262 | |||
| 263 | /*********************** | ||
| 264 | *** SincLin2 filter *** | ||
| 265 | ***********************/ | ||
| 266 | |||
| 267 | ✗ | SincLin2Filter::SincLin2Filter(int _taps) | |
| 268 | { | ||
| 269 | ✗ | taps = (double)clamp(_taps, 1, 30); | |
| 270 | ✗ | } | |
| 271 | |||
| 272 | ✗ | double SincLin2Filter::sinc(double value) | |
| 273 | { | ||
| 274 | ✗ | if (value > 0.000001) | |
| 275 | { | ||
| 276 | ✗ | value *= M_PI; | |
| 277 | ✗ | return sin(value) / value; | |
| 278 | } | ||
| 279 | ✗ | else return 1.0; | |
| 280 | } | ||
| 281 | |||
| 282 | ✗ | double SincLin2Filter::f(double value) | |
| 283 | { | ||
| 284 | ✗ | value = fabs(value); | |
| 285 | |||
| 286 | ✗ | if (value < (taps / 2.0)) return sinc(value); | |
| 287 | ✗ | else return sinc(value) * ((2.0 - (2.0 * value / taps))); | |
| 288 | |||
| 289 | } | ||
| 290 | |||
| 291 | |||
| 292 | /********************************* | ||
| 293 | *** UserDefined2 filter *** | ||
| 294 | *********************************/ | ||
| 295 | |||
| 296 | ✗ | UserDefined2Filter::UserDefined2Filter(double _b, double _c, double _s) | |
| 297 | { | ||
| 298 | ✗ | a = 1.0; // 0 sample = 1 | |
| 299 | ✗ | b = (double)clamp(_b, -50.0, 250.0); // 1 and -1 sample | |
| 300 | ✗ | c = (double)clamp(_c, -50.0, 250.0); // 2 and -2 sample | |
| 301 | ✗ | b = (b - 16.0) / 219.0; | |
| 302 | ✗ | c = (c - 16.0) / 219.0; | |
| 303 | ✗ | s = (double)clamp(_s, 1.5, 15.0); // filter support for resampler | |
| 304 | ✗ | } | |
| 305 | |||
| 306 | ✗ | double UserDefined2Filter::sinc(double value) | |
| 307 | { | ||
| 308 | |||
| 309 | ✗ | if (fabs(value) > 0.000001) | |
| 310 | { | ||
| 311 | ✗ | value *= M_PI; | |
| 312 | ✗ | return sin(value) / value; | |
| 313 | } | ||
| 314 | ✗ | else return 1.0; | |
| 315 | } | ||
| 316 | |||
| 317 | ✗ | double UserDefined2Filter::f(double x) | |
| 318 | { | ||
| 319 | ✗ | x = fabs(x); | |
| 320 | |||
| 321 | ✗ | return c * sinc(x + 2) + b * sinc(x + 1) + a * sinc(x) + b * sinc(x - 1) + c * sinc(x - 2); | |
| 322 | } | ||
| 323 | |||
| 324 | /* | ||
| 325 | * OPTIMAL SCANLINE CALCULATION NOTES (L2 CACHE BLOCKING) | ||
| 326 | * | ||
| 327 | * This function calculates the optimal vertical strip size (max_scanline) | ||
| 328 | * to be processed in a cache-blocked horizontal resizing operation. | ||
| 329 | * | ||
| 330 | * CONTEXT: Single-threaded, high-throughput workload with private L2 cache. | ||
| 331 | * The high FPS target justifies a more aggressive cache reservation factor. | ||
| 332 | * | ||
| 333 | * 1. COEFFICIENT TABLE EXCLUSION (Horizontal 2x resize of fullhd content): | ||
| 334 | * The coefficient table (~245 KB) is excluded from the calculation. It is | ||
| 335 | * treated as a streamed resource due to its size relative to L2 (512 KB). | ||
| 336 | * The hardware prefetcher is expected to handle its sequential access pattern | ||
| 337 | * efficiently without residing fully in the reserved L2 working set. | ||
| 338 | * | ||
| 339 | * 2. CACHE RESERVATION HEURISTIC: | ||
| 340 | * A factor of 0.75 (3/4) is reserved for the working set. This is an aggressive | ||
| 341 | * approach for a single-threaded task with a private L2 cache, minimizing | ||
| 342 | * cache thrashing risk from OS context switches while maximizing block size. | ||
| 343 | * | ||
| 344 | * 3. CONCRETE SCENARIO (512 KB L2, 1920->3840 upscale): | ||
| 345 | * Reserved L2 space: 512 KB * 0.75 = 384 KB (393,216 bytes). | ||
| 346 | * One scanline strip (src + tgt) is 23,040 bytes. | ||
| 347 | * The resulting max_scanline is 17 (393,216 / 23,040 ~ 17.06). | ||
| 348 | */ | ||
| 349 | 290 | int ResamplingProgram::resampler_h_detect_optimal_scanline(int src_width, int tgt_width, size_t l2_cache_size_bytes, size_t pixel_size) { | |
| 350 | 290 | constexpr double CACHE_RESERVE_FACTOR = 0.75; | |
| 351 | |||
| 352 | // Calculate the bytes needed for one (Source + Destination) scanline strip | ||
| 353 | 290 | size_t scanline_bytes = (static_cast<size_t>(src_width) + static_cast<size_t>(tgt_width)) * pixel_size; | |
| 354 | |||
| 355 | // Calculate the reserved bytes based on the aggressive factor | ||
| 356 | // Use floating point math for precision, then cast to size_t | ||
| 357 | 290 | size_t reserved_l2_bytes = static_cast<size_t>( | |
| 358 | 290 | static_cast<double>(l2_cache_size_bytes) * CACHE_RESERVE_FACTOR | |
| 359 | ); | ||
| 360 | |||
| 361 | // Calculate max_scanline (integer division for floor) | ||
| 362 | 290 | int max_scanline = static_cast<int>(reserved_l2_bytes / scanline_bytes); | |
| 363 | |||
| 364 | // Clamp to practical bounds (4 to 64 is typical range for strip size) | ||
| 365 | // Dynamic by sample_size. For float32 (size=4) was 4 min and 64 max, so for uint8_t size=1 it will be 4x times more. | ||
| 366 | // Heuristic limits to avoid too small or too large strip sizes. | ||
| 367 | int iMinLimit; | ||
| 368 | int iMaxLimit; | ||
| 369 | |||
| 370 |
3/4✓ Branch 2 → 3 taken 25 times.
✓ Branch 2 → 4 taken 66 times.
✓ Branch 2 → 5 taken 199 times.
✗ Branch 2 → 6 not taken.
|
290 | switch (pixel_size) |
| 371 | { | ||
| 372 | 25 | case 4: // float | |
| 373 | 25 | iMinLimit = 4; | |
| 374 | 25 | iMaxLimit = 64; | |
| 375 | 25 | break; | |
| 376 | 66 | case 2: // uint16_t | |
| 377 | 66 | iMinLimit = 8; | |
| 378 | 66 | iMaxLimit = 128; | |
| 379 | 66 | break; | |
| 380 | 199 | case 1: // uint8_t | |
| 381 | 199 | iMinLimit = 16; | |
| 382 | 199 | iMaxLimit = 256; | |
| 383 | 199 | break; | |
| 384 | ✗ | default: // should never happen | |
| 385 | ✗ | iMinLimit = 4; | |
| 386 | ✗ | iMaxLimit = 64; | |
| 387 | ✗ | break; | |
| 388 | } | ||
| 389 | |||
| 390 | 290 | max_scanline = std::min(std::max(max_scanline, iMinLimit), iMaxLimit); | |
| 391 | |||
| 392 | 290 | return max_scanline; | |
| 393 | } | ||
| 394 | |||
| 395 | /** | ||
| 396 | * @brief Checks if the data access pattern for horizontal resampling requires the slower Transpose method, | ||
| 397 | * or if the faster Permutex approach can be used. | ||
| 398 | * | ||
| 399 | * This function determines the feasibility of using vector instructions (like AVX2/AVX-512 Permutex) | ||
| 400 | * to load the source samples required for a group of output samples. This method is preferred when | ||
| 401 | * the total spread of necessary source samples is small. | ||
| 402 | * | ||
| 403 | * @note This check is performed for methods that process multiple output sample groups per loop pass | ||
| 404 | * (e.g., loading 2 groups of 8 samples for AVX-512). | ||
| 405 | * | ||
| 406 | * @param iSamplesInTheGroup The number of output samples processed in one vector iteration (e.g., 8, 16, 64). | ||
| 407 | * (Usually referred to as PIXELS_AT_A_TIME). | ||
| 408 | * @param permutex_index_diff_limit The maximum byte/element difference allowed between the earliest and latest | ||
| 409 | * required source sample for Permutex to be viable. This limit is dictated | ||
| 410 | * by the specific Permutex intrinsic used (e.g., 8 for _mm256_permutevar8x32_ps). | ||
| 411 | * @param kernel_size The size of the resampling filter kernel (number of coefficients). | ||
| 412 | * @return true If the source sample spread is too large, meaning only the Transpose method is allowed. | ||
| 413 | * @return false If the source sample spread is small enough, meaning the Permutex method can be used. | ||
| 414 | */ | ||
| 415 | 91 | bool ResamplingProgram::resize_h_planar_gather_permutex_vstripe_check(int iSamplesInTheGroup, int permutex_index_diff_limit, int kernel_size) | |
| 416 | { | ||
| 417 | // iSamplesInTheGroup is usually denoted as PIXELS_AT_A_TIME in H resampler code | ||
| 418 | // permutex_index_diff_limit is like iAccessibleSourceSamplesToGroup | ||
| 419 | |||
| 420 | // Alignment checks ensure safe access to pre-calculated arrays for the entire vector block. | ||
| 421 | |||
| 422 | // 'target_size_alignment' ensures safe access for the entire group (x to x + iSamplesInTheGroup - 1) | ||
| 423 | // Example: for iSamplesInTheGroup = 8, we need to ensure that | ||
| 424 | // - program->pixel_offset[x + 0] to program->pixel_offset[x + 7] and | ||
| 425 | // - corresponding coefficients (spread over coeff-strides like current_coeff + filter_size*0 to filter_size*7), where filter_size is the aligned coefficient stride. | ||
| 426 | // are valid if iSamplesInTheGroup = 8. | ||
| 427 |
1/2✗ Branch 2 → 3 not taken.
✓ Branch 2 → 4 taken 91 times.
|
91 | assert(target_size_alignment >= iSamplesInTheGroup); |
| 428 | |||
| 429 | // Ensure that coefficient loading is safe for "kernel_size" element loads. | ||
| 430 | // But this is not true. For float pixel types, we can keep the alignment to 8, which can be less than kernel_size. | ||
| 431 | // It's because it depends on the filter's implementation. E.g. resize_h_planar_float_avx512_permutex_vstripe_ks16 has kernel_size=16. | ||
| 432 | // However, it uses gather loads from coeffs, which do not need to be aligned to kernel_size. | ||
| 433 | // uint8_t avx512 versions with ks16 do need 32 byte alignment, the 'short' coefficients stride is aligned to 32 bytes, that is 16 coeffs. | ||
| 434 | // So in this case, we need the alignment. | ||
| 435 | // The check cannot be generalized here, it is put in the specific resampler implementations if needed. | ||
| 436 | // assert(filter_size_alignment >= kernel_size); | ||
| 437 | |||
| 438 |
2/2✓ Branch 10 → 5 taken 91 times.
✓ Branch 10 → 11 taken 90 times.
|
181 | for (int x = 0; x < target_size; x += iSamplesInTheGroup) // check each group |
| 439 | { | ||
| 440 | // Get the index of the first required source sample for this group. | ||
| 441 | 91 | int start_off = pixel_offset[x + 0]; | |
| 442 | |||
| 443 | // Get the index of the last required source sample. This is the offset for the last | ||
| 444 | // output pixel in the group (x + iSamplesInTheGroup - 1) plus the last kernel tap (kernel_size - 1). | ||
| 445 | // Note: pixel_offset[] values for x >= target_size are pre-padded to match target_size-1 (ensured by resize_prepare_coeffs()). | ||
| 446 | 91 | const int end_off = pixel_offset[x + (iSamplesInTheGroup - 1)] + (kernel_size - 1); | |
| 447 | |||
| 448 | // Check the total spread (difference) in source sample indices. | ||
| 449 | // This difference must be less than the limit imposed by the Permutex intrinsic's addressing capability. | ||
| 450 | // Examples of permutex_index_diff_limit: | ||
| 451 | // - 8 for _mm256_permutevar8x32_ps (float, avx2) | ||
| 452 | // - 32 for _mm512_permutex2var_ps (float, avx512) | ||
| 453 | // - 64 for _mm512_permutex2var_epi16 (uint16_t, avx512) | ||
| 454 | // - 128 for _mm512_permutex2var_epi8 (uint8_t, avx512) | ||
| 455 |
2/2✓ Branch 7 → 8 taken 1 time.
✓ Branch 7 → 9 taken 90 times.
|
91 | if ((end_off - start_off) >= permutex_index_diff_limit) { |
| 456 | 1 | return true; // spread is too wide; only the transpose method is allowed. | |
| 457 | } | ||
| 458 | } | ||
| 459 | 90 | return false; // Spread is acceptable; Permutex is OK. | |
| 460 | } | ||
| 461 | |||
| 462 | |||
| 463 | /****************************** | ||
| 464 | **** Resampling Patterns **** | ||
| 465 | *****************************/ | ||
| 466 | |||
| 467 | 290 | ResamplingProgram* ResamplingFunction::GetResamplingProgram(int source_size, double crop_start, double crop_size, int target_size, int bits_per_pixel, | |
| 468 | double center_pos_src, double center_pos_dst, | ||
| 469 | IScriptEnvironment* env) | ||
| 470 | { | ||
| 471 | // edge condition ideas from fmtconv, thanks. | ||
| 472 | 290 | double src_step = crop_size / double(target_size); // Distance between source pixels for adjacent dest pixels | |
| 473 | 290 | double zc_size = std::max(src_step, 1.0) / 1.0; // Size of filter unit step (kernel_scale=1.0 in our case) | |
| 474 | 290 | double imp_step = 1.0 / zc_size; // Corresponding distance in the impulse | |
| 475 |
1/2✓ Branch 3 → 4 taken 290 times.
✗ Branch 3 → 78 not taken.
|
290 | double filter_support = support() * zc_size; // Number of source pixels covered by the FIR |
| 476 | |||
| 477 | 290 | int fir_filter_size = std::max(int(std::ceil(filter_support * 2)), 1); | |
| 478 | 290 | int max_kernel_size = 0; | |
| 479 | |||
| 480 | 290 | const int last_line = source_size - 1; | |
| 481 | |||
| 482 |
3/8✓ Branch 5 → 6 taken 290 times.
✗ Branch 5 → 78 not taken.
✓ Branch 6 → 7 taken 290 times.
✗ Branch 6 → 72 not taken.
✗ Branch 7 → 8 not taken.
✓ Branch 7 → 9 taken 290 times.
✗ Branch 72 → 73 not taken.
✗ Branch 72 → 74 not taken.
|
290 | ResamplingProgram* program = new ResamplingProgram(fir_filter_size, source_size, target_size, crop_start, crop_size, bits_per_pixel, env); |
| 483 | |||
| 484 | // Initial position calculation | ||
| 485 | |||
| 486 | 290 | double pos = crop_start; | |
| 487 | |||
| 488 | /* | ||
| 489 | pre 3.7.4 logic: | ||
| 490 | |||
| 491 | Now in 2025, let's fact-check this comment. | ||
| 492 | |||
| 493 | pos = crop_start + ((crop_size - target_size) / (target_size*2)); // TODO this look wrong, gotta check | ||
| 494 | ==> | ||
| 495 | pos = crop_start + 1/2 * (crop_size / target_size - 1) | ||
| 496 | ==> | ||
| 497 | pos = crop_start + src_step * 0.5 - 1 * 0.5 | ||
| 498 | |||
| 499 | fmtconv generic formula: | ||
| 500 | |||
| 501 | pos = crop_start + src_step * center_pos_dst - 1 * center_pos_src; // 3.7.4- fmtconv | ||
| 502 | |||
| 503 | Solved: center_pos_dst = 0.5, center_pos_src = 0.5 in old Avisynth | ||
| 504 | |||
| 505 | */ | ||
| 506 | |||
| 507 | // Introduces an offset because samples are located at the center of the | ||
| 508 | // pixels, not on their boundaries. Excepted for pointresize. | ||
| 509 |
2/2✓ Branch 9 → 10 taken 224 times.
✓ Branch 9 → 11 taken 66 times.
|
290 | if (filter_support > 0) |
| 510 | { | ||
| 511 | // Pre 3.7.4 Avisynth worked with fixed center_pos_dst = center_pos_src = 0.5 | ||
| 512 | // Now it's externally configurable. In our use case they are always the same. | ||
| 513 | 224 | pos += src_step * center_pos_dst - 1 * center_pos_src; | |
| 514 | } | ||
| 515 | else | ||
| 516 | { | ||
| 517 | // In case of PointResize(), which now returns real 0 for support(). | ||
| 518 | // Avisynth heritage. | ||
| 519 | 66 | filter_support = 0.0001; | |
| 520 | } | ||
| 521 | |||
| 522 |
4/4✓ Branch 12 → 13 taken 91 times.
✓ Branch 12 → 15 taken 199 times.
✓ Branch 13 → 14 taken 66 times.
✓ Branch 13 → 15 taken 25 times.
|
290 | const int current_FPScale = (bits_per_pixel > 8 && bits_per_pixel <= 16) ? FPScale16 : FPScale; |
| 523 | |||
| 524 | 290 | std::vector<double> coef_tmp; | |
| 525 |
2/2✓ Branch 68 → 18 taken 10931 times.
✓ Branch 68 → 69 taken 290 times.
|
11221 | for (int i = 0; i < target_size; ++i) { |
| 526 | 10931 | coef_tmp.clear(); | |
| 527 | |||
| 528 |
1/2✗ Branch 19 → 20 not taken.
✓ Branch 19 → 21 taken 10931 times.
|
10931 | int start_pos = (int)(pos + filter_support) - fir_filter_size + 1; |
| 529 | 10931 | program->pixel_offset[i] = clamp(start_pos, 0, last_line); | |
| 530 | |||
| 531 | // First pass: Accumulate all coefficients for weighting | ||
| 532 | 10931 | double total = 0.0; | |
| 533 |
2/2✓ Branch 32 → 29 taken 107474 times.
✓ Branch 32 → 33 taken 10931 times.
|
118405 | for (int k = 0; k < fir_filter_size; ++k) { |
| 534 | 107474 | const int p = start_pos + k; | |
| 535 |
1/2✓ Branch 29 → 30 taken 107474 times.
✗ Branch 29 → 75 not taken.
|
107474 | double val = f((pos - p) * imp_step); |
| 536 |
1/2✓ Branch 30 → 31 taken 107474 times.
✗ Branch 30 → 75 not taken.
|
107474 | coef_tmp.push_back(val); |
| 537 | 107474 | total += val; | |
| 538 | } | ||
| 539 | |||
| 540 |
1/2✗ Branch 33 → 34 not taken.
✓ Branch 33 → 35 taken 10931 times.
|
10931 | if (total == 0.0) { |
| 541 | // Shouldn't happen for valid positions. | ||
| 542 | ✗ | total = 1.0; | |
| 543 | } | ||
| 544 | |||
| 545 | 10931 | const int coeff_arr_base_index = i * fir_filter_size; | |
| 546 | |||
| 547 | // Second pass: Generate real coefficients, handling edge conditions | ||
| 548 | 10931 | double accu = 0.0; | |
| 549 | 10931 | double prev_value = 0.0; | |
| 550 | |||
| 551 | 10931 | int kernel_size = 0; | |
| 552 | |||
| 553 |
2/2✓ Branch 35 → 36 taken 1130 times.
✓ Branch 35 → 43 taken 9801 times.
|
10931 | if (bits_per_pixel == 32) { |
| 554 | // Float version | ||
| 555 |
2/2✓ Branch 42 → 37 taken 8948 times.
✓ Branch 42 → 50 taken 1130 times.
|
10078 | for (int k = 0; k < fir_filter_size; ++k) { |
| 556 | 8948 | const int p = start_pos + k; | |
| 557 | 8948 | double val = coef_tmp[k]; | |
| 558 | 8948 | accu += val; | |
| 559 |
4/4✓ Branch 38 → 39 taken 8789 times.
✓ Branch 38 → 41 taken 159 times.
✓ Branch 39 → 40 taken 8627 times.
✓ Branch 39 → 41 taken 162 times.
|
8948 | if (p >= 0 && p <= last_line) { |
| 560 | 8627 | program->pixel_coefficient_float[coeff_arr_base_index + kernel_size] = float(accu / total); | |
| 561 | 8627 | ++kernel_size; | |
| 562 | 8627 | accu = 0; | |
| 563 | } | ||
| 564 | } | ||
| 565 | } | ||
| 566 | else { | ||
| 567 | // Integer version - using upscaled integer arithmetic (FPScale/FPScale16) | ||
| 568 |
2/2✓ Branch 49 → 44 taken 98526 times.
✓ Branch 49 → 50 taken 9801 times.
|
108327 | for (int k = 0; k < fir_filter_size; ++k) { |
| 569 | 98526 | const int p = start_pos + k; | |
| 570 | 98526 | double val = coef_tmp[k]; | |
| 571 | 98526 | accu += val; | |
| 572 |
4/4✓ Branch 45 → 46 taken 96183 times.
✓ Branch 45 → 48 taken 2343 times.
✓ Branch 46 → 47 taken 93726 times.
✓ Branch 46 → 48 taken 2457 times.
|
98526 | if (p >= 0 && p <= last_line) { |
| 573 | 93726 | double new_value = prev_value + accu / total; | |
| 574 | // differential approach ensures the filter coefficients sum to exactly FPScale) | ||
| 575 | // The subtraction method guarantees that no matter how many terms we add, the | ||
| 576 | // final sum will be exactly equal to the fixed-point representation of 1.0. | ||
| 577 | 93726 | program->pixel_coefficient[coeff_arr_base_index + kernel_size] = (short)((int)(new_value * current_FPScale + 0.5) - int(prev_value * current_FPScale + 0.5)); | |
| 578 | 93726 | prev_value = new_value; | |
| 579 | 93726 | ++kernel_size; | |
| 580 | 93726 | accu = 0; | |
| 581 | } | ||
| 582 | } | ||
| 583 | } | ||
| 584 | |||
| 585 | // We even haven't reached any valid line, | ||
| 586 | // or gathered accu values from past last line. | ||
| 587 |
2/2✓ Branch 50 → 51 taken 701 times.
✓ Branch 50 → 59 taken 10230 times.
|
10931 | if (accu != 0) |
| 588 | { | ||
| 589 |
2/2✓ Branch 51 → 52 taken 691 times.
✓ Branch 51 → 55 taken 10 times.
|
701 | if (kernel_size > 0) { |
| 590 | // Assign the remaining accumulator to the last line, just like we put | ||
| 591 | // the accumulator before the first valid line to the first line. | ||
| 592 |
2/2✓ Branch 52 → 53 taken 38 times.
✓ Branch 52 → 54 taken 653 times.
|
691 | if (bits_per_pixel == 32) |
| 593 | 38 | program->pixel_coefficient_float[coeff_arr_base_index + kernel_size - 1] += float(accu / total); | |
| 594 | else { | ||
| 595 | 653 | double new_value = prev_value + accu / total; | |
| 596 | 653 | program->pixel_coefficient[coeff_arr_base_index + kernel_size - 1] += (short)((int)(new_value * current_FPScale + 0.5) - int(prev_value * current_FPScale + 0.5)); | |
| 597 | } | ||
| 598 | // no change in kernel_size | ||
| 599 | } | ||
| 600 | else | ||
| 601 | { | ||
| 602 | // new entry, accu/total must be 1.0 here (we always normalize) | ||
| 603 |
1/2✗ Branch 55 → 56 not taken.
✓ Branch 55 → 57 taken 10 times.
|
10 | if (bits_per_pixel == 32) |
| 604 | ✗ | program->pixel_coefficient_float[coeff_arr_base_index + kernel_size] = float(accu / total); | |
| 605 | else | ||
| 606 | 10 | program->pixel_coefficient[coeff_arr_base_index + kernel_size] = (short)((int)(accu / total * current_FPScale + 0.5)); | |
| 607 | 10 | ++kernel_size; | |
| 608 | } | ||
| 609 | } | ||
| 610 | |||
| 611 |
1/2✗ Branch 59 → 60 not taken.
✓ Branch 59 → 64 taken 10931 times.
|
10931 | if (kernel_size == 0) { |
| 612 | // write a single 1.0 coeff entry | ||
| 613 | ✗ | if (bits_per_pixel == 32) | |
| 614 | ✗ | program->pixel_coefficient_float[coeff_arr_base_index + kernel_size] = 1.0f; | |
| 615 | else | ||
| 616 | ✗ | program->pixel_coefficient[coeff_arr_base_index + kernel_size] = (short)((int)(1.0 * current_FPScale + 0.5)); | |
| 617 | ✗ | ++kernel_size; | |
| 618 | } | ||
| 619 | |||
| 620 | 10931 | program->kernel_sizes[i] = kernel_size; | |
| 621 |
2/2✓ Branch 65 → 66 taken 679 times.
✓ Branch 65 → 67 taken 10252 times.
|
10931 | if (kernel_size > max_kernel_size) max_kernel_size = kernel_size; |
| 622 | |||
| 623 | 10931 | pos += src_step; | |
| 624 | } | ||
| 625 | |||
| 626 | // the different kernel sizes and coeff table will be later postprocessed | ||
| 627 | // to have aligned and equally sized coefficients. | ||
| 628 | |||
| 629 | 290 | program->filter_size_real = max_kernel_size; | |
| 630 | // can be less than original filter size if source dimensions are small | ||
| 631 | |||
| 632 | 290 | return program; | |
| 633 | 290 | } | |
| 634 | |||
| 635 | |||
| 636 | #if 0 | ||
| 637 | // old pre 3.7.4, kept for reference | ||
| 638 | ResamplingProgram* ResamplingFunction::GetResamplingProgram(int source_size, double crop_start, double crop_size, int target_size, int bits_per_pixel, IScriptEnvironment* env) | ||
| 639 | { | ||
| 640 | double filter_scale = double(target_size) / crop_size; | ||
| 641 | double filter_step = min(filter_scale, 1.0); | ||
| 642 | double filter_support = support() / filter_step; | ||
| 643 | int fir_filter_size = int(ceil(filter_support*2)); | ||
| 644 | |||
| 645 | ResamplingProgram* program = new ResamplingProgram(fir_filter_size, source_size, target_size, crop_start, crop_size, bits_per_pixel, env); | ||
| 646 | |||
| 647 | // this variable translates such that the image center remains fixed | ||
| 648 | double pos; | ||
| 649 | double pos_step = crop_size / target_size; | ||
| 650 | |||
| 651 | if (source_size <= filter_support) { | ||
| 652 | env->ThrowError("Resize: Source image too small for this resize method. Width=%d, Support=%d", source_size, int(ceil(filter_support))); | ||
| 653 | } | ||
| 654 | |||
| 655 | if (fir_filter_size == 1) // PointResize | ||
| 656 | pos = crop_start; | ||
| 657 | else | ||
| 658 | pos = crop_start + ((crop_size - target_size) / (target_size*2)); // TODO this look wrong, gotta check | ||
| 659 | |||
| 660 | const int current_FPScale = (bits_per_pixel > 8 && bits_per_pixel <= 16) ? FPScale16 : FPScale; | ||
| 661 | |||
| 662 | for (int i = 0; i < target_size; ++i) { | ||
| 663 | // Clamp start and end position such that it does not exceed frame size | ||
| 664 | int end_pos = int(pos + filter_support); | ||
| 665 | |||
| 666 | if (end_pos > source_size-1) | ||
| 667 | end_pos = source_size-1; | ||
| 668 | |||
| 669 | int start_pos = end_pos - fir_filter_size + 1; | ||
| 670 | |||
| 671 | if (start_pos < 0) | ||
| 672 | start_pos = 0; | ||
| 673 | |||
| 674 | program->pixel_offset[i] = start_pos; | ||
| 675 | |||
| 676 | // check the simd-optimized (8 pixels and filter coefficients at a time) limit to not reach beyond the last pixel | ||
| 677 | // in order not to have NaN floats | ||
| 678 | if (start_pos + AlignNumber(fir_filter_size, ALIGN_FLOAT_RESIZER_COEFF_SIZE) - 1 > source_size - 1) | ||
| 679 | { | ||
| 680 | if (!program->overread_possible_filter_size_aligned) { | ||
| 681 | // register the first occurance | ||
| 682 | program->overread_possible_filter_size_aligned = true; | ||
| 683 | program->source_overread_offset = start_pos; | ||
| 684 | program->source_overread_beyond_targetx = i; | ||
| 685 | } | ||
| 686 | } | ||
| 687 | |||
| 688 | // the following code ensures that the coefficients add to exactly FPScale or FPScale16 | ||
| 689 | double total = 0.0; | ||
| 690 | |||
| 691 | // Ensure that we have a valid position | ||
| 692 | double ok_pos = clamp(pos, 0.0, (double)(source_size-1)); | ||
| 693 | |||
| 694 | // Accumulate all coefficients for weighting | ||
| 695 | for (int j = 0; j < fir_filter_size; ++j) { | ||
| 696 | total += f((start_pos+j - ok_pos) * filter_step); | ||
| 697 | } | ||
| 698 | |||
| 699 | if (total == 0.0) { | ||
| 700 | // Shouldn't happened for valid positions. | ||
| 701 | total = 1.0; | ||
| 702 | } | ||
| 703 | |||
| 704 | double value = 0.0; | ||
| 705 | |||
| 706 | // Now we generate real coefficient | ||
| 707 | if (bits_per_pixel == 32) { | ||
| 708 | // float | ||
| 709 | for (int k = 0; k < fir_filter_size; ++k) { | ||
| 710 | double new_value = value + f((start_pos + k - ok_pos) * filter_step) / total; | ||
| 711 | program->pixel_coefficient_float[i*fir_filter_size + k] = float(new_value - value); // no scaling for float | ||
| 712 | value = new_value; | ||
| 713 | } | ||
| 714 | } | ||
| 715 | else { | ||
| 716 | for (int k = 0; k < fir_filter_size; ++k) { | ||
| 717 | double new_value = value + f((start_pos + k - ok_pos) * filter_step) / total; | ||
| 718 | // FIXME: is it correct to round negative values upwards? | ||
| 719 | // Answer : No with int cast, yes with floor() instead. | ||
| 720 | program->pixel_coefficient[i*fir_filter_size + k] = (short)((int)(new_value*current_FPScale + 0.5) - int(value*current_FPScale + 0.5)); // to make it round across pixels | ||
| 721 | value = new_value; | ||
| 722 | } | ||
| 723 | } | ||
| 724 | |||
| 725 | pos += pos_step; | ||
| 726 | } | ||
| 727 | |||
| 728 | // aligned as 8, now fill with safe values for 8 pixels/cycle simd loop | ||
| 729 | for (int i = target_size; i < AlignNumber(target_size, ALIGN_RESIZER_TARGET_SIZE); ++i) | ||
| 730 | program->pixel_offset[i] = source_size - fir_filter_size; | ||
| 731 | |||
| 732 | return program; | ||
| 733 | } | ||
| 734 | #endif | ||
| 735 |