CODEWITHSUNDEEP
juliadistance
leventov
leventov

Reputation: 15283

How to compute distance between two lines (sequence of 2-coordinate points) in Julia

I have two functions, y = f(x) and y = g(x), defined through Array{Float64,2} -- a sequence of (x, y) points.

I want to find such f'(x) = f(x+dx) + dy that the distance between the functions on the chart is minimal. For example, here, the distance is already quite small, but it could be even smaller.

enter image description here

In my particular case, the number of points in two arrays is identical, so I tried to minimise the following function (assuming p = [dx, dy]):

loss(p) = sum(abs2, mapslices(x -> x .- p, f .- g, dims=2))

But, apparently, this doesn't work well. Also, it would be nice to be able to work with arrays of different sizes.

Upvotes: 4

Views: 237

Answers (3)

Rafael_Guerra
Rafael_Guerra

Reputation: 130

@leventov, I have thought about this solution using Optim, which minimizes the vector along one curve with the minimum distances to the second. Arbitrary array lengths allowed and there is no need to interpolate/extrapolate.

using Optim, LinearAlgebra, Plots, Printf

function costf(x0)
    n, m = length(x1), length(x2)
    d = fill(Inf,n,1)
    for (ux, uy, i) in zip(x1, y1, 1:n)
        for (vx, vy) in zip(x2, y2)
            dij = norm([ux - vx - x0[1], uy - vy - x0[2]])
            dij < d[i] ? d[i] = dij : nothing
        end
    end
    return sum(d .^0.5)   # better metric for this problem
end

# x1, y1 = [-2.; -1; 0; 1; 2; 3], [4.; 1; 0; 1; 4; 9]
# x2, y2 = [0.; 0.5; 1; 2; 1], [0.; -0.75; -1; 0; 3]
x1 = -2:0.1:2;  y1 = x1 .^2
z = [exp(im*i/36*2pi) for i in 0:36]
x2, y2 = pi .+ 0.5*real(z), exp(1) .+ 0.5*(1.0 .+ imag(z))
x0 = [-2.0, -1.0]
res = optimize(costf, x0)
x0 = Optim.minimizer(res)

str = "Optim shift: dx= " * @sprintf("%.2f",x0[1]) * ", dy= " * @sprintf("%.2f",x0[2])
plot(x1,y1, color=:cyan, labels= "Input curve-1", legend=:topleft, aspectratio=1)
plot!(x2,y2, color=:blue, labels= "Input curve-2")
plot!(x2 .+ x0[1], y2 .+ x0[2], color=:red, labels= str)

Fig.1 Optim-shifted curve Fig.2 Optim-shifted circle

Upvotes: 3

Bogumił Kamiński
Bogumił Kamiński

Reputation: 69879

So the steps you would need to take is:

  1. find the range of dx that will ensure at least 50% overlap of the curves
  2. define both curves, call them f and g, using interpolate (without extrapolate) from Interpolations.jl using the points you have
  3. define a function that will calculate the aggregate distance between them; the easiest would be to use QuadGK.jl with the range of integration based on dx and the range of curves (so that we integrate without extrapolating, call the ends of this range xmin and xmax): quadgk(x -> (f(x)-g(x+dy)+dy)^2, xmin, xmax)); make it a function of dx and dy
  4. then use e.g. Optim.jl to find optimal dx and dy (with the constraint on dx)

Upvotes: 4

leventov
leventov

Reputation: 15283

One solution that I came up with (not sure how optimal is it):

using Interpolations, LsqFit

etp = extrapolate(interpolate((fx,), fy, Gridded(Linear())), Line())

model(x, p) = collect(map(i -> etp[i], x.+p[1])).+p[2]

fit = curve_fit(model, gx, gy, [0.0, 0.0])
res = coef(fit) # outputs an array like [dx, dy]

Also, this solution actually does extrapolate outside of the lines. I didn't find a way to change the fitting range depending on the parameters.

Upvotes: 1

Related Questions