filters: update ewa_lanczos4sharpest blur value
With little bit of bikeshedding, based on haasn code I made two arbitrary precision solvers, in Python using mpmath and in C using MPFR. They all agree on this value, even standard glibc double solver.
For reference code below. MPFR version is not worth sharing as it doesn't bring anything to the table and it is pain to read it anyway.
import mpmath as mp
mp.mp.dps = 1000
print(mp.mp)
R1 = mp.besseljzero(1, 1) / mp.pi
R2 = mp.besseljzero(1, 2) / mp.pi
R3 = mp.besseljzero(1, 3) / mp.pi
R4 = mp.besseljzero(1, 4) / mp.pi
R = R4
def jinc(x):
if x < mp.mpf(0):
return mp.mpf(1)
x *= mp.pi
return mp.mpf(2) * mp.besselj(1, x) / x
def lanczos(x):
if x >= R:
return mp.mpf(0)
return jinc(x) * jinc(x * R1/R)
def residual(scale):
return mp.nsum(lambda x, y: lanczos(scale * mp.sqrt(x**2 + y**2)), [1, R+1], [0, R+1])
print(f'R1: {mp.nstr(R1, 50)}')
print(f'R2: {mp.nstr(R2, 50)}')
print(f'R3: {mp.nstr(R3, 50)}')
print(f'R4: {mp.nstr(R4, 50)}')
print(f'R: {mp.nstr(R, 50)}')
root = mp.findroot(residual, 1)
blur = mp.mpf(1) / root
print(f'root: {mp.nstr(root, 50)}')
print(f'radius: {mp.nstr(R, 50)}')
print(f'blur: {mp.nstr(blur, 50)}')
print(f'new radius: {mp.nstr(R * blur, 50)}')Mathematica for the memes
jinc[x_] := BesselJ[1, x] / x;
j[x_] := If[x < 0, 1, 2 * jinc[x * Pi]];
lanczos[x_, r1_, r4_] := If[x >= r4, 0, j[x] * j[x * r1 / r4]];
residual[scale_, r1_, r_] := Sum[lanczos[scale * Sqrt[x^2 + y^2], r1, r], {x, 1,Floor[r], 1}, {y, 0, Floor[r], 1}];
R[r_]= BesselJZero[1, r] / Pi;
root = First[p /. Solve[residual[p, R[1], R[4]] == 0 && p > 1 && p < 2]]
N[1 / root, 50]Edited by Kacper Michajłow
