Least-Squares Fitting of an Ellipse with General Ellipse Equation (including tilt and translation)

Question

I would like to utilize the following code to fit some 2D data to an ellipse, which I got from this post.

import numpy as np
import matplotlib.pyplot as plt
alpha = 5
beta = 3
N = 500
DIM = 2

np.random.seed(2)

# Generate random points on the unit circle by sampling uniform angles
theta = np.random.uniform(0, 2*np.pi, (N,1))
eps_noise = 0.2 * np.random.normal(size=[N,1])
circle = np.hstack([np.cos(theta), np.sin(theta)])

# Stretch and rotate circle to an ellipse with random linear tranformation
B = np.random.randint(-3, 3, (DIM, DIM))
noisy_ellipse = circle.dot(B) + eps_noise

# Extract x coords and y coords of the ellipse as column vectors
X = noisy_ellipse[:,0:1]
Y = noisy_ellipse[:,1:]

# Formulate and solve the least squares problem ||Mx - b ||^2
M = np.hstack([X**2, X * Y, Y**2, X, Y])
b = np.ones_like(X)
x = np.linalg.lstsq(M, b)[0].squeeze()

# Print the equation of the ellipse in standard form
print('The ellipse is given by {0:.3}x^2 + {1:.3}xy+{2:.3}y^2+{3:.3}x+{4:.3}y = 1'.format(x[0], x[1],x[2],x[3],x[4]))

# Plot the noisy data
plt.scatter(X, Y, label='Data Points')

# Plot the original ellipse from which the data was generated
phi = np.linspace(0, 2*np.pi, 1000).reshape((1000,1))
c = np.hstack([np.cos(phi), np.sin(phi)])
ground_truth_ellipse = c.dot(B)
plt.plot(ground_truth_ellipse[:,0], ground_truth_ellipse[:,1], 'k--', label='Generating Ellipse')

# Plot the least squares ellipse
x_coord = np.linspace(-5,5,300)
y_coord = np.linspace(-5,5,300)
X_coord, Y_coord = np.meshgrid(x_coord, y_coord)
Z_coord = x[0] * X_coord ** 2 + x[1] * X_coord * Y_coord + x[2] * Y_coord**2 + x[3] * X_coord + x[4] * Y_coord
plt.contour(X_coord, Y_coord, Z_coord, levels=[1], colors=('r'), linewidths=2)

plt.legend()
plt.xlabel('X')
plt.ylabel('Y')
plt.show()

However, I need to use the full general equation: Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 1, which I got from here and includes tilt and translation of the ellipse.

In the above code, the M matrix:

M = np.hstack([X**2, X * Y, Y**2, X, Y])

assumes I'm fitting Ax^2 + Bxy + Cy^2 + Dx + Ey = 1.

Question: How do I "add" in the F coefficient, or is that not possible using this method?

Answer 1

Indeed it is possible with a little knowledge in optimization. We know that the trivial solution for F is certainly F=1 and other parameters equal to 0, but we don't want that and there is probably another solution other than the trivial solution. So let's get into some theory.

In the standard least squares method we have

Standard least squares

But there are several ways to extend this method so that it can meet the demands of each problem. In your case, what we know is that the trivial solution resets the parameters A, B, C, D and E, that is, the simplest solution is to think that these for meters cannot be Zero

So, a solution is to add a regularizer that adds to our objective function, a term that adds an error smaller and smaller if our parameters distance from zero. One of the ways is: to obtain the inverse of the sum of the parameters we have (the higher the sum of the parameters, the smaller the penalty from the regularizer will be.), in the way that Regularizer = alpha/sum(|x|). And sum(|x|) is L1 norm in Linear Algebra. Increase alpha to increase the strength of the regularizer

So our LSTSQ is

Our modified least squares

but be careful, because the norm of x cannot be zero, we will deal with this in code.

So, let's code

Step 1

Import minimize from scipy

from scipy.optimize import minimize

Step 2

Create your own minimize function

(the adjusted variable always comes before the other parameters, that's why x comes before)

def target_func(x, M, b):
    norm = np.linalg.norm(x, ord=1)     # get L1 norm of x
    norm = norm if norm != 0 else 1e-6  # Verify Zero, in case Zero, place a tolerance (1e-6 in this example)
    reg = 100/(norm)                    # Evaluate regularizer (using alpha= 100 for this exemple)
    return np.linalg.norm(M@x - b, ord=2) + reg # Write de minimize function based on lstsq

Step 3

Create parameters and apply on minimize function

M = np.hstack([X**2, X * Y, Y**2, X, Y, np.ones_like(X)]) # Add ones for param F
b = np.ones_like(X).flatten()                             # Create b (as flatten for prevent problems)
x0 = np.random.rand(6)                                    # Create Initial solution for our problem

res_log = minimize(target_func,                           # Our minimize target
                   x0,                                    # Initial solution
                   args=(M, b),                           # Args for problem (M and b)
)
x = res_log.x # Get the x 

Step 4

Execute the same code for plot, now with the x[5] (F param).

# Print the equation of the ellipse in standard form
print('The ellipse is given by {0:.3}x^2 + {1:.3}xy + {2:.3}y^2 + {3:.3}x + {4:.3}y + {5:.3} = 1'.format(*x))

# Plot the noisy data
plt.scatter(X, Y, label='Data Points')

# Plot the original ellipse from which the data was generated
phi = np.linspace(0, 2*np.pi, 1000).reshape((1000,1))
c = np.hstack([np.cos(phi), np.sin(phi)])
ground_truth_ellipse = c.dot(B)
plt.plot(ground_truth_ellipse[:,0], ground_truth_ellipse[:,1], 'k--', label='Generating Ellipse')

# Plot the least squares ellipse
x_coord = np.linspace(-5,5,300)
y_coord = np.linspace(-5,5,300)
X_coord, Y_coord = np.meshgrid(x_coord, y_coord)
Z_coord = x[0] * X_coord ** 2 + x[1] * X_coord * Y_coord + x[2] * Y_coord**2 + x[3] * X_coord + x[4] * Y_coord + x[5] # Add x[5]
plt.contour(X_coord, Y_coord, Z_coord, levels=[1], colors=('r'), linewidths=2)

plt.legend()
plt.xlabel('X')
plt.ylabel('Y')
plt.show()

Result

The ellipse is given by -0.337x^2 + -0.143xy + -0.543y^2 + -0.0325x + -0.0331y + 5.33 = 1

Minimize of Scipy

Basically the "minimize" function works by iterations, where it calculates where the cost function has the lowest value, so it adjusts the parameters at each iteration until it finds the best ones given a tolerance and maximum number of iterations allowed.

See more on docs: https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html