Fix/improve construction of NoInterp objects (#177)

This might seem strange, but it's a convenient way of doing *extrapolation* on arrays at integer coordinates.

Author: timholy

src/b-splines/indexing.jl

@@ -94,6 +94,10 @@ function getindex_return_type(::Type{BSplineInterpolation{T,N,TCoefs,IT,GT,Pad}}
     _promote_mul(eltype(TCoefs), I)
 end
 
+@noinline function getindex_return_type(::Type{Tel}, ::Type{BSplineInterpolation{T,N,TCoefs,IT,GT,Pad}}, argtypes::Tuple) where {Tel,T,N,TCoefs,IT<:DimSpec{BSpline},GT<:DimSpec{GridType},Pad}
+    reduce(_promote_mul, eltype(TCoefs), (Tel, argtypes...))
+end
+
 @generated function gradient!(g::AbstractVector, itp::BSplineInterpolation{T,N}, xs::Number...) where {T,N}
     length(xs) == N || error("Can only be called with $N indexes")
     gradient_impl(itp)

src/extrapolation/filled.jl

@@ -6,31 +6,32 @@ mutable struct FilledExtrapolation{T,N,ITP<:AbstractInterpolation,IT,GT,FT} <: A
 end
 
 function FilledExtrapolation(itp::AbstractInterpolation{T,N,IT,GT}, fillvalue) where {T,N,IT,GT}
-    FilledExtrapolation{T,N,typeof(itp),IT,GT,typeof(fillvalue)}(itp, fillvalue)
+    Te = promote_type(T,typeof(fillvalue))
+    FilledExtrapolation{Te,N,typeof(itp),IT,GT,typeof(fillvalue)}(itp, fillvalue)
 end
 
 Base.parent(A::FilledExtrapolation) = A.itp
 
 """
 `extrapolate(itp, fillvalue)` creates an extrapolation object that returns the `fillvalue` any time the indexes in `itp[x1,x2,...]` are out-of-bounds.
 """
-extrapolate(itp::AbstractInterpolation{T,N,IT,GT}, fillvalue) where {T,N,IT,GT} = FilledExtrapolation(itp, convert(eltype(itp), fillvalue))
+extrapolate(itp::AbstractInterpolation{T,N,IT,GT}, fillvalue) where {T,N,IT,GT} = FilledExtrapolation(itp, fillvalue)
 
 function getindex_impl(fitp::Type{FilledExtrapolation{T,N,ITP,IT,GT,FT}}, args) where {T,N,ITP,IT,GT,FT}
-   n = length(args)
-   n == N || return error("Must index $(N)-dimensional interpolation objects with $(nindexes(N))")
-
-   Tret = FT<:Number ? getindex_return_type(ITP, args) : FT
-   meta = Expr(:meta, :inline)
-   quote
-       $meta
-       # Check to see if we're in the extrapolation region, i.e.,
-       # out-of-bounds in an index
-       inds_etp = indices(fitp)
-       @nexprs $N d->((args[d] < lbound(fitp, d, inds_etp[d]) || args[d] > ubound(fitp, d, inds_etp[d]))) && return convert($Tret, fitp.fillvalue)::$Tret
-       # In the interpolation region
-       return convert($Tret, getindex(fitp.itp,args...))::$Tret
-   end
+    n = length(args)
+    n == N || return error("Must index $(N)-dimensional interpolation objects with $(nindexes(N))")
+
+    Tret = FT<:Number ? getindex_return_type(T, ITP, args) : FT
+    meta = Expr(:meta, :inline)
+    quote
+        $meta
+        # Check to see if we're in the extrapolation region, i.e.,
+        # out-of-bounds in an index
+        inds_etp = indices(fitp)
+        @nexprs $N d->((args[d] < lbound(fitp, d, inds_etp[d]) || args[d] > ubound(fitp, d, inds_etp[d]))) && return convert($Tret, fitp.fillvalue)::$Tret
+        # In the interpolation region
+        return convert($Tret, getindex(fitp.itp,args...))::$Tret
+    end
 end
 
 

src/nointerp/nointerp.jl

@@ -1,3 +1,9 @@
+function interpolate(A::AbstractArray, ::NoInterp, gt::GT) where {GT<:DimSpec{GridType}}
+    interpolate(Int, eltype(A), A, NoInterp(), gt)
+end
+
+iextract(::Type{NoInterp}, d) = NoInterp
+
 function define_indices_d(::Type{NoInterp}, d, pad)
     symix, symx = Symbol("ix_",d), Symbol("x_",d)
     :($symix = convert(Int, $symx))
@@ -27,4 +33,6 @@ function count_interp_dims(it::Type{IT}, N) where IT<:Tuple{Vararg{Interpolation
         n += count_interp_dims(p, 1)
     end
     n
-end
+end
+
+prefilter(::Type{TWeights}, ::Type{TC}, A, ::Type{IT},::Type{GT}) where {TWeights, TC, IT<:NoInterp, GT<:GridType} = A, Val{0}()

test/nointerp.jl

@@ -0,0 +1,19 @@
+module NoInterpTests
+println("Testing NoInterp...")
+using Interpolations, Base.Test
+
+a = reshape(1:12, 3, 4)
+ai = interpolate(a, NoInterp(), OnGrid())
+@test eltype(ai) == Int
+@test ai[1,1] == 1
+@test ai[3, 3] == 9
+@test_throws InexactError ai[2.2, 2]
+@test_throws InexactError ai[2, 2.2]
+
+ae = extrapolate(ai, NaN)
+@test eltype(ae) == Float64
+@test ae[1,1] === 1.0
+@test ae[0,1] === NaN
+@test_throws InexactError ae[1.5,2]
+
+end

test/runtests.jl

@@ -8,6 +8,7 @@ using Interpolations
 include("extrapolation/runtests.jl")
 # b-spline interpolation tests
 include("b-splines/runtests.jl")
+include("nointerp.jl")
 
 # scaling tests
 include("scaling/runtests.jl")