|
| 1 | +// |
| 2 | +// ScaleInterpolator.cpp |
| 3 | +// JxclCoder [https://github.com/awxkee/jxl-coder-swift] |
| 4 | +// |
| 5 | +// Created by Radzivon Bartoshyk on 03/10/2023. |
| 6 | +// |
| 7 | +// Permission is hereby granted, free of charge, to any person obtaining a copy |
| 8 | +// of this software and associated documentation files (the "Software"), to deal |
| 9 | +// in the Software without restriction, including without limitation the rights |
| 10 | +// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell |
| 11 | +// copies of the Software, and to permit persons to whom the Software is |
| 12 | +// furnished to do so, subject to the following conditions: |
| 13 | +// |
| 14 | +// The above copyright notice and this permission notice shall be included in |
| 15 | +// all copies or substantial portions of the Software. |
| 16 | +// |
| 17 | +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR |
| 18 | +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, |
| 19 | +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE |
| 20 | +// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER |
| 21 | +// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, |
| 22 | +// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN |
| 23 | +// THE SOFTWARE. |
| 24 | +// |
| 25 | + |
| 26 | +#include <stdio.h> |
| 27 | +#include "ScaleInterpolator.h" |
| 28 | +#import "half.hpp" |
| 29 | +#include <algorithm> |
| 30 | + |
| 31 | +using namespace half_float; |
| 32 | +using namespace std; |
| 33 | + |
| 34 | +// P Found using maxima |
| 35 | +// |
| 36 | +// y(x) := 4 * x * (%pi-x) / (%pi^2) ; |
| 37 | +// z(x) := (1-p)*y(x) + p * y(x)^2; |
| 38 | +// e(x) := z(x) - sin(x); |
| 39 | +// solve( diff( integrate( e(x)^2, x, 0, %pi/2 ), p ) = 0, p ),numer; |
| 40 | +// |
| 41 | +// [p = .2248391013559941] |
| 42 | +template <typename T> |
| 43 | +inline T fastSin1(T x) |
| 44 | +{ |
| 45 | + // x = fmod(x + M_PI, M_PI * 2) - M_PI; |
| 46 | + constexpr T A = T(4.0)/(T(M_PI)*T(M_PI)); |
| 47 | + constexpr T P = 0.2248391013559941; |
| 48 | + T y = A * x * ( T(M_PI) - x ); |
| 49 | + return y * ( (1-P) + y * P ); |
| 50 | +} |
| 51 | + |
| 52 | +// P and Q found using maxima |
| 53 | +// |
| 54 | +// y(x) := 4 * x * (%pi-x) / (%pi^2) ; |
| 55 | +// zz(x) := (1-p-q)*y(x) + p * y(x)^2 + q * y(x)^3 |
| 56 | +// ee(x) := zz(x) - sin(x) |
| 57 | +// solve( [ integrate( diff(ee(x)^2, p ), x, 0, %pi/2 ) = 0, integrate( diff(ee(x)^2,q), x, 0, %pi/2 ) = 0 ] , [p,q] ),numer; |
| 58 | +// |
| 59 | +// [[p = .1952403377008734, q = .01915214119105392]] |
| 60 | +template <typename T> |
| 61 | +inline T fastSin2(T x) |
| 62 | +{ |
| 63 | + constexpr T A = T(4.0)/(T(M_PI)*T(M_PI)); |
| 64 | + constexpr T P = 0.1952403377008734; |
| 65 | + constexpr T Q = 0.01915214119105392; |
| 66 | + |
| 67 | + T y = A * x * ( T(M_PI) - x ); |
| 68 | + |
| 69 | + return y*((1-P-Q) + y*( P + y * Q ) ) ; |
| 70 | +} |
| 71 | + |
| 72 | +template <typename T> |
| 73 | +inline T fastCos(T x) { |
| 74 | + constexpr T C0 = 0.99940307; |
| 75 | + constexpr T C1 = -0.49558072; |
| 76 | + constexpr T C2 = 0.03679168; |
| 77 | + constexpr T C3 = -0.00434102; |
| 78 | + |
| 79 | + // Map x to the range [-pi, pi] |
| 80 | + while (x < -2*M_PI) { |
| 81 | + x += 2.0 * M_PI; |
| 82 | + } |
| 83 | + while (x > 2*M_PI) { |
| 84 | + x -= 2.0 * M_PI; |
| 85 | + } |
| 86 | + |
| 87 | + // Calculate cos(x) using Chebyshev polynomial approximation |
| 88 | + T x2 = x * x; |
| 89 | + T result = C0 + x2 * (C1 + x2 * (C2 + x2 * C3)); |
| 90 | + return result; |
| 91 | +} |
| 92 | + |
| 93 | +template <typename T> |
| 94 | +inline half_float::half PromoteToHalf(T t, float maxColors) { |
| 95 | + half_float::half result((float)t / maxColors); |
| 96 | + return result; |
| 97 | +} |
| 98 | + |
| 99 | +template <typename D, typename T> |
| 100 | +inline D PromoteTo(T t, float maxColors) { |
| 101 | + D result = static_cast<D>((float)t / maxColors); |
| 102 | + return result; |
| 103 | +} |
| 104 | + |
| 105 | +template <typename T> |
| 106 | +inline T DemoteHalfTo(half t, float maxColors) { |
| 107 | + return (T)clamp(((float)t * (float)maxColors), 0.0f, (float)maxColors); |
| 108 | +} |
| 109 | + |
| 110 | +template <typename D, typename T> |
| 111 | +inline D DemoteTo(T t, float maxColors) { |
| 112 | + return (D)clamp(((float)t * (float)maxColors), 0.0f, (float)maxColors); |
| 113 | +} |
| 114 | + |
| 115 | +template <typename T> |
| 116 | +inline T CubicBSpline(T t) { |
| 117 | + T absX = abs(t); |
| 118 | + if (absX <= 1.0) { |
| 119 | + T doubled = absX * absX; |
| 120 | + T tripled = doubled * absX; |
| 121 | + return T(2.0 / 3.0) - doubled + (T(0.5) * tripled); |
| 122 | + } else if (absX <= 2.0) { |
| 123 | + return ((T(2.0) - absX) * (T(2.0) - absX) * (T(2.0) - absX)) / T(6.0); |
| 124 | + } else { |
| 125 | + return T(0.0); |
| 126 | + } |
| 127 | +} |
| 128 | + |
| 129 | +template <typename T> |
| 130 | +T SimpleCubic(T t, T A, T B, T C, T D) { |
| 131 | + T duplet = t * t; |
| 132 | + T triplet = duplet * t; |
| 133 | + T a = -A / T(2.0) + (T(3.0) * B) / T(2.0) - (T(3.0) * C) / T(2.0) + D / T(2.0); |
| 134 | + T b = A - (T(5.0) * B) / T(2.0) + T(2.0) * C - D / T(2.0); |
| 135 | + T c = -A / T(2.0) + C / T(2.0); |
| 136 | + T d = B; |
| 137 | + return a * triplet * T(3.0) + b * duplet + c * t + d; |
| 138 | +} |
| 139 | + |
| 140 | +template <typename T> |
| 141 | +inline T CubicHermite(T d, T p0, T p1, T p2, T p3) { |
| 142 | + constexpr T C = T(0.0); |
| 143 | + constexpr T B = T(0.0); |
| 144 | + T duplet = d * d; |
| 145 | + T triplet = duplet * d; |
| 146 | + T firstRow = ((T(-1/6.0)*B - C)*p0 + (T(-1.5)*B - C + T(2.0))*p1 + (T(1.5)*B + C - T(2.0))*p2 + (T(1/6.0)*B + C)*p3)*triplet; |
| 147 | + T secondRow = ((T(0.5)*B + 2*C)*p0 + (T(2.0)*B + C - T(3.0))*p1 + (T(-2.5)*B-T(2.0)*C + T(3.0))*p2 - C*p3)*duplet; |
| 148 | + T thirdRow = ((T(-0.5)*B - C)*p0 + (T(0.5)*B+C)*p2)*d; |
| 149 | + T fourthRow = (T(1.0/6.0)*B)*p0 + (T(-1.0/3.0)*B+T(1))*p1 + (T(1.0/6.0)*B)*p2; |
| 150 | + return firstRow + secondRow + thirdRow + fourthRow; |
| 151 | +} |
| 152 | + |
| 153 | +template <typename T> |
| 154 | +inline T CubicBSpline(T d, T p0, T p1, T p2, T p3) { |
| 155 | + // T t2 = t * t; |
| 156 | + // T a0 = -T(0.5) * p0 + T(1.5) * p1 - T(1.5) * p2 + T(0.5) * p3; |
| 157 | + // T a1 = p0 - T(2.5) * p1 + T(2.0f) * p2 - T(0.5) * p3; |
| 158 | + // T a2 = -T(0.5) * p0 + T(0.5) * p2; |
| 159 | + // T a3 = p1; |
| 160 | + // return (a0 * t * t2 + a1 * t2 + a2 * t + a3); |
| 161 | + constexpr T C = T(0.0); |
| 162 | + constexpr T B = T(1.0); |
| 163 | + T duplet = d * d; |
| 164 | + T triplet = duplet * d; |
| 165 | + T firstRow = ((T(-1/6.0)*B - C)*p0 + (T(-1.5)*B - C + T(2.0))*p1 + (T(1.5)*B + C - T(2.0))*p2 + (T(1/6.0)*B + C)*p3)*triplet; |
| 166 | + T secondRow = ((T(0.5)*B + 2*C)*p0 + (T(2.0)*B + C - T(3.0))*p1 + (T(-2.5)*B-T(2.0)*C + T(3.0))*p2 - C*p3)*duplet; |
| 167 | + T thirdRow = ((T(-0.5)*B - C)*p0 + (T(0.5)*B+C)*p2)*d; |
| 168 | + T fourthRow = (T(1.0/6.0)*B)*p0 + (T(-1.0/3.0)*B+T(1))*p1 + (T(1.0/6.0)*B)*p2; |
| 169 | + return firstRow + secondRow + thirdRow + fourthRow; |
| 170 | +} |
| 171 | + |
| 172 | +template <typename T> |
| 173 | +inline T FilterMitchell(T t) { |
| 174 | + T x = abs(t); |
| 175 | + |
| 176 | + if (x < 1.0f) |
| 177 | + return (T(16) + x*x*(T(21) * x - T(36)))/T(18); |
| 178 | + else if (x < 2.0f) |
| 179 | + return (T(32) + x*(T(-60) + x*(T(36) - T(7)*x)))/T(18); |
| 180 | + |
| 181 | + return T(0.0f); |
| 182 | +} |
| 183 | + |
| 184 | +template <typename T> |
| 185 | +T MitchellNetravali(T d, T p0, T p1, T p2, T p3) { |
| 186 | + constexpr T C = T(1.0/3.0); |
| 187 | + constexpr T B = T(1.0/3.0); |
| 188 | + T duplet = d * d; |
| 189 | + T triplet = duplet * d; |
| 190 | + T firstRow = ((T(-1/6.0)*B - C)*p0 + (T(-1.5)*B - C + T(2.0))*p1 + (T(1.5)*B + C - T(2.0))*p2 + (T(1/6.0)*B + C)*p3)*triplet; |
| 191 | + T secondRow = ((T(0.5)*B + 2*C)*p0 + (T(2.0)*B + C - T(3.0))*p1 + (T(-2.5)*B-T(2.0)*C + T(3.0))*p2 - C*p3)*duplet; |
| 192 | + T thirdRow = ((T(-0.5)*B - C)*p0 + (T(0.5)*B+C)*p2)*d; |
| 193 | + T fourthRow = (T(1.0/6.0)*B)*p0 + (T(-1.0/3.0)*B+T(1))*p1 + (T(1.0/6.0)*B)*p2; |
| 194 | + return firstRow + secondRow + thirdRow + fourthRow; |
| 195 | +} |
| 196 | + |
| 197 | +template <typename T> |
| 198 | +inline T sinc(T x) { |
| 199 | + if (x == 0.0) { |
| 200 | + return T(1.0); |
| 201 | + } else { |
| 202 | + return fastSin2(x) / x; |
| 203 | + } |
| 204 | +} |
| 205 | + |
| 206 | +template <typename T> |
| 207 | +inline T LanczosWindow(T x, const T a) { |
| 208 | + if (abs(x) < a) { |
| 209 | + return sinc(T(M_PI) * x) * sinc(T(M_PI) * x / a); |
| 210 | + } |
| 211 | + return T(0.0); |
| 212 | +} |
| 213 | + |
| 214 | +template <typename T> |
| 215 | +inline T HannWindow(T x, const T length) { |
| 216 | + if (abs(x) <= length / 2) { |
| 217 | + T cx = fastCos(T(M_PI)*x / length); |
| 218 | + return (1/length)*cx*cx; |
| 219 | + } else { |
| 220 | + return T(0); |
| 221 | + } |
| 222 | +} |
| 223 | + |
| 224 | +template <typename T> |
| 225 | +inline T CatmullRom(T x) { |
| 226 | + x = (T)abs(x); |
| 227 | + |
| 228 | + if (x < 1.0f) |
| 229 | + return T(1) - x*x*(T(2.5f) - T(1.5f)*x); |
| 230 | + else if (x < 2.0f) |
| 231 | + return T(2) - x*(T(4) + x*(T(0.5f)*x - T(2.5f))); |
| 232 | + |
| 233 | + return T(0.0f); |
| 234 | +} |
| 235 | + |
| 236 | +template <typename T> |
| 237 | +inline T CatmullRom(T x, T p0, T p1, T p2, T p3) { |
| 238 | + x = abs(x); |
| 239 | + |
| 240 | + if (x < T(1.0)) { |
| 241 | + T doublePower = x * x; |
| 242 | + T triplePower = doublePower * x; |
| 243 | + return T(0.5) * (((T(2.0) * p1) + |
| 244 | + (-p0 + p2) * x + |
| 245 | + (T(2.0) * p0 - T(5.0) * p1 + T(4.0) * p2 - p3) * doublePower + |
| 246 | + (-p0 + T(3.0) * p1 - T(3.0) * p2 + p3) * triplePower)); |
| 247 | + } |
| 248 | + return T(0.0); |
| 249 | +} |
| 250 | + |
| 251 | +template float fastSin1(float x); |
| 252 | +template float fastSin2(float x); |
| 253 | +template float fastCos(float x); |
| 254 | +template float PromoteTo<float, uint16_t>(uint16_t, float); |
| 255 | +template float PromoteTo<float, uint8_t>(uint8_t, float); |
| 256 | +template uint8_t DemoteTo<uint8_t, float>(float, float); |
| 257 | +template uint16_t DemoteTo<uint16_t, float>(float, float); |
| 258 | +template float CatmullRom(float x, float p0, float p1, float p2, float p3); |
| 259 | +template float LanczosWindow(float, const float); |
| 260 | +template float HannWindow(float, const float); |
| 261 | +template float SimpleCubic(float, float, float, float, float); |
| 262 | +template float CubicBSpline(float, float, float, float, float); |
| 263 | +template float CubicHermite(float, float, float, float, float); |
| 264 | +template float MitchellNetravali(float, float, float, float, float); |
0 commit comments