How can I score how well an ellipse fits my points?

Question

There are dozens of answers on here about how to fit an ellipse to your points, but none seem to suggest ways to quantify the fit. Here's some points and the ellipse fitted from this answer:

import numpy as np
from skimage.measure import EllipseModel
from matplotlib.patches import Ellipse
import matplotlib.pyplot as plt
x1 = [710.804, 710.117, 709.565, 709.036, 707.839, 707.424, 706.889, 705.913, 705.037, 704.58, 703.758,704.105, 704.552, 704.833, 705.204, 706.027, 706.932, 708.041, 708.849, 709.379, 709.797, 710.272,710.494, 710.871, 711.033, 711.018, 710.804]
y1 = [493.076, 491.902, 491.409, 490.947, 490.396, 491.456, 491.887, 492.917, 494.022, 494.882, 496.085,496.934, 497.723, 498.17, 498.656, 498.929, 499.248, 499.156, 498.768, 498.487, 497.853, 497.212,496.753, 495.957, 495.003, 493.997, 493.076]

points = list(zip(x1,y1))
print(points)
a_points = np.array(points)
x = a_points[:, 0]
y = a_points[:, 1]
ell = EllipseModel()
ell.estimate(a_points)
residuals = ell.residuals(a_points)
xc, yc, a, b, theta = ell.params

print("center = ",  (xc, yc))
print("angle of rotation = ",  theta)
print("axes = ", (a,b))
print('residuals= ', residuals)
fig, axs = plt.subplots(2, 1, sharex=True, sharey=True)
axs[0].scatter(x,y)
axs[1].scatter(x, y)
axs[1].scatter(xc, yc, color='red', s=100)
axs[1].set_xlim(x.min(), x.max())
axs[1].set_ylim(y.min(), y.max())
ell_patch = Ellipse((xc, yc), 2*a, 2*b, theta*180/np.pi, 
edgecolor='red', facecolor='none')
axs[1].add_patch(ell_patch)
plt.show()

center = (707.7248344321281, 495.2907644989805)

angle of rotation = 2.1051122302898926

axes = (4.783618808884847, 2.92938238506285)

residuals= [0.4031933 0.60237637 0.44816215 0.03196744 0.77430711 0.00493033 0.04060345 0.07892221 0.03482433 0.08585788 0.55088573 0.11629398 0.2794519 0.43588704 0.52172484 0.71671995 0.32973538 0.0389334 0.19938749 0.35342937 0.2419931 0.24903068 0.1998942 0.18446699 0.02005871 0.18610954 0.4031933 ]

Once I have created a fit, I want to calculate how well it fits the points. I can calculate the residuals using an EllipseModel function. I'm not sure how to use them for this purpose though. I can't use Pearson's correlation as it compares 1D lists. I thought I could calculate the closest point on the ellipse to each of my points by using the ratio of the residual to the distance from point to the ellipse centre of mass and transforming X and Y accordingly. I then thought to try dot product but that seems to require vectors, not points.

Can someone suggest a way to score the fit of the ellipse model?

Answer 1

I had a similar problem that required comparing the fitting errors of ellipses of different dimensions. The following code calculates the angle between each observed point and the center of the ellipse. Then, it calculates the distance between the observed point and a point on the ellipse at the same angle. The distances are added for all the points and normalized by dividing by the number of points. The average distance is then divided by the minor axis of the ellipse to calculate the error distance relative to the size of the ellipse.

import math
import numpy as np
from skimage.measure import EllipseModel
from matplotlib.patches import Ellipse
import matplotlib.pyplot as plt

# * * * * * * *
# Calculate the coordinates of the intersection point of the ellipse and the line
# * * * * * * *
def ellipse_line_intersection(Xc, Yc, L, W, phi, theta):
    # Get semimajor and semiminor axes
    a = L / 2
    b = W / 2
    
    # Calculate the parametric equations of the line
    dx = a * np.cos(theta)
    dy = a * np.sin(theta)
    
    # Rotate the line to align the ellipse with the coordinate axes
    cos_phi = np.cos(phi)
    sin_phi = np.sin(phi)
    x1_rot = (dx * cos_phi + dy * sin_phi)
    y1_rot = (-dx * sin_phi + dy * cos_phi)
    
    # Coefficients of the quadratic equation
    A = (x1_rot / a) ** 2 + (y1_rot / b) ** 2
    B = 2 * ((Xc * cos_phi + Yc * sin_phi - Xc * cos_phi - Yc * sin_phi) * x1_rot / a ** 2 + 
             (Xc * sin_phi - Yc * cos_phi - Xc * sin_phi + Yc * cos_phi) * y1_rot / b ** 2)
    C = ((Xc * cos_phi + Yc * sin_phi - Xc * cos_phi - Yc * sin_phi) ** 2 / a ** 2 +
         (Xc * sin_phi - Yc * cos_phi - Xc * sin_phi + Yc * cos_phi) ** 2 / b ** 2) - 1
    
    # Solve the quadratic equation
    discriminant = B ** 2 - 4 * A * C
    if discriminant < 0:
        return []  # No intersection
    
    t1 = (-B + np.sqrt(discriminant)) / (2 * A)
    t2 = (-B - np.sqrt(discriminant)) / (2 * A)
    
    # Calculate the intersection points
    x1 = Xc + dx * t1  # positive root
    y1 = Yc + dy * t1
    x2 = Xc + dx * t2  # negative root
    y2 = Yc + dy * t2
    
    return [(x1, y1), (x2, y2)]

