Merge pull request #9 from tlycken/linear-extrap

More boundary conditions for quadratic B-splines

src/Interpolations.jl

@@ -3,7 +3,7 @@ module Interpolations
 using Base.Cartesian
 using Compat
 
-import Base: size, eltype, getindex
+import Base: size, eltype, getindex, ndims
 
 export 
     Interpolation,
@@ -13,10 +13,19 @@ export
     ExtrapError,
     ExtrapNaN,
     ExtrapConstant,
+    ExtrapLinear,
+    ExtrapReflect,
+    ExtrapPeriodic,
     OnCell,
     OnGrid,
-    ExtendInner,
-    Flat
+    Flat,
+    Line,
+    Free,
+    Periodic,
+
+    degree,
+    boundarycondition,
+    gridrepresentation
 
 abstract Degree{N}
 
@@ -26,8 +35,11 @@ type OnCell <: GridRepresentation end
 
 abstract BoundaryCondition
 type None <: BoundaryCondition end
-type ExtendInner <: BoundaryCondition end
 type Flat <: BoundaryCondition end
+type Line <: BoundaryCondition end
+typealias Natural Line
+type Free <: BoundaryCondition end
+type Periodic <: BoundaryCondition end
 
 abstract InterpolationType{D<:Degree,BC<:BoundaryCondition,GR<:GridRepresentation}
 
@@ -49,27 +61,72 @@ end
 # However, all prefilters copy the array, so do that here as well
 prefilter{T,N,IT<:InterpolationType}(A::AbstractArray{T,N}, ::IT) = copy(A)
 
-size(itp::Interpolation, d::Integer) = size(itp.coefs, d)
-size(itp::Interpolation) = size(itp.coefs)
+size{T,N,IT<:InterpolationType}(itp::Interpolation{T,N,IT}, d::Integer) =
+    size(itp.coefs, d) - 2*padding(IT())
+size(itp::Interpolation) = tuple([size(itp,i) for i in 1:ndims(itp)]...)
+ndims(itp::Interpolation) = ndims(itp.coefs)
 eltype(itp::Interpolation) = eltype(itp.coefs)
 
 offsetsym(off, d) = off == -1 ? symbol(string("ixm_", d)) :
                     off ==  0 ? symbol(string("ix_", d)) :
                     off ==  1 ? symbol(string("ixp_", d)) :
                     off ==  2 ? symbol(string("ixpp_", d)) : error("offset $off not recognized")
 
+boundarycondition{D,BC<:BoundaryCondition}(::InterpolationType{D,BC}) = BC()
+boundarycondition{T,N,IT}(::Interpolation{T,N,IT}) = boundarycondition(IT())
 gridrepresentation{D,BC,GR<:GridRepresentation}(::InterpolationType{D,BC,GR}) = GR()
+gridrepresentation{T,N,IT}(::Interpolation{T,N,IT}) = gridrepresentation(IT())
 degree{D<:Degree,BC,GR}(::InterpolationType{D,BC,GR}) = D()
+degree{T,N,IT}(::Interpolation{T,N,IT}) = degree(IT())
 
 include("constant.jl")
 include("linear.jl")
 include("quadratic.jl")
 
+# If nothing else is specified, don't pad at all
+padding(::InterpolationType) = 0
+padding{T,N,IT<:InterpolationType}(::Interpolation{T,N,IT}) = padding(IT())
 
-# This creates getindex methods for all supported combinations
-for IT in (Constant{OnCell},Linear{OnGrid},Quadratic{ExtendInner,OnCell},Quadratic{Flat,OnCell})
-    for EB in (ExtrapError,ExtrapNaN,ExtrapConstant)
+function pad_size_and_index(sz::Tuple, pad)
+    sz = Int[s+2pad for s in sz]
+    N = length(sz)
+    ind = cell(N)
+    for i in 1:N
+        ind[i] = 1+pad:sz[i]-pad
+    end
+    sz, ind
+end
+function copy_with_padding(A, it::InterpolationType)
+    pad = padding(it)
+    sz,ind = pad_size_and_index(size(A), pad)
+    coefs = fill(convert(eltype(A), 0), sz...)
+    coefs[ind...] = A
+    coefs, pad
+end
 
