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:

define residuals and jacobian of residuals
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))

Have a look [here](https://stackoverflow.com/a/61607790/803359). As stated by @jjacquelin , no need for non-linear regression. — mikuszefski, May 17 '21 at 07:09

score 1 · Answer 1 · answered Aug 26 '22 at 06:49

I cannot comment just yet, but I wanted to share the python implementation of the accepted answer. I have just wasted a few hours with this. The method does get the angle right, but it can lead to very weird fits. If anyone needs the code, it's in here along with two other ellipse fitting methods:

https://github.com/thelampire/flap_nstx/blob/unstable/tools/fit_objects.py

(I cross-checked probably three times, please feel free to criticize it or submit an issue/suggestion).

JJacquelin · Accepted Answer · 2021-05-25T18:27:50.863

0

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.

edited May 25 '21 at 18:27

answered May 14 '21 at 22:19

JJacquelin

1,529
1
9
11

Thank you @JJacquelin. Do you happen to have the code you used to produce the figure? – Sorade May 20 '21 at 17:59
There is no code in fact. This is a so simple numerical calculus that Mathcad is sufficient. The above figure shows the sceen copy. There is nothing more. The figure also is drawn with the graphing tool implemented in Mathcad. – JJacquelin May 20 '21 at 21:19
I'm sorry, but this is not so simple for me. Did you convert the coordinates from polar to cartesian first? What is k, and n? – Sorade May 25 '21 at 09:34
All the calculus are in Cartesian coordinates, not in polar. n is the number of points. The points are numbered from k=1 (point number 1, cooedinates x_1,y_1) to k=n (point number n , coordinates x_n,y_n). – JJacquelin May 25 '21 at 09:46
JJacquelin, the equation for plotting the graphs appears to work in most cases but in some instances the number inside the square-root is negative and as such throws an error. Is that normal? For example with an x of 84.48. – Sorade May 25 '21 at 10:23
Yes, of course. In ploting the ellipse point after point with successives values of x, if x is outside the range of definition of the ellipse the value of y is not real. Before to compute the root one have to check if the argument is positif. If not the root must not be computed (there is no point to draw corrresponding to this value of x). – JJacquelin May 25 '21 at 11:27
Another manner to draw the ellipse is added at the end of the answer. – JJacquelin May 25 '21 at 18:29

Ellipse fitting to determine rotation (Python)

2 Answers2