# * * * * * * *
# Main program
# * * * * * * *

x1 = [710.804, 710.117, 709.565, 709.036, 707.839, 707.424, 706.889, 705.913, 705.037, 704.58, 703.758,704.105, 704.552, 704.833, 705.204, 706.027, 706.932, 708.041, 708.849, 709.379, 709.797, 710.272,710.494, 710.871, 711.033, 711.018, 710.804]
y1 = [493.076, 491.902, 491.409, 490.947, 490.396, 491.456, 491.887, 492.917, 494.022, 494.882, 496.085,496.934, 497.723, 498.17, 498.656, 498.929, 499.248, 499.156, 498.768, 498.487, 497.853, 497.212,496.753, 495.957, 495.003, 493.997, 493.076]

points = list(zip(x1,y1))
print(points)
a_points = np.array(points)
n1 = len(a_points)
calc_x = [0]*n1   # coordinates of points along the elliptical curve
calc_y = [0]*n1
x = a_points[:, 0]
y = a_points[:, 1]
ell = EllipseModel()
ell.estimate(a_points)
residuals = ell.residuals(a_points)
xc, yc, a, b, phi = ell.params  # a and b are semimajor and semiminor axes, respectively

print(f' n1={n1}')
ddt = 1  # >0 for debugging
fitting_error = 0

n = 0
distance = 0
sum_of_error_distances = 0
while n < n1:
  x_o = x[n]  # observed x
  y_o = y[n]  # observed y
  # Calculate theta (angle to the observed point) = arctan( (y-yc)/(x-xc) )
  theta = math.atan2( (y_o - yc), (x_o - xc) )
  if theta < 0 :
    theta = math.radians(360) + theta    # correct for quadrants III and IV
  if ddt > 0 :
    print(f'x_o[{n}]={x_o:,.4f} y_o[{n}]={y_o:,.4f} theta={theta:.4f} ({np.rad2deg(theta):.3f} deg.)')

  # calculate predicted intersection point on elliptical curve based on theta
  intersection_points = ellipse_line_intersection(xc, yc, 2*a, 2*b, phi, theta)
  calc_x[n], calc_y[n] = intersection_points[0]  # select only the positive solution
  if ddt > 0: 
    print(f'   calc_x[{n}]={calc_x[n]:.3f}, calc_y[{n}]={calc_y[n]:.3f}')
  # distance between observed and predicted points
  distance = math.sqrt((x_o - calc_x[n])**2 + (y_o - calc_y[n])**2)  
  sum_of_error_distances += distance

  n = n + 1

average_error_distance = sum_of_error_distances/n
print(f'Average error distance = {average_error_distance:.4f}')
fitting_error = average_error_distance*100/b
# Percent average error distance relative to minor axis 
print(f'Fitting error = {fitting_error:.4f}%')

print("center = ",  (xc, yc))
print("angle of rotation = ", phi)
print("axes = ", (a,b))
print('residuals= ', residuals)
fig, axs = plt.subplots(2, 1, sharex=True, sharey=True)
axs[0].scatter(x,y)
axs[1].scatter(x, y)
axs[1].plot(calc_x, calc_y, 'x')
axs[1].scatter(xc, yc, color='red', s=100)
axs[1].set_xlim(x.min(), x.max())
axs[1].set_ylim(y.min(), y.max())
axs[0].set_title(f' points: {n1}, fitting error: {fitting_error:.3f}%')
ell_patch = Ellipse(xy=(xc, yc), width=2*a, height=2*b, angle=phi*180/np.pi, edgecolor='red', facecolor='none')
axs[1].add_patch(ell_patch)
plt.show()

ellipse with Fitting error

Answer 2

A meaningful rating of the fit is given by one of

• the average distance of the points to the ellipse (can be computed using the residual function),

• the RMS distance, which turns out to coincide with the least-squares residual (normalized).

Depending on the goal of assessing the fit, you may consider an absolute value as above, or a relative one, for instance by dividing by the length of the major axis.

The average estimator is more robust than the RMS one.