When interpolating Float16 data, the interpolator reports its element type as Float64 even if it correctly returns Float16 data when called.
julia> using Interpolations
julia> itp = interpolate(rand(Float16, 3, 3), BSpline(Linear()));
julia> isa(itp, Interpolations.BSplineInterpolation{Float16})
false
julia> eltype(itp)
Float64
julia> itp(Float16(1.1), Float16(1.1))
Float16(0.9077)
It looks to be due to the tweight function defaulting to Float64 and I can apparently solve the problem by specializing tweight for Float16 or changing its default value to be the same type as the input array. Since I am not very familiar with the codebase I am not sure if that's all it should be done.