# 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:

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:

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â€¦

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 `y`s?

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.

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:

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