I’m struggling with an OutOfMemory() error while trying to index into an interpolation object with an image size set of indices. I know enough to think that this means that there is probably something wrong with my implementation.
dx and dy are 1320x1000 arrays that I’ve generated (via Beier-Neely)[Beier–Neely morphing algorithm - Wikipedia] as maps for warping a 1320x1000 pixel image. Each entry defines an offset that that particular pixel needs to “move” in the x and y directions respectfully.
I use these offsets and a meshgrid to figure out where I need to sample the interpolated image’s surface, and this is where I run out of memory.
I’m sure I’ve just missed something, does anyone have experience here?
Interpolations uses the same rules for itp[i, j] that Julia (and Matlab, which this seems to mimic) uses for general indexing: when i and 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[1]+dyi, rng[1]...), clamp(I[2]+dxi, rng[2]...)
imw[I] = itp[y, x]
end
return imw
end
## Usage demo
using TestImages, ImageView
img = testimage("mandrill");
imshow(img)
# Let's use a random warp. We blur it with a Gaussian kernel to make it reasonably smooth.
using ImageFiltering
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);
imshow(imgw)
This implementation of imWarp has many advantages:
it’s much shorter (it doesn’t need meshgrid at all)
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.
Sorry @mcbaron, I had some typos in the above (and probably knots hanging around as a global variable). I edited the post and included a runnable example.
There is no such thing Many important contributions started with simple implementations. The Images.jl maintainers are always willing to give good feedback. They are awesome.
All you need to do is provide an initial implementation that works with simple Julia arrays. Later, people will help you wrap the arrays in the Images.jl types. That is another beautiful thing about multiple dispatch, start simple and build up types later.
Some time has passed since @tim.holy posted his original answer and since I wanted to use the code, I had to update it to work with Julia and its modules as of now.
Here it is updated without deprecated code if its useful to anyone:
using Interpolations, ColorVectorSpace, Images
function imWarp(img, dx, dy)
itp = interpolate(img, BSpline(Linear()))
inds = indices_spatial(img)
rng = extrema.(inds)
imw = similar(img, eltype(itp))
for I in CartesianIndices(inds)
dxi, dyi = dx[I], dy[I]
y, x = clamp(I[1]+dyi, rng[1]...), clamp(I[2]+dxi, rng[2]...)
imw[I] = itp(y, x)
end
return imw
end
# Usage demo
using TestImages, ImageView
img = testimage("mandril");
imshow(img)
using ImageFiltering
# 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);
imshow(imgw)
