From a127d7d39129ee648c7676d3c736b0700aa92a12 Mon Sep 17 00:00:00 2001 From: Ben Arthur Date: Fri, 24 Oct 2025 17:13:47 -0400 Subject: [PATCH 1/6] reduced precision for lanczos_opencv --- src/lanczos/lanczos_opencv.jl | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/src/lanczos/lanczos_opencv.jl b/src/lanczos/lanczos_opencv.jl index 62f350df..556847a2 100644 --- a/src/lanczos/lanczos_opencv.jl +++ b/src/lanczos/lanczos_opencv.jl @@ -65,20 +65,21 @@ const s45 = 0.70710678118654752440084436210485 const l4_2d_cs = SA[1 0; -s45 -s45; 0 1; s45 -s45; -1 0; s45 s45; 0 -1; -s45 s45] -function _lanczos4_opencv(δx) - p_4 = π / 4 +function _lanczos4_opencv(δx::T) where T + F = float(T) + p_4 = π / F(4) y0 = -(δx + 3) * p_4 - s0, c0 = sincos(y0) + s0::F, c0::F = sincos(y0) cs = ntuple(8) do i y = (δx + 4 - i) * p_4 # Improve precision of Lanczos OpenCV4 #451, avoid NaN if iszero(y) - y = eps(oneunit(y))/8 + y = eps(oneunit(y))/F(8) end # Numerator is the sin subtraction identity # It is equivalent to the following # f(δx,i) = sin( π/4*( 5*(i-1)-δx-3 ) ) - (l4_2d_cs[i, 1] * s0 + l4_2d_cs[i, 2] * c0) / y^2 + (F(l4_2d_cs[i, 1]) * s0 + F(l4_2d_cs[i, 2]) * c0) / y^2 end sum_cs = sum(cs) normed_cs = ntuple(i -> cs[i] / sum_cs, Val(8)) From 5f69ead662058083dbb42c635f342e5ba2721275 Mon Sep 17 00:00:00 2001 From: Ben Arthur Date: Fri, 24 Oct 2025 17:41:06 -0400 Subject: [PATCH 2/6] function barrier for float(F) --- src/lanczos/lanczos_opencv.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/lanczos/lanczos_opencv.jl b/src/lanczos/lanczos_opencv.jl index 556847a2..fd04c34d 100644 --- a/src/lanczos/lanczos_opencv.jl +++ b/src/lanczos/lanczos_opencv.jl @@ -54,7 +54,8 @@ struct Lanczos4OpenCV <: AbstractLanczos end degree(::Lanczos4OpenCV) = 4 -value_weights(::Lanczos4OpenCV, δx) = _lanczos4_opencv(δx) +value_weights(::Lanczos4OpenCV, δx::T) where {T} = + _lanczos4_opencv(float(T), float(T).(l4_2d_cs), δx) # s45 = sqrt(2)/2 const s45 = 0.70710678118654752440084436210485 @@ -65,8 +66,7 @@ const s45 = 0.70710678118654752440084436210485 const l4_2d_cs = SA[1 0; -s45 -s45; 0 1; s45 -s45; -1 0; s45 s45; 0 -1; -s45 s45] -function _lanczos4_opencv(δx::T) where T - F = float(T) +function _lanczos4_opencv(::Type{F}, l4_2d_cs, δx) where {F} p_4 = π / F(4) y0 = -(δx + 3) * p_4 s0::F, c0::F = sincos(y0) From 6dfc81456e441899bec414711a8e5967df0fce2f Mon Sep 17 00:00:00 2001 From: Ben Arthur Date: Mon, 27 Oct 2025 18:20:26 -0400 Subject: [PATCH 3/6] specify lanczos precision in type parameter --- src/lanczos/lanczos_opencv.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/lanczos/lanczos_opencv.jl b/src/lanczos/lanczos_opencv.jl index fd04c34d..7433b08d 100644 --- a/src/lanczos/lanczos_opencv.jl +++ b/src/lanczos/lanczos_opencv.jl @@ -50,11 +50,12 @@ export Lanczos4OpenCV Alternative implementation of Lanczos resampling using algorithm `lanczos4` function of OpenCV: https://github.com/opencv/opencv/blob/de15636724967faf62c2d1bce26f4335e4b359e5/modules/imgproc/src/resize.cpp#L917-L946 """ -struct Lanczos4OpenCV <: AbstractLanczos end +struct Lanczos4OpenCV{T} <: AbstractLanczos end +Lanczos4OpenCV() = Lanczos4OpenCV{Float64}() degree(::Lanczos4OpenCV) = 4 -value_weights(::Lanczos4OpenCV, δx::T) where {T} = +value_weights(::Lanczos4OpenCV{T}, δx) where {T} = _lanczos4_opencv(float(T), float(T).(l4_2d_cs), δx) # s45 = sqrt(2)/2 From 740fa93b6baa857fb7007a931cd0c57f11095fdb Mon Sep 17 00:00:00 2001 From: Ben Arthur Date: Mon, 27 Oct 2025 18:21:21 -0400 Subject: [PATCH 4/6] more faithfully follow opencv's lanczos code --- src/lanczos/lanczos_opencv.jl | 20 +++++++++----------- 1 file changed, 9 insertions(+), 11 deletions(-) diff --git a/src/lanczos/lanczos_opencv.jl b/src/lanczos/lanczos_opencv.jl index 7433b08d..e1a3d11d 100644 --- a/src/lanczos/lanczos_opencv.jl +++ b/src/lanczos/lanczos_opencv.jl @@ -68,19 +68,17 @@ const l4_2d_cs = SA[1 0; -s45 -s45; 0 1; s45 -s45; -1 0; s45 s45; 0 -1; -s45 s45 function _lanczos4_opencv(::Type{F}, l4_2d_cs, δx) where {F} - p_4 = π / F(4) - y0 = -(δx + 3) * p_4 - s0::F, c0::F = sincos(y0) + p_4 = π * F(0.25) + y0 = -(δx + F(3)) * p_4 + s0, c0 = sincos(y0) cs = ntuple(8) do i - y = (δx + 4 - i) * p_4 - # Improve precision of Lanczos OpenCV4 #451, avoid NaN - if iszero(y) - y = eps(oneunit(y))/F(8) + y = δx + 4 - i + if abs(y) >= F(1e-6) + y *= p_4 + (F(l4_2d_cs[i, 1]) * s0 + F(l4_2d_cs[i, 2]) * c0) / y^2 + else + F(1e30) end - # Numerator is the sin subtraction identity - # It is equivalent to the following - # f(δx,i) = sin( π/4*( 5*(i-1)-δx-3 ) ) - (F(l4_2d_cs[i, 1]) * s0 + F(l4_2d_cs[i, 2]) * c0) / y^2 end sum_cs = sum(cs) normed_cs = ntuple(i -> cs[i] / sum_cs, Val(8)) From 04dc3a8fd86540f0286399264078d82e850f2aaf Mon Sep 17 00:00:00 2001 From: Ben Arthur Date: Fri, 8 May 2026 16:14:29 -0400 Subject: [PATCH 5/6] preserve former lanczos opencv code --- src/lanczos/lanczos_opencv.jl | 73 ++++++++++++++++++++++++++++------- test/lanczos/runtests.jl | 9 ++++- 2 files changed, 67 insertions(+), 15 deletions(-) diff --git a/src/lanczos/lanczos_opencv.jl b/src/lanczos/lanczos_opencv.jl index e1a3d11d..b9b0a2c1 100644 --- a/src/lanczos/lanczos_opencv.jl +++ b/src/lanczos/lanczos_opencv.jl @@ -42,32 +42,79 @@ the use of this software, even if advised of the possibility of such damage. =# using StaticArrays -export Lanczos4OpenCV +export Lanczos4OpenCV, Lanczos4OpenCVFaithful + +# s45 = sqrt(2)/2 +const s45 = 0.70710678118654752440084436210485 + +# l4_2d_cs is a lookup table that could be generated by +# x = (0:7)*45*5 +# l4_2d_cs = [cosd.(x) sind.(x)]' +const l4_2d_cs = SA[1 0; -s45 -s45; 0 1; s45 -s45; -1 0; s45 s45; 0 -1; -s45 s45] + """ Lanczos4OpenCV() -Alternative implementation of Lanczos resampling using algorithm `lanczos4` function of OpenCV: +Alternative implementation of Lanczos resampling using algorithm `lanczos4` +function of OpenCV: + https://github.com/opencv/opencv/blob/de15636724967faf62c2d1bce26f4335e4b359e5/modules/imgproc/src/resize.cpp#L917-L946 + +If your data are Float32, consider using `Lanczos4OpenCVFaithful` instead. """ -struct Lanczos4OpenCV{T} <: AbstractLanczos end -Lanczos4OpenCV() = Lanczos4OpenCV{Float64}() +struct Lanczos4OpenCV <: AbstractLanczos end degree(::Lanczos4OpenCV) = 4 -value_weights(::Lanczos4OpenCV{T}, δx) where {T} = - _lanczos4_opencv(float(T), float(T).(l4_2d_cs), δx) +value_weights(::Lanczos4OpenCV, δx) = _lanczos4_opencv(δx) -# s45 = sqrt(2)/2 -const s45 = 0.70710678118654752440084436210485 +# a simple refactoring of `_lanczos4_opencv` to use Float32 throughout results +# in s0 == c0 and so the numerator of cs[4] is 0 as is normed_cs[4] (when it +# should be 1) -# l4_2d_cs is a lookup table that could be generated by -# x = (0:7)*45*5 -# l4_2d_cs = [cosd.(x) sind.(x)]' -const l4_2d_cs = SA[1 0; -s45 -s45; 0 1; s45 -s45; -1 0; s45 s45; 0 -1; -s45 s45] +function _lanczos4_opencv(δx) + p_4 = π / 4 + y0 = -(δx + 3) * p_4 + s0, c0 = sincos(y0) + cs = ntuple(8) do i + y = (δx + 4 - i) * p_4 + # Improve precision of Lanczos OpenCV4 #451, avoid NaN + if iszero(y) + y = eps(oneunit(y))/8 + end + # Numerator is the sin subtraction identity + # It is equivalent to the following + # f(δx,i) = sin( π/4*( 5*(i-1)-δx-3 ) ) + (l4_2d_cs[i, 1] * s0 + l4_2d_cs[i, 2] * c0) / y^2 + end + sum_cs = sum(cs) + normed_cs = ntuple(i -> cs[i] / sum_cs, Val(8)) + return normed_cs +end + + +""" + Lanczos4OpenCVFaithful{T}() + +Similar to Lanczos4OpenCV(), but mirrors the calculations more closely to +reduce numerical discrepancies with the OpenCV implementation, and supports +reduced precision (e.g. Float32, the default). +""" +struct Lanczos4OpenCVFaithful{T} <: AbstractLanczos end +Lanczos4OpenCVFaithful() = Lanczos4OpenCVFaithful{Float32}() + +degree(::Lanczos4OpenCVFaithful) = 4 + +value_weights(::Lanczos4OpenCVFaithful{T}, δx) where {T} = + _lanczos4_opencv_faithful(float(T), float(T).(l4_2d_cs), δx) +# main differences with `_lanczos4_opencv` are (1) the criterion for preventing a +# division by zero is `< pi/4 * 1e-6` (instead of `iszero`) and (2) the resulting +# coefficient is replaced with `1e30` (instead of using `(eps()/8)^2` in the +# denominator) -function _lanczos4_opencv(::Type{F}, l4_2d_cs, δx) where {F} +function _lanczos4_opencv_faithful(::Type{F}, l4_2d_cs, δx) where {F} p_4 = π * F(0.25) y0 = -(δx + F(3)) * p_4 s0, c0 = sincos(y0) diff --git a/test/lanczos/runtests.jl b/test/lanczos/runtests.jl index ffd3f614..6e495da1 100644 --- a/test/lanczos/runtests.jl +++ b/test/lanczos/runtests.jl @@ -28,11 +28,9 @@ using Interpolations: degree, @test 1 < itp(1.5) < 2 @test 99 < itp(99.5) < 100 - # symmetry check interpolant = itp.(X) @test interpolant ≈ reverse(interpolant) - end @testset "Lanczos OpenCV4" begin @@ -60,7 +58,14 @@ end # symmetry check interpolant = itp.(X) @test interpolant ≈ reverse(interpolant) +end +@testset "Lanczos OpenCV4 Faithful" begin + @test 1f0 == Interpolations._lanczos4_opencv_faithful(Float32, Interpolations.l4_2d_cs, 0.0f0)[4] + @test 1.0 == Interpolations._lanczos4_opencv_faithful(Float64, Interpolations.l4_2d_cs, 0.0)[4] + @test 1.0 == Interpolations._lanczos4_opencv(0.0)[4] + # were _lanczos4_opencv to be refactored to use Float32 throughout, + # the 4th element returned would be 0.0 end end From 9ba58a097323cd90054d595d58a22a21d0f12a62 Mon Sep 17 00:00:00 2001 From: Ben Arthur Date: Fri, 8 May 2026 16:40:35 -0400 Subject: [PATCH 6/6] increase test coverage --- test/lanczos/runtests.jl | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/test/lanczos/runtests.jl b/test/lanczos/runtests.jl index 6e495da1..f683c324 100644 --- a/test/lanczos/runtests.jl +++ b/test/lanczos/runtests.jl @@ -61,9 +61,12 @@ end end @testset "Lanczos OpenCV4 Faithful" begin - @test 1f0 == Interpolations._lanczos4_opencv_faithful(Float32, Interpolations.l4_2d_cs, 0.0f0)[4] - @test 1.0 == Interpolations._lanczos4_opencv_faithful(Float64, Interpolations.l4_2d_cs, 0.0)[4] - @test 1.0 == Interpolations._lanczos4_opencv(0.0)[4] + idx = Interpolations.degree(Lanczos4OpenCVFaithful{Float32}()) + @test 1f0 == Interpolations.value_weights(Lanczos4OpenCVFaithful{Float32}(), 0.0f0)[idx] + idx = Interpolations.degree(Lanczos4OpenCVFaithful{Float64}()) + @test 1.0 == Interpolations.value_weights(Lanczos4OpenCVFaithful{Float64}(), 0.0)[idx] + idx = Interpolations.degree(Lanczos4OpenCV()) + @test 1.0 == Interpolations.value_weights(Lanczos4OpenCV(), 0.0)[idx] # were _lanczos4_opencv to be refactored to use Float32 throughout, # the 4th element returned would be 0.0 end