Ellipse fitting to determine rotation (Python)

Question

I am looking to fit an ellipse to some data points I have.

What I want: to determine the rotation angle of my data using an ellipse

The data: the data I have is in polar coordinates (θ, r)

theta = [0.0, 0.103, 0.206, 0.309, 0.412, 0.515, 0.618, 0.721, 0.824, 0.927, 1.03, 1.133, 1.236, 1.339, 1.442, 1.545, 1.648, 1.751, 1.854, 1.957, 2.06, 2.163, 2.266, 2.369, 2.472, 2.575, 2.678, 2.781, 2.884, 2.987, 3.09, 3.193, 3.296, 3.399, 3.502, 3.605, 3.708, 3.811, 3.914, 4.017, 4.12, 4.223, 4.326, 4.429, 4.532, 4.635, 4.738, 4.841, 4.944, 5.047, 5.15, 5.253, 5.356, 5.459, 5.562, 5.665, 5.768, 5.871, 5.974, 6.077, 6.18]

r = [84.48, 83.11, 77.67, 76.62, 90.12, 89.64, 84.07, 95.21, 104.63, 119.19, 125.19, 140.25, 146.33, 145.11, 164.0, 202.87, 214.81, 258.5, 281.94, 268.5, 224.76, 238.61, 270.08, 245.86, 220.04, 179.98, 181.51, 189.53, 172.87, 153.29, 138.32, 156.67, 146.21, 129.28, 139.76, 132.12, 138.73, 133.83, 136.15, 172.02, 163.2, 157.6, 142.73, 130.79, 130.24, 128.88, 124.7, 119.37, 115.28, 118.02, 117.89, 121.73, 115.13, 103.02, 84.43, 83.69, 82.26, 87.87, 88.84, 92.53, 94.67]

The algorithm so far:

0. define residuals and jacobian of residuals
1. use the scipy optimize.leastsq

(here is the walkthrough for those interested https://scipython.com/book/chapter-8-scipy/examples/non-linear-fitting-to-an-ellipse/ )

On my dataset however, the excentricity comes out negative, which it shoudln't be if it is an ellipse (0 < e < 1).

I've tried to add a rotation term that would depend on theta but without any luck so far. Here is the code to fit the ellipse without my extra term that messes everything up:

import numpy as np
from scipy import optimize
import pylab

def f(theta, p):
    a, e = p
    return a * (1 - e**2)/(1 - e*np.cos(theta))

def residuals(p, r, theta):
    """ Return the observed - calculated residuals using f(theta, p). """
    return r - f(theta, p)

def jac(p, r, theta):
    """ Calculate and return the Jacobian of residuals. """
    a, e = p
    da = (1 - e**2)/(1 - e*np.cos(theta))
    de = (-2*a*e*(1-e*np.cos(theta)) + a*(1-e**2)*np.cos(theta))/(1 -
                                                        e*np.cos(theta))**2
    return -da,  -de
    return np.array((-da, -de)).T

def fit_ellipse(theta, r, p0 = (1,0.5)):
    # Initial guesses for a, e
    p0 = (1, 0.5)
    plsq = optimize.leastsq(residuals, p0, Dfun=jac, args=(r, theta), col_deriv=True)
    #return plsq
    print(plsq)
    
    pylab.polar(theta, r, 'o')
    theta_grid = np.linspace(0, 2*np.pi, 200)
    pylab.polar(theta_grid, f(theta_grid, plsq[0]), lw=2)
    pylab.show()

fit_ellipse(theta, r, p0 = (1,0.5))

Answer 1

There is no need for non-linear regression (if no particular criteria of fitting is requiered). A simple linear regression leads to the result below :

The symbols and notations are consistent with Eqs. (15-23) from https://mathworld.wolfram.com/Ellipse.ht

In addition : Answer to a comment.

The equations used to draw the graph of the ellipse are :

Another manner to draw the ellipse (which avoids complex roots) is :

The parameter theta has to go from 0 to 2pi.