diff --git a/src/lanczos/lanczos_opencv.jl b/src/lanczos/lanczos_opencv.jl index 62f350df..b9b0a2c1 100644 --- a/src/lanczos/lanczos_opencv.jl +++ b/src/lanczos/lanczos_opencv.jl @@ -42,13 +42,26 @@ 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 <: AbstractLanczos end @@ -56,14 +69,9 @@ degree(::Lanczos4OpenCV) = 4 value_weights(::Lanczos4OpenCV, δx) = _lanczos4_opencv(δx) -# 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] - +# 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) function _lanczos4_opencv(δx) p_4 = π / 4 @@ -84,3 +92,42 @@ function _lanczos4_opencv(δx) 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_faithful(::Type{F}, l4_2d_cs, δx) where {F} + p_4 = π * F(0.25) + y0 = -(δx + F(3)) * p_4 + s0, c0 = sincos(y0) + cs = ntuple(8) do i + 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 + end + sum_cs = sum(cs) + normed_cs = ntuple(i -> cs[i] / sum_cs, Val(8)) + return normed_cs +end diff --git a/test/lanczos/runtests.jl b/test/lanczos/runtests.jl index ffd3f614..f683c324 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,17 @@ end # symmetry check interpolant = itp.(X) @test interpolant ≈ reverse(interpolant) +end +@testset "Lanczos OpenCV4 Faithful" begin + 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 end