How to transform this grid?

The following question is partly relevant to camera calibrations, but is relevant to any situation where you want to transform a set of coordinates from one map to another:

example

Say I have a bunch of known coordinates, the black dots above, and some new points, the red dots.
I know that these wanky black dots should really be spaced nicely on a different coordinate system:

example1

How can I transform one system to the other, including the red dots?

1 Like

have a look at https://github.com/JuliaMath/Interpolations.jl

1 Like

I know, but here’s my problem:
In my case the input is a set of x and y pairs and the output is also a set of x and y pairs. In Interpolations.jl the output is one number (e.g. z), except in the case of Parametric splines but that doesn’t work here.
Can you see how I can use Interpolations.jl here?

Seems like a good case for Finite Element interpolation?

Just use 2 interpolants one for X one for Y.

Edit: Alternatively try using [x,y] (Vector not Tuple) as the scalar Z data. It will probably just work.

Indeed it does! Awesome. Thanks @TsurHerman!

Actually, this does not work, here’s why:

Interpolations.jl works with multidimensional knots. These knots can be 1:n [1], scaled e.g. 1.:2.:40. [2], or gridded e.g. [x^2 for x = 1:8] [3].

This means that the individual levels of the knots are independent from the levels of knots from a higher dimension. Here’s what I mean:
In the case of:

knots = ([x^2 for x = 1:8], [0.2y for y = 1:20])

x[2] is equal to 4 no matter what y is equal to. I need knots that are closer to this:

x = [i + (rand()-.5)/7 for i in 1:10, j in 1:10]
y = [j + (rand()-.5)/7 for i in 1:10, j in 1:10]
knots = (x,y)

but to quote the error message Interpolations.jl gives me:
“knot vectors must have the same number of elements as the corresponding dimension of the array”.

The reason I need such “independent” knots in this specific case is that I’m mapping from a coordinate system where each x and y are different from each other (see the black dots in the first figure of this post) to one where the x and y points are regular. In fact, if I could have done the inverse of that then the problem would have been solved. And indeed I can use Interpolations.jl to convert from the correct grid to the wanky one, but I need to go the other way around.

Hope I managed to explain myself… :pray:

I’ve been trying to figure out how to apply FEMBasis.jl to this problem. @kristoffer.carlsson, could you give me an example of how to use Finite Element interpolation in this specific case: noisy x and y pairs to regular x and ys?

Thanks!!!

Use 2 interpolants , let the knots or “interpolation points” be your wanky grid. Let Z of the first interpolant be the corrected X (the matching regular grid) and interpolate , repeat for Y.

Hmm… It seems to “work”, but I can’t help feeling this can’t be right.

I’m essentially feeding the interpolation knots that are one set of x values per y value, where each respective level of x is very similar across the levels of y:

 0.914293   1.12394   0.888184  0.939205  …   0.977576  0.950751  0.878769   0.963788
 2.07355    2.05893   2.11997   2.10999       1.91213   1.92552   1.95643    2.00277 
 2.96772    2.99958   3.07324   3.05772       3.03218   3.0367    3.03536    2.93699 
 4.10374    3.95213   3.9412    3.99084       4.0744    3.90924   3.913      3.93404 
 4.92337    5.04418   5.11515   5.03542       4.9369    4.95768   5.06762    5.08227 
 5.98328    5.95238   6.01513   6.12016   …   6.07736   5.9214    6.0        6.09871 
 6.93604    7.0379    6.97642   6.94647       7.02528   6.99787   6.95569    7.08504 
 7.92821    7.91549   8.08373   8.01549       8.06387   8.09904   8.02936    7.89755 
 8.9328     9.06075   9.022     8.94114       9.05516   8.97034   9.10086    9.03877 
 9.93574   10.0749   10.016     9.89389      10.0083    9.88933   9.98516   10.1002  

