Skip to content

Commit cfb4d13

Browse files
mkittimkitti
authored
Merge pull request #338 from ohmsweetohm1/patch-1
Add interp1-like interpolation modes
2 parents 7719e98 + cd2a812 commit cfb4d13

File tree

8 files changed

+154
-42
lines changed

8 files changed

+154
-42
lines changed

‎perf/benchmarks.jl‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -30,14 +30,14 @@ end
3030
strip_prefix(str::AbstractString) = replace(str, "Interpolations."=>"")
3131
benchstr(::Type{T}) where {T<:Interpolations.GridType} = strip_prefix(string(T))
3232

33-
benchstr(::Type{Constant}) = "Constant()"
33+
benchstr(::Type{<:Constant}) = "Constant()"
3434
benchstr(::Type{Linear}) = "Linear()"
3535
benchstr(::Type{Quadratic{BC}}, ::Type{GT}) where {BC<:Interpolations.BoundaryCondition,GT<:Interpolations.GridType} =
3636
string("Quadratic(", strip_prefix(string(BC)), "(", strip_prefix(string(GT)), "()))")
3737
benchstr(::Type{Cubic{BC}}, ::Type{GT}) where {BC<:Interpolations.BoundaryCondition,GT<:Interpolations.GridType} =
3838
string("Cubic(", strip_prefix(string(BC)), "(", strip_prefix(string(GT)), "()))")
3939

40-
groupstr(::Type{Constant}) = "constant"
40+
groupstr(::Type{<:Constant}) = "constant"
4141
groupstr(::Type{Linear}) = "linear"
4242
groupstr(::Type{Quadratic}) = "quadratic"
4343
groupstr(::Type{Cubic}) = "cubic"

‎src/b-splines/b-splines.jl‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -252,6 +252,6 @@ include("indexing.jl")
252252
include("prefiltering.jl")
253253
include("../filter1d.jl")
254254

255-
Base.parent(A::BSplineInterpolation{T,N,TCoefs,UT}) where {T,N,TCoefs,UT<:Union{BSpline{Linear},BSpline{Constant}}} = A.coefs
255+
Base.parent(A::BSplineInterpolation{T,N,TCoefs,UT}) where {T,N,TCoefs,UT<:Union{BSpline{Linear},BSpline{<:Constant}}} = A.coefs
256256
Base.parent(A::BSplineInterpolation{T,N,TCoefs,UT}) where {T,N,TCoefs,UT} =
257257
throw(ArgumentError("The given BSplineInterpolation does not serve as a \"view\" for a parent array. This would only be true for Constant and Linear b-splines."))

‎src/b-splines/constant.jl‎

Lines changed: 24 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1,12 +1,33 @@
1-
struct Constant <: Degree{0} end
1+
export Nearest, Previous, Next
2+
3+
abstract type ConstantInterpType end
4+
struct Nearest <: ConstantInterpType end
5+
struct Previous <: ConstantInterpType end
6+
struct Next <: ConstantInterpType end
7+
8+
struct Constant{T<:ConstantInterpType} <: Degree{0}
9+
Constant() = new{Nearest}()
10+
Constant{Previous}() = new{Previous}()
11+
Constant{Next}() = new{Next}()
12+
end
213

314
"""
415
Constant b-splines are *nearest-neighbor* interpolations, and effectively
516
return `A[round(Int,x)]` when interpolating.
617
"""
718
Constant
819

9-
function positions(::Constant, ax, x) # discontinuity occurs at half-integer locations
20+
function positions(c::Constant{Previous}, ax, x) # discontinuity occurs at integer locations
21+
xm = floorbounds(x, ax)
22+
δx = x - xm
23+
fast_trunc(Int, xm), δx
24+
end
25+
function positions(c::Constant{Next}, ax, x) # discontinuity occurs at integer locations
26+
xm = ceilbounds(x, ax)
27+
δx = x - xm
28+
fast_trunc(Int, xm), δx
29+
end
30+
function positions(c::Constant{Nearest}, ax, x) # discontinuity occurs at half-integer locations
1031
xm = roundbounds(x, ax)
1132
δx = x - xm
1233
fast_trunc(Int, xm), δx
@@ -16,4 +37,4 @@ value_weights(::Constant, δx) = (1,)
1637
gradient_weights(::Constant, δx) = (0,)
1738
hessian_weights(::Constant, δx) = (0,)
1839

19-
padded_axis(ax::AbstractUnitRange, ::BSpline{Constant}) = ax
40+
padded_axis(ax::AbstractUnitRange, ::BSpline{<:Constant}) = ax

