CODEWITHSUNDEEP
pythonellipsedata-fitting
Sorade
Sorade

Reputation: 937

Ellipse fitting to determine rotation (Python)

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:

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

Upvotes: 3

Views: 1159

Answers (2)

thelampire
thelampire

Reputation: 31

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).

Upvotes: 1

JJacquelin
JJacquelin

Reputation: 1705

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 :

enter image description here

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 :

enter image description here

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

enter image description here

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

Upvotes: 0

Related Questions