Interpolations uses the same rules for
itp[i, j] that Julia (and Matlab, which this seems to mimic) uses for general indexing: when
j are vectors/arrays,
A[i, j] computes the result over the cartesian product
i⊗j. In this case that would mean you’re trying to compute more than 10^12 points, which is why it runs out of memory.
Just as a tip, idiomatic Julia+Images+Interpolations makes this function considerably simpler and faster than the matlab version:
using Interpolations, ColorVectorSpace
function imWarp(img, dx, dy)
itp = interpolate(img, BSpline(Linear()), OnGrid())
inds = indices(img)
rng = extrema.(inds)
imw = similar(img, eltype(itp))
for I in CartesianRange(inds)
dxi, dyi = dx[I], dy[I]
y, x = clamp(I+dyi, rng...), clamp(I+dxi, rng...)
imw[I] = itp[y, x]
## Usage demo
using TestImages, ImageView
img = testimage("mandrill");
# Let's use a random warp. We blur it with a Gaussian kernel to make it reasonably smooth.
kern = KernelFactors.IIRGaussian((10,10)) # IIRGaussian is fast even for very large σ
dx, dy = imfilter(100*randn(size(img)), kern), imfilter(100*randn(size(img)), kern);
imgw = imWarp(img, dx, dy);
This implementation of
imWarp has many advantages:
- it’s much shorter (it doesn’t need
meshgrid at all)
- it works for color images as well as grayscale images (you don’t have to have a separate color channel, see http://juliaimages.github.io/latest/arrays_colors.html)
- for interpolation, the slow step is usually computation of the coefficients. Here you only have to do that once per pixel even for color images. If you use the channel-based approach, for RGB images you have to do that 3 times per pixel.
- you don’t need to create all those temporary arrays
- once you understand the syntax, it’s much more obvious what this does (in particular it avoids the trap that your version fell into)
Vectorization can sometimes be nice, but the fact that Matlab requires vectorization in order to get decent performance often seems to make code a bit torturous.