‎src/b-splines/indexing.jl‎

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -199,6 +199,16 @@ function _floorbounds(x, ax::Union{Tuple{Real,Real}, AbstractUnitRange})
199199
ifelse(x < l, floor(x+h), floor(x+zero(h)))
200200
end
201201

202+
ceilbounds(x::Integer, ax::Tuple{Real,Real}) = x
203+
ceilbounds(x::Integer, ax::AbstractUnitRange) = x
204+
ceilbounds(x, ax::Tuple{Real,Real}) = _ceilbounds(x, ax)
205+
ceilbounds(x, ax::AbstractUnitRange) = _ceilbounds(x, ax)
206+
function _ceilbounds(x, ax::Union{Tuple{Real,Real}, AbstractUnitRange})
207+
u = last(ax)
208+
h = half(x)
209+
ifelse(x > u, ceil(x+h), ceil(x+zero(h)))
210+
end
211+
202212
half(x) = oneunit(x)/2
203213

204214
symmatrix(h::NTuple{1,Any}) = SMatrix{1,1}(h)

‎src/gridded/indexing.jl‎

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -40,6 +40,14 @@ function gridded_floorbounds(x, knotvec::AbstractVector)
4040
max(i, first(axes1(knotvec)))
4141
end
4242

43+
ceilbounds(x::Integer, knotvec::AbstractVector) = gridded_ceilbounds(x, knotvec)
44+
ceilbounds(x, knotvec::AbstractVector) = gridded_ceilbounds(x, knotvec)
45+
function gridded_ceilbounds(x, knotvec::AbstractVector)
46+
i = find_knot_index(knotvec, x)
47+
iclamp = max(i, first(axes1(knotvec)))
48+
min(iclamp+1, last(axes1(knotvec)))
49+
end
50+
4351
@inline find_knot_index(knotv, x) = searchsortedfirst(knotv, x, Base.Order.ForwardOrdering()) - 1
4452
@inline find_knot_index(knotv, x::AbstractVector) = searchsortedfirst_vec(knotv, x) .- 1
4553

‎test/b-splines/constant.jl‎

Lines changed: 99 additions & 26 deletions
Original file line numberDiff line numberDiff line change
@@ -5,35 +5,108 @@
55
A2 = rand(Float64, N1, N1) * 100
66
A3 = rand(Float64, N1, N1, N1) * 100
77

8-
for (constructor, copier) in ((interpolate, x->x), (interpolate!, copy))
9-
isinplace = constructor == interpolate!
10-
itp1 = @inferred(constructor(copier(A1), BSpline(Constant())))
11-
itp2 = @inferred(constructor(copier(A2), BSpline(Constant())))
12-
itp3 = @inferred(constructor(copier(A3), BSpline(Constant())))
13-
14-
@test parent(itp1) === itp1.coefs
15-
@test Interpolations.lbounds(itp1) == (1,)
16-
@test Interpolations.ubounds(itp1) == (N1,)
17-
18-
# Evaluation on provided data points
19-
for (itp, A) in ((itp1, A1), (itp2, A2), (itp3, A3))
20-
check_axes(itp, A, isinplace)
21-
check_inbounds_values(itp, A)
22-
check_oob(itp)
23-
can_eval_near_boundaries(itp)
24-
end
8+
@testset "Constant{Nearest}" begin
9+
for (constructor, copier) in ((interpolate, x -> x), (interpolate!, copy))
10+
isinplace = constructor == interpolate!
11+
itp1 = @inferred(constructor(copier(A1), BSpline(Constant())))
12+
itp2 = @inferred(constructor(copier(A2), BSpline(Constant())))
13+
itp3 = @inferred(constructor(copier(A3), BSpline(Constant())))
14+
15+
@test parent(itp1) === itp1.coefs
16+
@test Interpolations.lbounds(itp1) == (1,)
17+
@test Interpolations.ubounds(itp1) == (N1,)
2518

