Merge branch 'master' into master

Author: mkitti

.gitignore

@@ -26,3 +26,4 @@
 ## Build folder of Documenter.jl
 docs/build/
 Manifest.toml
+.vscode

.travis.yml

@@ -19,3 +19,6 @@ jobs:
                                                Pkg.instantiate()'
         - julia --project=docs/ docs/make.jl
       after_success: skip
+  allow_failures:
+    - julia: nightly
+

Project.toml

@@ -1,6 +1,6 @@
 name = "Interpolations"
 uuid = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
-version = "0.12.10"
+version = "0.13.0"
 
 [deps]
 AxisAlgorithms = "13072b0f-2c55-5437-9ae7-d433b7a33950"

src/Interpolations.jl

@@ -473,5 +473,6 @@ include("utils.jl")
 include("io.jl")
 include("convenience-constructors.jl")
 include("deprecations.jl")
+include("lanczos/lanczos.jl")
 
 end # module

src/b-splines/cubic.jl

@@ -7,7 +7,7 @@ end
 """
     Cubic(bc::BoundaryCondition)
 
-Indicate that the corresponding axis should use quadratic interpolation.
+Indicate that the corresponding axis should use cubic interpolation.
 
 # Extended help
 

src/lanczos/lanczos.jl

@@ -0,0 +1,132 @@
+using Base.Cartesian
+using StaticArrays
+
+export Lanczos, Lanczos4OpenCV
+
+abstract type AbstractLanczos <: InterpolationType end
+
+"""
+    Lanczos{N}(a=4)
+
+Lanczos resampling via a kernel with scale parameter `a` and support over `N` neighbors.
+
+This form of interpolation is merely the discrete convolution of the samples with a Lanczos kernel of size `a`. The size is directly related to how "far" the interpolation will reach for information, and has `O(N^2)` impact on runtime. An alternative implementation matching `lanczos4` from OpenCV is available as Lanczos4OpenCV.
+"""
+struct Lanczos{N} <: AbstractLanczos
+    a::Int
+    
+    function Lanczos{N}(a) where N
+        N < a && @warn "Using a smaller support than scale for Lanczos window. Proceed with caution."
+        new{N}(a)
+    end
+end
+
+Lanczos(a=4) = Lanczos{a}(a)
+
+"""
+    LanczosInterpolation
+"""
+struct LanczosInterpolation{T,N,IT <: DimSpec{AbstractLanczos},A <: AbstractArray{T,N},P <: Tuple{Vararg{AbstractArray,N}}} <: AbstractInterpolation{T,N,IT}
+    coefs::A
+    parentaxes::P
+    it::IT
+end
+
+@generated degree(::Lanczos{N}) where {N} = :($N)
+
+getknots(itp::LanczosInterpolation) = axes(itp)
+coefficients(itp::LanczosInterpolation) = itp.coefs
+itpflag(itp::LanczosInterpolation) = itp.it
+
+size(itp::LanczosInterpolation) = map(length, itp.parentaxes)
+axes(itp::LanczosInterpolation) = itp.parentaxes
+lbounds(itp::LanczosInterpolation) = map(first, itp.parentaxes)
+ubounds(itp::LanczosInterpolation) = map(last, itp.parentaxes)
+
+function interpolate(A::AbstractArray{T}, it::AbstractLanczos) where T
+    Apad = copy_with_padding(float(T), A, it)
+    return LanczosInterpolation(Apad, axes(A), it)
+end
+
+@inline function (itp::LanczosInterpolation{T,N})(x::Vararg{<:Number,N}) where {T,N}
+    @boundscheck (checkbounds(Bool, itp, x...) || Base.throw_boundserror(itp, x))
+    wis = weightedindexes((value_weights,), itpinfo(itp)..., x)
+    itp.coefs[wis...]
+end
+
+function weightedindex_parts(fs, it::AbstractLanczos, ax::AbstractUnitRange{<:Integer}, x)
+    pos, δx = positions(it, ax, x)
+    (position = pos, coefs = fmap(fs, it, δx))
+end
+
+function positions(it::AbstractLanczos, ax, x)
+    xf = floorbounds(x, ax)
+    δx = x - xf
+    fast_trunc(Int, xf) - degree(it) + 1, δx
+end
+
+function value_weights(it::Lanczos, δx::S) where S
+    N = degree(it)
+    # short-circuit if integral
+    isinteger(δx) && return ntuple(i->convert(float(S), i == N - δx), Val(2N))
+
+    # LUTs
+    #it.a === N === 4 && return _lanczos4(δx)
+
+    cs = ntuple(i -> lanczos(N - i + δx, it.a, N), Val(2N))
+    sum_cs = sum(cs)
+    normed_cs = ntuple(i -> cs[i] / sum_cs, Val(length(cs)))
+    return normed_cs
+end
+
+function padded_axis(ax::AbstractUnitRange, it::AbstractLanczos)
+    N = degree(it)
+    return first(ax) - N + 1:last(ax) + N
+end
+
+# precise implementations for fast evaluation of common kernels
+
+"""
+    lanczos(x, a, n=a)
+
+Implementation of the [Lanczos kernel](https://en.wikipedia.org/wiki/Lanczos_resampling)
+"""
+lanczos(x::T, a::Integer, n=a) where {T} = abs(x) < n ? T(sinc(x) * sinc(x / a)) : zero(T)
+
+"""
+    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
+
+degree(::Lanczos4OpenCV) = 4
+
+value_weights(::Lanczos4OpenCV, δx::S) where S = ifelse(isinteger(δx),ntuple(i->convert(float(S), i == 4 - δx), Val(8)) ,_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]
+
+
+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
+        # 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

test/lanczos/runtests.jl

@@ -0,0 +1,64 @@
+using Interpolations: degree,
+                      getknots,
+                      coefficients,
+                      itpflag,
+                      DimSpec,
+                      AbstractLanczos,
+                      lbounds,
+                      ubounds
+
+@testset "Lanczos" begin
+
+@testset "Lanczos($N)" for N in 2:4
+    X = 1:100
+    X = [X; reverse(X)[2:end]]
+    itp = interpolate(X, Lanczos(N))
+
+    # properties
+    @test getknots(itp) == axes(itp) == axes(X)
+    @test coefficients(itp)[axes(itp)...] == X
+    @test itpflag(itp) isa Lanczos{N}
+    @test size(itp) == size(X)
+    @test lbounds(itp) == (firstindex(X),)
+    @test ubounds(itp) == (lastindex(X),)
+
+    # recovers input points exactly
+    @test itp.(X) == X
+
+    @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
+    X = 1:100
+    X = [X; reverse(X)[2:end]]
+    itp = interpolate(X, Lanczos4OpenCV())
+
+    # properties
+    @test getknots(itp) == axes(itp) == axes(X)
+    @test coefficients(itp)[axes(itp)...] == X
+    @test itpflag(itp) isa Lanczos4OpenCV
+    @test size(itp) == size(X)
+    @test lbounds(itp) == (firstindex(X),)
+    @test ubounds(itp) == (lastindex(X),)
+
+    # recovers input points exactly
+    @test itp.(X) == X
+
+    @test 1 < itp(1.5) < 2
+    @test 99 < itp(99.5) < 100
+
+
+    # symmetry check
+    interpolant = itp.(X)
+    @test interpolant ≈ reverse(interpolant)
+
+end
+
+end

test/runtests.jl

@@ -28,6 +28,9 @@ const isci = get(ENV, "CI", "") in ("true", "True")
     # monotonic tests
     include("monotonic/runtests.jl")
 
+    # lanczos tests
+    include("lanczos/runtests.jl")
+
     # test gradient evaluation
     include("gradient.jl")
     isci && println("finished gradient")