WIP: reduced precision for lanczos_opencv

src/lanczos/lanczos_opencv.jl

@@ -42,28 +42,36 @@ 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
 
 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

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