Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Simple Julia code that can fit ellipses to 2D data.
"""
fit_ellipse(x, y) → a, b, θ, x_0, y_0, p
Fit an ellipse into the 2D data defined by `(x, y)`.
Return: semi-major axis length, semi-minor axis length, ellipse rotation,
center coordinates and a parameter container for quadratic form of ellipse,
which is just `M = hcat(x.^2, x.*y, y.^2, x, y); M*p = 1`.
Code modified from:
https://www.matecdev.com/posts/julia-least-squares-qr.html
using a lot of stuff from:
https://en.wikipedia.org/wiki/Ellipse#General_ellipse
"""
function fit_ellipse(x, y)
@assert length(x) == length(y)
# design matrix
M = hcat(x.^2, x.*y, y.^2, x, y)
p = M\ones(length(x)) # best fit parameters for the ellipse
A, B, C, D, E = p
F = -1.0
# do some clamping beause sometimes we calculate stuff like
# `sqrt(-7.559810648707566e-26)`
Δ = clamp(B^2 - 4*A*C, 0, Inf)
Λ = clamp((A-C)^2 + B^2, 0, Inf)
a, b = [-sqrt(
clamp(
2 * (A*E^2 + C*D^2 - B*D*E + (Δ)*F) *
( (A+C) + op(sqrt(Λ)) ), 0, Inf)) / (B^2 - 4A*C) for op in (+, -)
]
θ = atan((C - A - sqrt(Λ))/B)
x_0 = (2*C*D - B*E)/Δ
y_0 = (2*A*E - B*D)/Δ
return a, b, θ, x_0, y_0, p
end
function test_fit_ellipse(;θ = π/4, a = 3, b = 1.5, x_0 = 3, y_0 = -1,)
# Dummy data:
# Notice that parametric form of ellipse is just
# (x, y) = (a*cos(t), b*sin(t))
# so we just rotate and add center...
fx(t) = cos(θ)*a*cos(t) - sin(θ)*b*sin(t) + x_0
fy(t) = sin(θ)*a*cos(t) + cos(θ)*b*sin(t) + y_0
N = 200;
ts = LinRange(0,2π,N);
x = fx.(ts) + randn(N)*0.01;
y = fy.(ts) + randn(N)*0.01;
fig, ax = scatter(x, y; axis = (aspect = DataAspect(),))
fit_ellipse_and_plot!(ax, x, y; color = :red)
return fig
end
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment