Skip to content

filters: update ewa_lanczos4sharpest blur value

Kacper Michajłow requested to merge kasper93/libplacebo:m into master

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

Merge request reports