3

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))
Improve this question
1
  • Have a look here. As stated by @jjacquelin , no need for non-linear regression. Commented May 17, 2021 at 7:09

2 Answers 2

Reset to default
1

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

Improve this answer
Sign up to request clarification or add additional context in comments.

Comments

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 :

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.

Improve this answer

7 Comments

Thank you @JJacquelin. Do you happen to have the code you used to produce the figure?
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.
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?
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, 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.
|

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.