Simple Julia code that can fit ellipses to 2D data.

fit_ellipse.jl

"""
    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