For a Matrix{T} as input array, the function map_coordinates
is defined as follows:
using Interpolations
import StaticArrays: SVector
function map_coordinates(input::Matrix{T}, coordinates::Vector{SVector{N, T}}; method=BSpline(Linear())) where {N,T}
output = T[]
itp = interpolate(input, method, OnCell())
for coord in coordinates
push!(output, itp[coord...])
end
output
end
input
is a Matrix{T} of size (m,n), and its elements
are the the values of a real function, f, defined on the planar rectangle [1, m] x [1,n], at the integer coordinates (i, j), i=1…m, j=1, …n.
i.e. input[i,j]=f(i, j)
.
coordinates
is a vector of planar coordinates (x,y), with x\in[1,m], y\in [1, n].
map_coordinates
evaluates the function f at these coordinates and pushes it values to output.
const Sv=SVector
input = [1.65 -1.35 -0.98;
2.4 1.2 0.7;
4.3 3.36 2.54;
3.2 2.86 1.74]
Let us get the values at the indices for second column
julia> map_coordinates(input, [Sv(1.0,2), Sv(2,2.0), Sv(3,2), Sv(4,2)])
4-element Vector{Float64}:
-1.35
1.2
3.36
2.86
Retrieve the values at coordinates inside the rectangle on which the function is defined:
julia> coordinates=[Sv(1.1, 1.5), Sv(1.2, 1.89), Sv(3,2)]
julia> map_coordinates(input, coordinates)
3-element Vector{Float64}:
0.31500000000000006
-0.5495999999999999
3.36
If at least a coordinate is outside this rectangle we get an error. The scipy function assigns a prescribed value, cval=0, to such coordinates. (we can add this to the above definition) .
julia> map_coordinates(input, [Sv(0.97, 1.2)])
ERROR: BoundsError: attempt to access 4×3 interpolate(::Matrix{Float64}, BSpline(Linear())) with element type Float64 at index [0.97, 1.2]