26-
# Evaluation between data points (tests constancy)
27-
for i in 2:N1-1
28-
@test A1[i] == itp1(i+.3) == itp1(i+.3) == itp1(i-.3) == itp1(i-.3)
19+
# Evaluation on provided data points
20+
for (itp, A) in ((itp1, A1), (itp2, A2), (itp3, A3))
21+
check_axes(itp, A, isinplace)
22+
check_inbounds_values(itp, A)
23+
check_oob(itp)
24+
can_eval_near_boundaries(itp)
25+
end
26+
27+
# Evaluation between data points (tests constancy)
28+
for i in 2:N1 - 1
29+
@test A1[i] == itp1(i + .3) == itp1(i - .3)
30+
end
31+
# 2D
32+
for i in 2:N1 - 1, j in 2:N1 - 1
33+
@test A2[i,j] == itp2(i + .4, j - .3) == itp2(i - .4, j + .3)
34+
end
35+
# 3D
36+
for i in 2:N1 - 1, j in 2:N1 - 1, k in 2:N1 - 1
37+
@test A3[i,j,k] == itp3(i + .4, j - .3, k + .1) == itp3(i + .4, j - .3, k + .2)
38+
end
2939
end
30-
# 2D
31-
for i in 2:N1-1, j in 2:N1-1
32-
@test A2[i,j] == itp2(i+.4,j-.3) == itp2(i+.4,j-.3)
40+
end
41+
42+
@testset "Constant{Previous}" begin
43+
for (constructor, copier) in ((interpolate, x -> x), (interpolate!, copy))
44+
isinplace = constructor == interpolate!
45+
itp1 = @inferred(constructor(copier(A1), BSpline(Constant{Previous}())))
46+
itp2 = @inferred(constructor(copier(A2), BSpline(Constant{Previous}())))
47+
itp3 = @inferred(constructor(copier(A3), BSpline(Constant{Previous}())))
48+
49+
@test parent(itp1) === itp1.coefs
50+
@test Interpolations.lbounds(itp1) == (1,)
51+
@test Interpolations.ubounds(itp1) == (N1,)
52+
53+
# Evaluation on provided data points
54+
for (itp, A) in ((itp1, A1), (itp2, A2), (itp3, A3))
55+
check_axes(itp, A, isinplace)
56+
check_inbounds_values(itp, A)
57+
check_oob(itp)
58+
can_eval_near_boundaries(itp)
59+
end
60+
61+
# Evaluation between data points (tests constancy)
62+
for i in 2:N1 - 1
63+
@test A1[i] == itp1(i + .3) == itp1(i + .6)
64+
@test A1[i - 1] == itp1(i - .3) == itp1(i - .6)
65+
end
66+
# 2D
67+
for i in 2:N1 - 1, j in 2:N1 - 1
68+
@test A2[i,j - 1] == itp2(i + .4, j - .3) == itp2(i + .4, j - .6)
69+
end
70+
# 3D
71+
for i in 2:N1 - 1, j in 2:N1 - 1, k in 2:N1 - 1
72+
@test A3[i,j - 1,k] == itp3(i + .4, j - .3, k + .1) == itp3(i + .4, j - .3, k + .2)
73+
end
3374
end
34-
# 3D
35-
for i in 2:N1-1, j in 2:N1-1, k in 2:N1-1
36-
@test A3[i,j,k] == itp3(i+.4,j-.3,k+.1) == itp3(i+.4,j-.3,k+.2)
75+
end
76+
77+
@testset "Constant{Next}" begin
78+
for (constructor, copier) in ((interpolate, x -> x), (interpolate!, copy))
79+
isinplace = constructor == interpolate!
80+
itp1 = @inferred(constructor(copier(A1), BSpline(Constant{Next}())))
81+
itp2 = @inferred(constructor(copier(A2), BSpline(Constant{Next}())))
82+
itp3 = @inferred(constructor(copier(A3), BSpline(Constant{Next}())))
83+
84+
@test parent(itp1) === itp1.coefs
85+
@test Interpolations.lbounds(itp1) == (1,)
86+
@test Interpolations.ubounds(itp1) == (N1,)
87+
88+
# Evaluation on provided data points
89+
for (itp, A) in ((itp1, A1), (itp2, A2), (itp3, A3))
90+
check_axes(itp, A, isinplace)
91+
check_inbounds_values(itp, A)
92+
check_oob(itp)
93+
can_eval_near_boundaries(itp)
94+
end
95+
96+
# Evaluation between data points (tests constancy)
97+
for i in 2:N1 - 1
98+
@test A1[i + 1] == itp1(i + .3) == itp1(i + .6)
99+
@test A1[i] == itp1(i - .3) == itp1(i - .6)
100+
end
101+
# 2D
102+
for i in 2:N1 - 1, j in 2:N1 - 1
103+
@test A2[i + 1,j] == itp2(i + .4, j - .3) == itp2(i + .4, j - .6)
104+
end
105+
# 3D
106+
for i in 2:N1 - 1, j in 2:N1 - 1, k in 2:N1 - 1
107+
@test A3[i + 1,j,k + 1] == itp3(i + .4, j - .3, k + .1) == itp3(i + .4, j - .3, k + .2)
108+
end
37109
end
38110
end
39111
end
112+

