A Python script to regress an ellipse from a list of (x, y) coordinate pairs by solving systems of linear equations and using monte carlo

ellipse_regress.py

# This is a Python script to regress an ellipse from a list of (x, y) coordinate pairs
# it heuristically optimises the coefficients of the following equation:
# Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0
# It tries to minimise the MSE between observed points and points predicted by the model

# Method:
# It may not be the most optimal method but it's reasonably fast, was quick to implement and works
# - Randomly picks 5 points out of the list of points
# - Solves for the coefficients by solving the corresponding system of equations
# - The solution is evaluated by randomly selecting coordinate pairs and evaluating their MSE
# - This is iterated many times and the solution with the lowest MSE is taken to be an approximation for the solution

from sympy import *
import random

F = 1

# X = [789, 532, 195, ...]
# Y = [452, 115, 415, ...]

def get_random_sol():
	eqns = []
	A, B, C, D, E = symbols("A B C D E")

	for _ in range(5):
		i = random.randint(0, len(X)-1)
		x = X[i]
		y = Y[i]
		eqns.append(A*x**2 + B*x*y + C*y**2 + D*x + E*y + F)

	coeffs = tuple(linsolve(eqns, [A, B, C, D, E]))[0]

	if E in coeffs[-1].free_symbols:
		# coeffs[-1] is symbolic
		# => wasn't able to solve (e.g because of repeated point(s))
		return None
	return coeffs

def eval_dist(x, y_actual, coeffs):
	A, B, C, D, E = coeffs
	y = Symbol('y')
	y_preds = map(lambda a: a.evalf(), solve(A*x**2 + B*x*y + C*y**2 + D*x + E*y + F, y))
	dists = [((y_actual - y_pred)**2).evalf() for y_pred in y_preds if im(y_pred) == 0]

	if len(dists) > 0:
		return min(dists)

def eval_sol(coeffs, iters=10):
	n = 0
	s = 0  # sum of square deviations
	for _ in range(iters):
		i = random.randint(0, len(X)-1)
		dist = eval_dist(X[i], Y[i], coeffs)

		if dist is not None:
			s += dist
			n += 1

	if n > 0:
		return s / n  # = mse

def get_sol(iters=50):
	min_mse = float('inf')
	min_coeffs = None

	for _ in range(iters):
		coeffs = get_random_sol()
		if coeffs is None: continue

		mse = eval_sol(coeffs)
		if mse is None: continue

		if mse < min_mse:
			min_mse = mse
			min_coeffs = coeffs

	print("MSE: ", min_mse)
	return min_coeffs + (F,)

print(get_sol())