# Using JuliaMath/Interpolations.jl on an image

**URL:** <https://discourse.julialang.org/t/using-juliamath-interpolations-jl-on-an-image/7113>\
**Category:** General Usage\
**Created:** [November 17, 2017, 1:38am UTC](https://discourse.julialang.org/t/using-juliamath-interpolations-jl-on-an-image/7113 "2017-11-17T01:38:30Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![mcbaron](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mcbaron/32/4027_2.png) [@mcbaron](https://discourse.julialang.org/u/mcbaron)\
**Post date:** [November 17, 2017, 1:38am UTC](https://discourse.julialang.org/t/using-juliamath-interpolations-jl-on-an-image/7113/1 "2017-11-17T01:38:30Z")

</div>

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](https://gist.github.com/mcbaron/a5a730513efd4b4b923a1b009728533b)

dx and dy are 1320x1000 arrays that I’ve generated (via Beier-Neely)[[Beier–Neely morphing algorithm - Wikipedia](https://en.wikipedia.org/wiki/Beier%E2%80%93Neely_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,

---

<div class="post-metadata">

**Author:** ![tim.holy](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tim.holy/32/52_2.png) [@tim.holy](https://discourse.julialang.org/u/tim.holy)\
**Post date:** [November 18, 2017, 4:56pm UTC](https://discourse.julialang.org/t/using-juliamath-interpolations-jl-on-an-image/7113/2 "2017-11-18T16:56:51Z")

</div>

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:

```julia
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](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.

---

<div class="post-metadata">

**Author:** ![mcbaron](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mcbaron/32/4027_2.png) [@mcbaron](https://discourse.julialang.org/u/mcbaron)\
**Post date:** [November 25, 2017, 9:45pm UTC](https://discourse.julialang.org/t/using-juliamath-interpolations-jl-on-an-image/7113/3 "2017-11-25T21:45:23Z")

</div>

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

The same way that I did, or differently?

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [November 26, 2017, 2:16am UTC](https://discourse.julialang.org/t/using-juliamath-interpolations-jl-on-an-image/7113/4 "2017-11-26T02:16:44Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![tim.holy](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tim.holy/32/52_2.png) [@tim.holy](https://discourse.julialang.org/u/tim.holy)\
**Post date:** [November 26, 2017, 10:31am UTC](https://discourse.julialang.org/t/using-juliamath-interpolations-jl-on-an-image/7113/5 "2017-11-26T10:31:23Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![mcbaron](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mcbaron/32/4027_2.png) [@mcbaron](https://discourse.julialang.org/u/mcbaron)\
**Post date:** [November 26, 2017, 6:11pm UTC](https://discourse.julialang.org/t/using-juliamath-interpolations-jl-on-an-image/7113/6 "2017-11-26T18:11:54Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [November 26, 2017, 9:49pm UTC](https://discourse.julialang.org/t/using-juliamath-interpolations-jl-on-an-image/7113/7 "2017-11-26T21:49:45Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Kunz-David](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kunz-david/32/24888_2.png) [@Kunz-David](https://discourse.julialang.org/u/Kunz-David)\
**Post date:** [April 4, 2020, 10:31pm UTC](https://discourse.julialang.org/t/using-juliamath-interpolations-jl-on-an-image/7113/8 "2020-04-04T22:31:36Z")

</div>

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:

```julia
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)

```

---

<div class="post-metadata">

**Author:** ![momo1](https://avatars.discourse-cdn.com/v4/letter/m/bbe5ce/32.png) [@momo1](https://discourse.julialang.org/u/momo1)\
**Post date:** [September 19, 2022, 6:18am UTC](https://discourse.julialang.org/t/using-juliamath-interpolations-jl-on-an-image/7113/9 "2022-09-19T06:18:39Z")

</div>

That gets an error now imgw = imWarp(img, dx, dy);

> ERROR: UndefVarError: indices\_spatial not defined  
> Stacktrace:  
> [1] imWarp(img::Matrix{ColorTypes.RGB{FixedPointNumbers.N0f8}}, dx::Matrix{Float64}, dy::Matrix{Float64})  
> @ Main ./REPL[62]:3  
> [2] top-level scope  
> @ REPL[69]:1

So using ImageCore will get that to work