image.jl

eval_neumann(img, x) = img[clamp.(Tuple(x), axes(img))...]

function interpolate_bilinear(img, x)
    x0 = floor.(Int, x)
    x1 = x0 .+ 1

    val = zero(eltype(img))
    for idx in CartesianIndices(ntuple(_->0:1, ndims(img)))
	cornerbool = Bool.(Tuple(idx))
	λ = ifelse.(cornerbool, x .- x0, x1 .- x)
	corner = ifelse.(cornerbool, x1, x0)
	val += prod(λ) * eval_neumann(img, CartesianIndex(corner))
    end
    return val
end

function warp_backwards(img, u)
    d = ndims(img)
    axes(img) == axes(u)[2:end] && d == size(u, 1) ||
	throw(ArgumentError("invalid dimensions"))

    res = similar(img)
    for I in CartesianIndices(img)
	x = Tuple(I) .+ ntuple(i -> u[i, I], d)
	res[I] = interpolate_bilinear(img, x)
    end
    return res
end