+# This creates getindex methods for all supported combinations
+for IT in (
+        Constant{OnGrid},
+        Constant{OnCell},
+        Linear{OnGrid},
+        Linear{OnCell},
+        Quadratic{Flat,OnCell},
+        Quadratic{Flat,OnGrid},
+        Quadratic{Line,OnGrid},
+        Quadratic{Line,OnCell},
+        Quadratic{Free,OnGrid},
+        Quadratic{Free,OnCell},
+        Quadratic{Periodic,OnGrid},
+        Quadratic{Periodic,OnCell},
+    )
+    for EB in (
+            ExtrapError,
+            ExtrapNaN,
+            ExtrapConstant,
+            ExtrapLinear,
+            ExtrapReflect,
+            ExtrapPeriodic,
+        )
         eval(:(function getindex{T}(itp::Interpolation{T,1,$IT,$EB}, x::Real, d)
             d == 1 || throw(BoundsError())
             itp[x]
@@ -108,24 +165,40 @@ for IT in (Constant{OnCell},Linear{OnGrid},Quadratic{ExtendInner,OnCell},Quadrat
 end
 
 # This creates prefilter specializations for all interpolation types that need them
-for IT in (Quadratic{ExtendInner,OnCell},Quadratic{Flat,OnCell})
+for IT in (
+        Quadratic{Flat,OnCell},
+        Quadratic{Flat,OnGrid},
+        Quadratic{Line,OnGrid},
+        Quadratic{Line,OnCell},
+        Quadratic{Free,OnGrid},
+        Quadratic{Free,OnCell},
+        Quadratic{Periodic,OnGrid},
+        Quadratic{Periodic,OnCell},
+    )
     @ngenerate N promote_type_grid(T, x...) function prefilter{T,N}(A::Array{T,N},it::IT)
-        ret = similar(A)
-        szs = collect(size(A))
-        strds = collect(strides(A))
+        ret, pad = copy_with_padding(A,it)
+
+        szs = collect(size(ret))
+        strds = collect(strides(ret))
 
         for dim in 1:N
             n = szs[dim]
             szs[dim] = 1
 
-            M = prefiltering_system_matrix(eltype(A), n, it)
+            M, b = prefiltering_system(eltype(A), n, it)
 
             @nloops N i d->1:szs[d] begin
                 cc = @ntuple N i
-                strt = 1 + sum([(cc[i]-1)*strds[i] for i in 1:length(cc)])
+
+                strt_diff = sum([(cc[i]-1)*strds[i] for i in 1:length(cc)])
+                strt = 1 + strt_diff
                 rng = range(strt, strds[dim], n)
-                ret[rng] = M \ vec(A[rng])
-            end            
+
+                bdiff = ret[rng]
+                b += bdiff
+                ret[rng] = M \ b
+                b -= bdiff
+            end
             szs[dim] = n
         end
         ret

src/constant.jl

@@ -1,24 +1,22 @@
 type ConstantDegree <: Degree{0} end
 type Constant{GR<:GridRepresentation} <: InterpolationType{ConstantDegree,None,GR} end
-
 Constant{GR<:GridRepresentation}(::GR) = Constant{GR}()
 
-function bc_gen{IT<:Constant}(::IT, N)
-    quote
-        @nexprs $N d->(ix_d = iround(x_d))
-    end
+function bc_gen(::Constant{OnGrid}, N)
+    :(@nexprs $N d->(ix_d = round(Int, x_d)))
+end
+function bc_gen(::Constant{OnCell}, N)
+    :(@nexprs $N d->(ix_d = clamp(round(Int,x_d), 1, size(itp,d))))
 end
 
 function indices(::Constant, N)
-    quote
-        # Constant interpolation doesn't need an fx_d
-    end
+    # Constant interpolation doesn't need an fx_d
+    # but we still need to return a quote
+    :()
 end
 
 function coefficients(::Constant, N)
-    quote
-        @nexprs $N d->(c_d = one(typeof(x_d)))
-    end
+    :(@nexprs $N d->(c_d = one(typeof(x_d))))
 end
 
 function index_gen(degree::ConstantDegree, N::Integer, offsets...)

src/cubic.jl

@@ -0,0 +1,13 @@
+# This assumes integral values ixm_d, ix_d, ixp_d, and ixpp_d (typically ixm_d = ix_d-1, ixp_d = ix_d+1, except at boundaries),
+# coefficients cm_d, c_d, cp_d, and cpp_d, and an array itp.coefs
+function index_gen(::Type{Cubic}, N::Integer, offsets...)
+    if length(offsets) < N
+        d = length(offsets)+1
+        symm, sym, symp, sympp =  symbol(string("cm_",d)), symbol(string("c_",d)), symbol(string("cp_",d)), symbol(string("cpp_",d))
+        return :($symm * $(index_gen(Cubic, N, offsets...,-1)) + $sym * $(index_gen(Cubic, N, offsets..., 0)) +
+                 $symp * $(index_gen(Cubic, N, offsets..., 1)) + $sympp * $(index_gen(Cubic, N, offsets..., 2)))
+    else
+        indices = [offsetsym(offsets[d], d) for d = 1:N]
+        return :(itp.coefs[$(indices...)])
+    end
+end

src/extrapolation.jl

@@ -1,27 +1,119 @@
-
 abstract ExtrapolationBehavior
-type ExtrapError <: ExtrapolationBehavior end
 
+type ExtrapError <: ExtrapolationBehavior end
 function extrap_gen(::OnGrid, ::ExtrapError, N)
     quote
         @nexprs $N d->(1 <= x_d <= size(itp,d) || throw(BoundsError()))
     end
 end
-extrap_gen(::OnCell, e::ExtrapError, N) = extrap_gen(OnGrid(), e, N)
+function extrap_gen(::OnCell, ::ExtrapError, N)
+    quote
+        @nexprs $N d->(.5 <= x_d <= size(itp,d) + .5 || throw(BoundsError()))
+    end
+end
 
 type ExtrapNaN <: ExtrapolationBehavior end
-
 function extrap_gen(::OnGrid, ::ExtrapNaN, N)
     quote
         @nexprs $N d->(1 <= x_d <= size(itp,d) || return convert(T, NaN))
     end
 end
-extrap_gen(::OnCell, e::ExtrapNaN, N) = extrap_gen(OnGrid(), e, N)
+function extrap_gen(::OnCell, ::ExtrapNaN, N)
+    quote
+        @nexprs $N d->(.5 <= x_d <= size(itp,d) + .5 || return convert(T, NaN))
+    end
+end
 
 type ExtrapConstant <: ExtrapolationBehavior end
 function extrap_gen(::OnGrid, ::ExtrapConstant, N)
     quote
         @nexprs $N d->(x_d = clamp(x_d, 1, size(itp,d)))
     end
 end
-extrap_gen(::OnCell, e::ExtrapConstant, N) = extrap_gen(OnGrid(), e, N)
+function extrap_gen(::OnCell, ::ExtrapConstant, N)
+    quote
+        @nexprs $N d->(x_d = clamp(x_d, .5, size(itp,d)+.5))
+    end
+end
+
+type ExtrapReflect <: ExtrapolationBehavior end
+function extrap_gen(::OnGrid, ::ExtrapReflect, N)
+    quote
+        @nexprs $N d->begin
+            # translate x_d to inside the domain, and count the translations
+            ntransl = 0
+            while x_d < 1
+                x_d += size(itp, d) - 1
+                ntransl += 1
+            end
+            while x_d > size(itp, d)
+                x_d -= size(itp, d) - 1
+                ntransl += 1
+            end
+
+            # if odd number of translations, also reflect inside the domain
+            if ntransl > 0 && mod(ntransl, 2) != 0
+                x_d = size(itp, d) + 1 - x_d
+            end
+        end
+    end
+end
+function extrap_gen(::OnCell, ::ExtrapReflect, N)
+    quote
+        @nexprs $N d->begin
+            # translate x_d to inside the domain, and count the translations
+            ntransl = 0
+            while x_d < .5
+                x_d += size(itp,d)
+                ntransl += 1
+            end
+            while x_d > size(itp,d) + .5
+                x_d -= size(itp,d)
+                ntransl += 1
+            end
+
+            # if odd number of translations, also reflect inside the domain
+            if ntransl > 0 && mod(ntransl, 2) != 0
+                x_d = size(itp, d) + 1 - x_d
+            end 
+        end
+    end
+end
+
+type ExtrapPeriodic <: ExtrapolationBehavior end
+function extrap_gen(::GridRepresentation, ::ExtrapPeriodic, N)
+    quote
+        @nexprs $N d->begin
+            # translate x_d to inside the domain
+            n = convert(typeof(x_d), size(itp,d))
+            while x_d < .5
+                x_d += n
+            end
+            while x_d >= n + .5
+                x_d -= n
+            end
+        end
+    end
+end
+
+type ExtrapLinear <: ExtrapolationBehavior end
+function extrap_gen(::OnGrid, ::ExtrapLinear, N)
+    quote
+        @nexprs $N d->begin
+            if x_d < 1
+                fx_d = x_d - convert(typeof(x_d), 1)
+
+                k = -4*itp.coefs[1]
+                return itp[1] - k * fx_d
+            end
+            if x_d > size(itp, d)
+                s_d = size(itp,d)
+                fx_d = x_d - convert(typeof(x_d), s_d)
+
+                k = itp[s_d] - itp[s_d - 1]
+                return itp[s_d] + k * fx_d
+            end
+        end
+    end
+end
+extrap_gen(::OnCell, e::ExtrapLinear, N) = extrap_gen(OnGrid(), e, N)

src/linear.jl

@@ -2,17 +2,11 @@ type LinearDegree <: Degree{1} end
 type Linear{GR<:GridRepresentation} <: InterpolationType{LinearDegree,None,GR} end
 Linear{GR<:GridRepresentation}(::GR) = Linear{GR}()
 
-Interpolation{T,N,EB<:ExtrapolationBehavior}(A::Array{T,N}, ::Linear, ::EB) = Interpolation{T,N,Linear{OnGrid},EB}(A)
-
-function bc_gen(::Linear{OnGrid}, N)
-    quote
-        # LinearOnGrid's boundary condition doesn't require any action;
-        # just return the cell index of the interpolation coordinate
-        @nexprs $N d->(ix_d = @compat floor(Integer, x_d))
-    end
+function bc_gen(::Linear, N)
+    :(@nexprs $N d->(ix_d = clamp(floor(Int,x_d), 1, size(itp,d)-1)))
 end
 
-function indices(::Linear{OnGrid}, N)
+function indices(::Linear, N)
     quote
         # fx_d is a parameter in [0,1] such that x_d = ix_d + fx_d
         @nexprs $N d->(fx_d = x_d - convert(typeof(x_d), ix_d))
@@ -21,7 +15,7 @@ function indices(::Linear{OnGrid}, N)
     end
 end
 
-function coefficients(::Linear{OnGrid}, N)
+function coefficients(::Linear, N)
     quote
         @nexprs $N d->begin
             c_d = one(typeof(fx_d)) - fx_d