‎test/io.jl‎

Lines changed: 6 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -8,12 +8,12 @@ using Test
88
A = rand(8,20)
99

1010
itp = interpolate(A, BSpline(Constant()))
11-
@test summary(itp) == "8×20 interpolate(::Matrix{Float64}, BSpline(Constant())) with element type Float64" ||
12-
summary(itp) == "8×20 interpolate(::Array{Float64,2}, BSpline(Constant())) with element type Float64"
11+
@test summary(itp) == "8×20 interpolate(::Matrix{Float64}, BSpline(Constant{Nearest}())) with element type Float64" ||
12+
summary(itp) == "8×20 interpolate(::Array{Float64,2}, BSpline(Constant{Nearest}())) with element type Float64"
1313

1414
itp = interpolate(A, BSpline(Constant()))
15-
@test summary(itp) == "8×20 interpolate(::Matrix{Float64}, BSpline(Constant())) with element type Float64" ||
16-
summary(itp) == "8×20 interpolate(::Array{Float64,2}, BSpline(Constant())) with element type Float64"
15+
@test summary(itp) == "8×20 interpolate(::Matrix{Float64}, BSpline(Constant{Nearest}())) with element type Float64" ||
16+
summary(itp) == "8×20 interpolate(::Array{Float64,2}, BSpline(Constant{Nearest}())) with element type Float64"
1717

1818
itp = interpolate(A, BSpline(Linear()))
1919
@test summary(itp) == "8×20 interpolate(::Matrix{Float64}, BSpline(Linear())) with element type Float64" ||
@@ -47,8 +47,8 @@ using Test
4747
summary(itp) == "8×20 interpolate((::Array{Int64,1},::Array{Float64,1}), ::Array{Float64,2}, Gridded(Linear())) with element type Float64"
4848

4949
itp = interpolate(knots, A, (Gridded(Linear()),Gridded(Constant())))
50-
@test summary(itp) == "8×20 interpolate((::Vector{Int64},::Vector{Float64}), ::Matrix{Float64}, (Gridded(Linear()), Gridded(Constant()))) with element type Float64" ||
51-
summary(itp) == "8×20 interpolate((::Array{Int64,1},::Array{Float64,1}), ::Array{Float64,2}, (Gridded(Linear()), Gridded(Constant()))) with element type Float64"
50+
@test summary(itp) == "8×20 interpolate((::Vector{Int64},::Vector{Float64}), ::Matrix{Float64}, (Gridded(Linear()), Gridded(Constant{Nearest}()))) with element type Float64" ||
51+
summary(itp) == "8×20 interpolate((::Array{Int64,1},::Array{Float64,1}), ::Array{Float64,2}, (Gridded(Linear()), Gridded(Constant{Nearest}()))) with element type Float64"
5252

5353
# issue #260
5454
A = (1:4)/4

‎test/type-instantiation.jl‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -68,10 +68,10 @@ for T in (
6868
end
6969
# sample one dimspec for each interpolation constructor to see that it
7070
# calls the construct_instance() function correctly
71-
@inferred(interpolate(A2, Tuple{BSpline{Linear},BSpline{Constant}}, Tuple{OnGrid,OnCell}))
72-
@inferred(interpolate!(copy(A2), Tuple{BSpline{Linear},BSpline{Constant}}, Tuple{OnGrid,OnCell}))
73-
@inferred(interpolate(A2, Tuple{BSpline{Linear},BSpline{Constant}}, (OnGrid(),OnCell())))
74-
@inferred(interpolate!(copy(A2), Tuple{BSpline{Linear},BSpline{Constant}}, (OnGrid(),OnCell())))
71+
@inferred(interpolate(A2, Tuple{BSpline{Linear},BSpline{<:Constant}}, Tuple{OnGrid,OnCell}))
72+
@inferred(interpolate!(copy(A2), Tuple{BSpline{Linear},BSpline{<:Constant}}, Tuple{OnGrid,OnCell}))
73+
@inferred(interpolate(A2, Tuple{BSpline{Linear},BSpline{<:Constant}}, (OnGrid(),OnCell())))
74+
@inferred(interpolate!(copy(A2), Tuple{BSpline{Linear},BSpline{<:Constant}}, (OnGrid(),OnCell())))
7575
@inferred(interpolate(A2, (BSpline(Linear()),BSpline(Constant())), Tuple{OnGrid,OnCell}))
7676
@inferred(interpolate!(copy(A2), (BSpline(Linear()),BSpline(Constant())), Tuple{OnGrid,OnCell}))
7777

0 commit comments

Comments
 (0)