and respective points that are identical across the set of y values:

 125.0    125.0    125.0    125.0    125.0    125.0    125.0    125.0    125.0    125.0  
 144.444  144.444  144.444  144.444  144.444  144.444  144.444  144.444  144.444  144.444
 163.889  163.889  163.889  163.889  163.889  163.889  163.889  163.889  163.889  163.889
 183.333  183.333  183.333  183.333  183.333  183.333  183.333  183.333  183.333  183.333
 202.778  202.778  202.778  202.778  202.778  202.778  202.778  202.778  202.778  202.778
 222.222  222.222  222.222  222.222  222.222  222.222  222.222  222.222  222.222  222.222
 241.667  241.667  241.667  241.667  241.667  241.667  241.667  241.667  241.667  241.667
 261.111  261.111  261.111  261.111  261.111  261.111  261.111  261.111  261.111  261.111
 280.556  280.556  280.556  280.556  280.556  280.556  280.556  280.556  280.556  280.556
 300.0    300.0    300.0    300.0    300.0    300.0    300.0    300.0    300.0    300.0  

This seems very wrong because the interpolation is ignorant of what the y value is for each of those columns (and in my case also, what the y value is for each row, simply put, each cell).

Here’s the code I’m using, just to be clear:

using Plots, Interpolations
image_x = [i + (rand()-.5)/4 for i in 1:10, j in 1:10]
image_y = [j + (rand()-.5)/4 for i in 1:10, j in 1:10]
plot(image_x,image_y, marker=:circle,color=:black)
plot!(image_x',image_y',color=:black,legend=false, grid=false)
x_in = [rand(1.5:9.5) + (rand()-.5)/2 for i in 1:10]
y_in = [rand(1.5:9.5) + (rand()-.5)/2 for i in 1:10]
scatter!(x_in,y_in,marker=:circle,color=:red)
png("a1.png")
world_x = linspace(125,300,10)
world_y = linspace(-5,53,10)
knots = vec(image_x)
A = vec(world_x .+ 0world_y')
o = sortperm(knots)
itx = interpolate((knots[o],), A[o], Gridded(Linear()))
knots = vec(image_y)
A = vec(0world_x .+ world_y')
o = sortperm(knots)
ity = interpolate((knots[o],), A[o], Gridded(Linear()))
x_out = itx[x_in]
y_out = ity[y_in]
x = repmat(world_x,1,10)
y = repmat(world_y',10)
plot(x,y, marker=:circle,color=:black)
plot!(x',y',color=:black,legend=false, grid=false)
scatter!(x_out, y_out, marker=:circle,color=:red)
png("a2.png")

In my lab’s internal deformable image registration code, we have this:

"""
`Ď• = GridDeformation(u::Array{<:SVector}, dims)` creates a
deformation `Ď•` for an array of size `dims`.  `u` specifies the
"pixel-wise" displacement at a series of control points that are
evenly-spaced over the domain specified by `dims` (i.e., using
knot-vectors `linspace(1,dims[d],size(u,d))`).  In particular, each
corner of the array is the site of one control point.

`Ď• = GridDeformation(u::Array{<:SVector}, knots)` specifies the
knot-vectors manually. `u` must have dimensions equal to
`(length(knots[1]), length(knots[2]), ...)`.

`Ď• = GridDeformation(u::Array{T<:Real}, ...)` constructs the
deformation from a "plain" array. For a deformation in `N` dimensions,
`u` must have `N+1` dimensions, where the first dimension corresponds
to the displacement along each axis (and therefore `size(u,1) == N`).

Finally, `Ď• = GridDeformation((u1, u2, ...), ...)` allows you to
construct the deformation using an `N`-tuple of shift-arrays, each
with `N` dimensions.
"""
immutable GridDeformation{T,N,A<:AbstractArray,L} <: AbstractDeformation{T,N}
    u::A
    knots::NTuple{N,L}

    function (::Type{GridDeformation{T,N,A,L}}){T,N,A,L,FV<:SVector}(u::AbstractArray{FV,N}, knots::NTuple{N,L})
        typeof(u) == A || error("typeof(u) = $(typeof(u)), which is different from $A")
        length(FV) == N || throw(DimensionMismatch("Dimensionality $(length(FV)) must match $N knot vectors"))
        for d = 1:N
            size(u, d) == length(knots[d]) || error("size(u) = $(size(u)), but the knots specify a grid of size $(map(length, knots))")
        end
        new{T,N,A,L}(u, knots)
    end
    function (::Type{GridDeformation{T,N,A,L}}){T,N,A,L,FV<:SVector}(u::ScaledInterpolation{FV,N})
        new{T,N,A,L}(u, u.ranges)
    end
end

I’m a bit swamped right now so didn’t look at this in detail, but I hope this helps get you started.

Thanks @tim.holy! Hmph, I can’t quite connect the dots here… It would be phenomenal if someone from your lab, that uses this code (and maybe the rest of it), could chime in…?
Or if the whole thing is on github I can get a look at its context? Unregistered and all…

Thanks for the help!

You are correct it doesn’t seem right… the interpolation is only for rectangular grids , yours is not a rectangular grid.

The quickest way then would be to triangulate the input data and to interpolate on the triangulation.
It would be piecewise linear continues but not in general differentiable on every point.

see
https://github.com/JuliaMath/Interpolations.jl/issues/118

Right, thanks for the connection. I ping:ed them on that issue, maybe I can help push it along. But I suspect/hope @tim.holy might have some code that could be ported into that.

https://github.com/kbarbary/Dierckx.jl does irregular grids.

1 Like

Thanks @mauro3, you right about Dierckx.jl capabilities in irregular grids (which I think Interpolations.jl can also do), but I can’t see how Dierckx.jl can be used to interpolate from a 2D coordinate to a 2D coordinate, it’s just from 2D to a number describing a curve, not a pair of number describing coordinates.

it’s just from 2D to a number describing a curve
A surface you mean

Now you can have 2 interpolates were the first one will give you X and the second Y.

from the docs of Dierckx,jl

Spline2D(x, y, z; w=ones(length(x)), kx=3, ky=3, s=0.0)
Spline2D(x, y, z; kx=3, ky=3, s=0.0)

Fit a 2-d spline to the input data. x and y must be 1-d arrays.

If z is also a 1-d array, the inputs are assumed to represent unstructured data, with z[i] being the function value at point (x[i], y[i]). In this case, the lengths of all inputs must match.

Thanks for the discussion and references , it might come in handy in my work

Yay! I think you’re right. This seems to work.
Here’s a final MWE:

using Plots, Dierckx
image_x = [i + (rand()-.5)/4 for i in 1:10, j in 1:10]
image_y = [j + (rand()-.5)/4 for i in 1:10, j in 1:10]
plot(image_x,image_y, marker=:circle,color=:black)
plot!(image_x',image_y',color=:black,legend=false, grid=false)
x_in = [rand(1.5:9.5) + (rand()-.5)/2 for i in 1:10]
y_in = [rand(1.5:9.5) + (rand()-.5)/2 for i in 1:10]
scatter!(x_in,y_in,marker=:circle,color=:red)
png("a1.png")
world_x = linspace(125,300,10)
world_y = linspace(-5,53,10)
A = vec(world_x .+ 0world_y')
s = Spline2D(vec(image_x), vec(image_y), A, s = 100.)
x_out = evaluate(s, x_in, y_in)
A = vec(0world_x .+ world_y')
s = Spline2D(vec(image_x), vec(image_y), A, s = 100.)
y_out = evaluate(s, x_in, y_in)
x = repmat(world_x,1,10)
y = repmat(world_y',10)
plot(x,y, marker=:circle,color=:black)
plot!(x',y',color=:black,legend=false, grid=false)
scatter!(x_out, y_out, marker=:circle,color=:red)
png("a2.png")

and the two resulting images:
a1
a2

Thanks @TsurHerman!

6 Likes

Excuse me for being curious, but are you tracking moving animals with a drone?

1 Like

By all means! No, I’m tracking dung beetles from a stationary video camera (on a tripod). Since the animals are only moving in 2D and the ground is flat, seems like it should be easier to just use this grid thing (than the full on camera calibration setup).

1 Like