39.3: LSF Example - Estimating the best Ellipses

%matplotlib inline
import matplotlib.pylab as plt
import numpy as np
import sympy as sym
sym.init_printing(use_unicode=True)

Now consider the following sets of points. Think of these as observations of planet moving in an elliptical orbit.

x=[0, 1.0, 1.1, -1.1, -1.2, 1.3]
y =[2*1.5, 2*1.0, 2*-0.99, 2*-1.02, 2*1.2, 2*0]

plt.scatter(x,y)
plt.axis('equal')

In this problem we want to try to fit an ellipse to the above data. First lets look at a general equation for an ellipse:

\[ \frac{(u+x)^2}{a^2} + \frac{(v+y)^2}{b^2} = 1 \]

Where \(u\) and \(v\) are the \(x\) and \(y\) coordinates for the center of the ellipse and \(a\) and \(b\) are the lengths of the axes sizes of the ellipse. A quick search on how to plot an ellipse in python comes up with the following example:

# Code from: https://stackoverflow.com/questions/10952060/plot-ellipse-with-matplotlib-pyplot-python

u=1.     #x-position of the center
v=0.5    #y-position of the center
a=2.     #radius on the x-axis
b=1.5    #radius on the y-axis

t = np.linspace(0, 2*np.pi, 100)
plt.plot( u+a*np.cos(t) , v+b*np.sin(t) )
plt.grid(color='lightgray',linestyle='--')
plt.show()

Notice this example uses equations of the form:

\[t = [0, \dots, 2\pi] \nonumber \]

\[x = u+a\cos(t) \nonumber \]

\[y = v+b\sin(t) \nonumber \]

Turns out that this form of the equation is easier to plot and the variables \(u\),\(v\),\(a\),\(b\) are the same as our original equation.

Now lets expand the original equation (equation \(\PageIndex{1}\) from above) and we get the following:

\[x^2−2ux-u^2+y^2−2vy+v^2=r^2 \]

Question

Why can’t we convert equation \(\PageIndex{2}\) into the form \(Ax=b\) and solve using Least Squares Fit? Discuss with your group and be prepared to share your thought with the class.

If we look at our data more closely we can simplify equations \(\PageIndex{1}\) and \(\PageIndex{2}\) by assuming the the centroid \((u,v)\) is at the origin. This assumption results in the following equation:

\[ \frac{x^2}{a^2} + \frac{y^2}{b^2} = 1 \]

Notice we can rewrite this into a standard linear set of equations by defining \(c_o = \frac{1}{a^2}\) and \(c_1 = \frac{1}{b^2}\) and rewriting the equation as follows:

\[ c_0x^2 + c_1y^2 = 1 \]

Do This

Given that we know the \(x\) and \(y\) values of our point observations, equation \(\PageIndex{4}\) is now linear and can be solved using Least Squares Fit. Using the observation points from above construct A and b as numpy matrixes for the overdefined system \(Ax=b\)

sym.Matrix(A)

sym.Matrix(b)

Do This

Solve the above over defined system of linear equations for \(c_0\) and \(c_1\) using LSF.

# Put your answer to the above question here

Assuming we have \(c\) in the correct format, we can now calculate \(a\) and \(b\) from the solution for \(c_0\) and \(c_1\) calculated in the previous step and plot using our plotting code:

c = 1/np.sqrt(np.abs(c))
b=c[1,0] 
a=c[0,0]
print(a,b)

u=0     #x-position of the center
v=0     #y-position of the center

t = np.linspace(0, 2*np.pi, 100)
plt.plot(u+a*np.cos(t) , v+b*np.sin(t) )
plt.scatter(x,y)
plt.grid(color='lightgray',linestyle='--')
plt.axis('equal');