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?
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.
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 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.
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!
Excuse me for being curious, but are you tracking moving animals with a drone?
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).