Instantly share code, notes, and snippets.

# Datseris/fit_ellipse.jl

Created July 4, 2022 13:01
Show Gist options
• Save Datseris/4f1bed174e1bd2f7214bdb3d92885988 to your computer and use it in GitHub Desktop.
Simple Julia code that can fit ellipses to 2D data.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 """ 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