Prompted by an mpv video player user
https://mpv.io/ (and software developer) who found the EWA LanczosSharpest promising and asked about the optimal values for higher number of lobes, I revised the computation of the "optimal" deblurs after finding suspicious stuff in the original versions. Not a bad thing, because I found mistakes: The values I gave for LanczosSharpest4 and LanczosSharpest5 were slightly off. My earlier code nailed the 2-lobe deblur value (the LanczosSharpest2 one), but produced values for the 3- and 4-lobe blurs which were off by 0.1% and 0.4%, respectively.
With the usual caveat that I re-did this rather fast, and that the numerical computations I perform have limited accuracy (when you and I use Bessel functions, we are generally relying on approximations which fall somewhat short of full double precision accuracy), here is the revisited set of blur values, from 2-lobes to 8-lobes:
Code: Select all
0.8882642150854034, 0.8845100233858514, 0.8880424053585475, 0.8845552452577251, 0.8876859533084833, 0.8851157217316564, 0.88670025197296
The corresponding disc radii are (when we deblur, the original disc shrinks):
Code: Select all
1.983609994601509, 3.7512626128371274, 4.655797087963903, 5.523093681748111, 6.4310715705557175, 7.298128670225191, 8.19833537345964
Finally, here is the rather user-unfriendly Axiom code that I used to nail the optimal values, and that can be used to check what I did:
Code: Select all
-- Author: Nicolas Robidoux
-- Computation of the deblur of the EWA Lanczos (Jinc-windowed
-- Jinckernel that makes the convolution, under no-op, the closest to
-- the identity in infinity linear operator norm. Equivalently, deblur
-- that makes the weight of the center pixel equal to 1. Equivalently,
-- deblur that makes the sum of the weights of the other pixels equal
-- to zero. One may argue that this gives a good approximation of the
-- sharpest resampling Jinc-windowed Jinc kernel.
-- Axiom Computer Algebra System code http://www.axiom-developer.org/
-- Run with )read pathToThisFile or by cut and pasting into the Axiom terminal window.
)cl a
digits(100)
-- Set the number of lobes that the code will check. You may change
-- this to any integer between 2 and 8.
numberOfLobes := 2
-- Jinc function roots from
-- http://cose.math.bas.bg/webMathematica/webComputing/BesselZeros.jsp
-- Note: The Bessel function J_1 has root at 0, which Jinc does not,
-- so the first root of Jinc is the second root of J_1 etc.
R1 := 3.8317059702075123156144358863081607665645452742878019287622989899188393095190114702141128747574231267244747330972398596077082935393429507463756833072973507947354388612251605091329310067675478669096098
R2 := 7.0155866698156187535370499814765247432763115029113138960553778269854960155020186630727149301794664578066071840984767118820251703769720326548673375642560309232023970190386543841800591275430452193373925
R3 := 13.323691936314223032393684126947876751216644731357865785477571526496567063347304782547119017948839951109445358057000268451481978642541401546368509755789722112502852743021277168603949713032474504688907
R4 := 16.470630050877632812552460470989551449438126822273125769944420007944442567640743899002931497029400154465903980274023428343447658407656485921770613343351595779252911253088717992488254261428837132098213
R5 := 19.615858510468242021125065884137509850247402661880544647351444764470801101727298271485389161157517698967796725446148833952670733337955890767096174193012715265947519553741418696730317913690178285283661
R6 := 22.760084380592771898053005152182257592905370738073226872005077130254273562589664194570326702300201120545899255931397351763974854463339052768274662068684438595699384628576045067718787233238726822557551
R7 := 25.903672087618382625495855445979874287905427031367247641367104494345056714163165058881472268367878871249811307040417420294816882799155844275939749558981998118306731999865474455453558204118663351897861
R8 := 29.046828534916855066647819883531961100414171793083875666040125573376578395536196209367825338086768931611841103790761895902010700540781204477017596137035335307002504273161322353313122231007878096884274
Roots := [ R1, R2, R3, R4, R5, R6, R7, R8 ]
-- Select the root that matches the number of lobes.
-- The EWA Lanczos (Jinc-windowed Jinc) that matches R2 is the two lobe one etc.
Root := Roots.numberOfLobes
Radius := Root / %pi
jinc x == ((besselJ(1,%pi*x))::Float)/x
R1overRadius := R1 / Radius
wind x == (besselJ(1,x*R1overRadius)::Float)/x
l x == if (x<Radius) then ( wind(x) * jinc(x) ) else 0.
z := l 1.e-128
-- The code assumes that we are deblurring (meaning that we are making the effective radius smaller, not larger).
-- As a result, we know that if we are farther away than Radius, the weight is 0.
integerRadius := (floor(Radius)) :: INT
quadrant r == reduce( +, [ reduce( +, [ l(r*sqrt((i^2+j^2)::Float)) for i in 1..integerRadius ] ) for j in 0..integerRadius ] )
residual r == quadrant(r) / ( z + 4. * quadrant(r) )
-----
-- Number of lobes = 2:
[ residual(1.12579115877571804+i*.000000000000000001) :: SF for i in 0..10 ]
-- Check: Residual at the estimated minimizer
minimizer2 := 1.125791158775718044
residual2 := residual(minimizer2)
-- Actual blur value:
blur2 := 1. / minimizer2
-----
-- Number of lobes = 3:
[ residual(1.13056943794945355+i*.000000000000000001) :: SF for i in 0..10 ]
minimizer3 := 1.130569437949453554
residual3 := residual(minimizer3)
blur3 := 1 / minimizer3
-----
-- Number of lobes = 4:
[ residual(1.12607235191234992+i*.000000000000000001) :: SF for i in 0..10 ]
minimizer4 := 1.126072351912349924
residual4 := residual(minimizer4)
blur4 := 1 / minimizer4
-----
-- Number of lobes = 5:
[ residual(1.13051163888428338+i*.000000000000000001) :: SF for i in 0..10 ]
minimizer5 := 1.13051163888428339
residual5 := residual(minimizer5)
blur5 := 1 / minimizer5
-----
-- Number of lobes = 6:
[ residual(1.12652452849221331+i*.000000000000000001) :: SF for i in 0..10 ]
minimizer6 := 1.126524528492213318
residual6 := residual(minimizer6)
blur6 := 1 / minimizer6
-----
-- Number of lobes = 7:
[ residual(1.12979577183826532+i*.000000000000000001) :: SF for i in 0..10 ]
minimizer7 := 1.129795771838265326
residual7 := residual(minimizer7)
blur7 := 1 / minimizer7
-----
-- Number of lobes = 8:
[ residual(1.12777683075531042+i*.000000000000000001) :: SF for i in 0..10 ]
minimizer8 := 1.127776830755310429
residual8 := residual(minimizer8)
blur8 := 1 / minimizer8
-----
[ residual2, residual3, residual4, residual5, residual6, residual7, residual8 ]
%.(numberOfLobes-1)
blur := [ blur2, blur3, blur4, blur5, blur6, blur7, blur8 ]
radius := [ blur.i * Roots.(i+1) / %pi for i in 1..7 ]
[ i for i in 2..8 ]
blur :: List SF
radius :: List SF
-- Final output:
-- Numbers of lobes for which the deblurs were computed:
-- [ 2, 3, 4, 5, 6, 7, 8]
-- Deblurs in the same order:
-- [ 0.8882642150854034, 0.8845100233858514, 0.8880424053585475, 0.8845552452577251, 0.8876859533084833, 0.8851157217316564, 0.88670025197296 ]
-- Nominal radii of the support of the corresponding Jinc-windowed Jincs:
-- [ 1.983609994601509, 3.7512626128371274, 4.655797087963903, 5.523093681748111, 6.4310715705557175, 7.298128670225191, 8.19833537345964 ]