Index

Goodness of fit for ellipses

This page describes a goodness of fit measure for ellipses that is independent of the number of sample points that is scaled to the size of the ellipse. A sample Python program is included. The strategy for calculating the goodness of fit measure is to calculate the angle of each point relative to the center of the ellipse, and then find the intersection point of a line at that angle on the ellipse.

Standard measures of fitness, like the Mean Squared Error, do not provide meaningful comparisons for ellipses of different sizes. To be useful, an error measure has to be independent of the number of sample points and the size of the ellipse. For this reason, an error distance is calculated between the observed points and points on the ellipse. First, the angle of an observed point relative to the center of the ellipse is calculated. Next, we calculate the coordinates of the intersection of the ellipse with a line at the same angle as the point. Then, the distances for all the points are calculated, and an average error is obtained for all the points. This is the sum of all the error distances divided by the number of points.

The fitting error is defined as the percentage of the average error relative to the semiminor axis of the ellipse, which is half the width of the ellipse. Dividing by the semiminor axis basically scales the average error to a proportion of the ellipse to allow comparisons of the fitting errors for ellipses of different sizes.

$$\text{error_distance} = c = \sqrt{(y_1-y_2)^2 + (x_1-x_2)^2} $$
$$\text{average_error} = \frac{1}{n} \sum_{i=0}^n c_i $$
$$\text{fitting_error} = \frac {\text{average_error}\times 100}{\text{semiminor_axis}}$$

The following image was produced by the Python program listed below. The points correspond to a Carolina Bay called Herndon Bay in North Carolina at Lat. 34.862002, Lon. -78.944509. Well-preserved Carolina Bays have a mathematically elliptical geometry. A Python program for fitting ellipses to the Carolina Bays by the least squares method is available at GitHub.

# * * * * * * * * * * * * * * 
# February 11, 2025
# A. Zamora - Python program to illustrate the calculation of goodness of fit for an ellipse
#  using the average error distance of the points relative to the semiminor axis
# * * * * * * * * * * * * * * 
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
# * * * * * * *

# Observed points
x1 = [86.3421, 8.6616, 0.0, 44.4931, 467.6329, 579.7773, 682.4396, 819.9304, 911.5605, 955.6889, 949.3979, 914.5693, 856.0354, 775.5285, 671.3163, 559.1719, 435.9955, 325.4922, 202.4069]
y1 = [813.2214, 701.1104, 566.1105, 445.3329, 71.111, 20.5555, 0.0, 30.0, 110.5554, 213.222, 330.1108, 418.5551, 525.8884, 625.4438, 713.8882, 789.5548, 844.7769, 868.4436, 865.2214]

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()  # fit an ellipse to the set of points using a least squares method
ell.estimate(a_points)
xc, yc, a, b, phi = ell.params  

print(f' Number of points = {n1}')
print(f' center (x,y): {xc:.3f}, {yc:.3f}')  # center of ellipse
print(f' semimajor axis = {a:,.4f}, semiminor axis = {b:,.4f}')
print(f' angle of rotation phi = {phi:.4f} ({np.rad2deg(phi):.3f} deg.)')

ddt = 0  # >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
print(f' Average error distance relative to semiminor axis = {fitting_error:.4f}%')

fig = plt.figure(figsize=(6, 6))
axs = plt.subplot()
axs.axis('equal')
axs.scatter(x, y, color='red')           # observed points are red dots
axs.plot(calc_x, calc_y, 'x')            # calculated points on ellipse are marked with x
axs.scatter(xc, yc, color='red', s=100)  # center is a red dot
axs.set_title(f' Number of points: {n1}\n Average error distance relative to semiminor axis: {fitting_error:.3f}%')
ell_patch = Ellipse(xy=(xc, yc), width=2*a, height=2*b, angle=phi*180/np.pi, edgecolor='b', facecolor='none')
axs.add_patch(ell_patch)
plt.show()

---