Using JuliaMath/Interpolations.jl on an image


#1

Hello all,

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.

I was hoping you could give me a pointer or two.
What’s going on:
I’ve thrown my code into a gist: https://gist.github.com/mcbaron/a5a730513efd4b4b923a1b009728533b

dx and dy are 1320x1000 arrays that I’ve generated (via Beier-Neely)[https://en.wikipedia.org/wiki/Beier–Neely_morphing_algorithm] 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?

Thanks,


#2

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)
  • 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.


#3

How are you defining knots, which is passed to interpolate, in this context?

The same way that I did, or differently?


#4

Very cool @mcbaron, by the way, are you planning to contribute the Beier-Neely algorithm to the Images.jl project at some point? :slight_smile:


#5

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.


#6

I’d love to!

My implementation is definitely not up to the caliber of the Images.jl package, and it’d have to be made more accepting and aware of types.


#7

There is no such thing :slight_smile: 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.