﻿ geometry - fit a ellipse in Python given a set of points xi=(xi,yi) - DeveloperSite- developersite.org

# geometry - fit a ellipse in Python given a set of points xi=(xi,yi)

Keywords：python

Question:

I am computing a series of index from a 2D points (x,y). One index is the ratio between minor and major axis. To fit the ellipse i am using the following post

when i run these function the final results looks strange because the center and the axis length are not in scale with the 2D points

``````center =  [  560415.53298363+0.j  6368878.84576771+0.j]
angle of rotation =  (-0.0528033467597-5.55111512313e-17j)
axes =  [0.00000000-557.21553487j  6817.76933256  +0.j]
``````

``````import numpy as np
from numpy.linalg import eig, inv

def fitEllipse(x,y):
x = x[:,np.newaxis]
y = y[:,np.newaxis]
D =  np.hstack((x*x, x*y, y*y, x, y, np.ones_like(x)))
S = np.dot(D.T,D)
C = np.zeros([6,6])
C[0,2] = C[2,0] = 2; C[1,1] = -1
E, V =  eig(np.dot(inv(S), C))
n = np.argmax(np.abs(E))
a = V[:,n]
return a

def ellipse_center(a):
b,c,d,f,g,a = a[1]/2, a[2], a[3]/2, a[4]/2, a[5], a[0]
num = b*b-a*c
x0=(c*d-b*f)/num
y0=(a*f-b*d)/num
return np.array([x0,y0])

def ellipse_angle_of_rotation( a ):
b,c,d,f,g,a = a[1]/2, a[2], a[3]/2, a[4]/2, a[5], a[0]
return 0.5*np.arctan(2*b/(a-c))

def ellipse_axis_length( a ):
b,c,d,f,g,a = a[1]/2, a[2], a[3]/2, a[4]/2, a[5], a[0]
up = 2*(a*f*f+c*d*d+g*b*b-2*b*d*f-a*c*g)
down1=(b*b-a*c)*( (c-a)*np.sqrt(1+4*b*b/((a-c)*(a-c)))-(c+a))
down2=(b*b-a*c)*( (a-c)*np.sqrt(1+4*b*b/((a-c)*(a-c)))-(c+a))
res1=np.sqrt(up/down1)
res2=np.sqrt(up/down2)
return np.array([res1, res2])

if __name__ == '__main__':

points = [(560036.4495758876, 6362071.890493258),
(560036.4495758876, 6362070.890493258),
(560036.9495758876, 6362070.890493258),
(560036.9495758876, 6362070.390493258),
(560037.4495758876, 6362070.390493258),
(560037.4495758876, 6362064.890493258),
(560036.4495758876, 6362064.890493258),
(560036.4495758876, 6362063.390493258),
(560035.4495758876, 6362063.390493258),
(560035.4495758876, 6362062.390493258),
(560034.9495758876, 6362062.390493258),
(560034.9495758876, 6362061.390493258),
(560032.9495758876, 6362061.390493258),
(560032.9495758876, 6362061.890493258),
(560030.4495758876, 6362061.890493258),
(560030.4495758876, 6362061.390493258),
(560029.9495758876, 6362061.390493258),
(560029.9495758876, 6362060.390493258),
(560029.4495758876, 6362060.390493258),
(560029.4495758876, 6362059.890493258),
(560028.9495758876, 6362059.890493258),
(560028.9495758876, 6362059.390493258),
(560028.4495758876, 6362059.390493258),
(560028.4495758876, 6362058.890493258),
(560027.4495758876, 6362058.890493258),
(560027.4495758876, 6362058.390493258),
(560026.9495758876, 6362058.390493258),
(560026.9495758876, 6362057.890493258),
(560025.4495758876, 6362057.890493258),
(560025.4495758876, 6362057.390493258),
(560023.4495758876, 6362057.390493258),
(560023.4495758876, 6362060.390493258),
(560023.9495758876, 6362060.390493258),
(560023.9495758876, 6362061.890493258),
(560024.4495758876, 6362061.890493258),
(560024.4495758876, 6362063.390493258),
(560024.9495758876, 6362063.390493258),
(560024.9495758876, 6362064.390493258),
(560025.4495758876, 6362064.390493258),
(560025.4495758876, 6362065.390493258),
(560025.9495758876, 6362065.390493258),
(560025.9495758876, 6362065.890493258),
(560026.4495758876, 6362065.890493258),
(560026.4495758876, 6362066.890493258),
(560026.9495758876, 6362066.890493258),
(560026.9495758876, 6362068.390493258),
(560027.4495758876, 6362068.390493258),
(560027.4495758876, 6362068.890493258),
(560027.9495758876, 6362068.890493258),
(560027.9495758876, 6362069.390493258),
(560028.4495758876, 6362069.390493258),
(560028.4495758876, 6362069.890493258),
(560033.4495758876, 6362069.890493258),
(560033.4495758876, 6362070.390493258),
(560033.9495758876, 6362070.390493258),
(560033.9495758876, 6362070.890493258),
(560034.4495758876, 6362070.890493258),
(560034.4495758876, 6362071.390493258),
(560034.9495758876, 6362071.390493258),
(560034.9495758876, 6362071.890493258),
(560036.4495758876, 6362071.890493258)]

a_points = np.array(points)
x = a_points[:, 0]
y = a_points[:, 1]
from pylab import *
plot(x,y)
show()
a = fitEllipse(x,y)
center = ellipse_center(a)
phi = ellipse_angle_of_rotation(a)
axes = ellipse_axis_length(a)

print "center = ",  center
print "angle of rotation = ",  phi
print "axes = ", axes

from pylab import *
plot(x,y)
plot(center[0:1],center[1:], color = 'red')
show()
``````

each vertex is a xi,y,i point

plot of 2D point and center of fit ellipse

using OpenCV i have the following result:

``````import cv

PointArray2D32f = cv.CreateMat(1, len(points), cv.CV_32FC2)
for (i, (x, y)) in enumerate(points):
PointArray2D32f[0, i] = (x, y)
# Fits ellipse to current contour.
(center, size, angle) = cv.FitEllipse2(PointArray2D32f)

(center, size, angle)
((560030.625, 6362066.5),(10.480490684509277, 17.20206642150879),144.34889221191406)
``````

The calculation for `fitEllipse` is returning bogus results because the values for `x` and `y` are very large compared to the variation between the values. If you try printing out the eigenvalues `E` for example, you see

``````array([  0.00000000e+00 +0.00000000e+00j,
0.00000000e+00 +0.00000000e+00j,
0.00000000e+00 +0.00000000e+00j,
-1.36159790e-12 +8.15049878e-12j,
-1.36159790e-12 -8.15049878e-12j,   1.18685632e-11 +0.00000000e+00j])
``````

They are all practically zero! Clearly there is some kind of numerical inaccuracy here.

You can fix the problem by moving the mean of your data closer to zero so that the values are more "normal"-sized and the variation between the numbers becomes more significant.

``````x = a_points[:, 0]
y = a_points[:, 1]
xmean = x.mean()
ymean = y.mean()
x = x-xmean
y = y-ymean
``````

You can then successfully find the center, phi and axes, and then re-shift the center back to (xmean, ymean):

``````center = ellipse_center(a)
center[0] += xmean
center[1] += ymean
``````

``````import numpy as np
import numpy.linalg as linalg
import matplotlib.pyplot as plt

def fitEllipse(x,y):
x = x[:,np.newaxis]
y = y[:,np.newaxis]
D =  np.hstack((x*x, x*y, y*y, x, y, np.ones_like(x)))
S = np.dot(D.T,D)
C = np.zeros([6,6])
C[0,2] = C[2,0] = 2; C[1,1] = -1
E, V =  linalg.eig(np.dot(linalg.inv(S), C))
n = np.argmax(np.abs(E))
a = V[:,n]
return a

def ellipse_center(a):
b,c,d,f,g,a = a[1]/2, a[2], a[3]/2, a[4]/2, a[5], a[0]
num = b*b-a*c
x0=(c*d-b*f)/num
y0=(a*f-b*d)/num
return np.array([x0,y0])

def ellipse_angle_of_rotation( a ):
b,c,d,f,g,a = a[1]/2, a[2], a[3]/2, a[4]/2, a[5], a[0]
return 0.5*np.arctan(2*b/(a-c))

def ellipse_axis_length( a ):
b,c,d,f,g,a = a[1]/2, a[2], a[3]/2, a[4]/2, a[5], a[0]
up = 2*(a*f*f+c*d*d+g*b*b-2*b*d*f-a*c*g)
down1=(b*b-a*c)*( (c-a)*np.sqrt(1+4*b*b/((a-c)*(a-c)))-(c+a))
down2=(b*b-a*c)*( (a-c)*np.sqrt(1+4*b*b/((a-c)*(a-c)))-(c+a))
res1=np.sqrt(up/down1)
res2=np.sqrt(up/down2)
return np.array([res1, res2])

def find_ellipse(x, y):
xmean = x.mean()
ymean = y.mean()
x -= xmean
y -= ymean
a = fitEllipse(x,y)
center = ellipse_center(a)
center[0] += xmean
center[1] += ymean
phi = ellipse_angle_of_rotation(a)
axes = ellipse_axis_length(a)
x += xmean
y += ymean
return center, phi, axes

if __name__ == '__main__':

points = [(560036.4495758876, 6362071.890493258),
(560036.4495758876, 6362070.890493258),
(560036.9495758876, 6362070.890493258),
(560036.9495758876, 6362070.390493258),
(560037.4495758876, 6362070.390493258),
(560037.4495758876, 6362064.890493258),
(560036.4495758876, 6362064.890493258),
(560036.4495758876, 6362063.390493258),
(560035.4495758876, 6362063.390493258),
(560035.4495758876, 6362062.390493258),
(560034.9495758876, 6362062.390493258),
(560034.9495758876, 6362061.390493258),
(560032.9495758876, 6362061.390493258),
(560032.9495758876, 6362061.890493258),
(560030.4495758876, 6362061.890493258),
(560030.4495758876, 6362061.390493258),
(560029.9495758876, 6362061.390493258),
(560029.9495758876, 6362060.390493258),
(560029.4495758876, 6362060.390493258),
(560029.4495758876, 6362059.890493258),
(560028.9495758876, 6362059.890493258),
(560028.9495758876, 6362059.390493258),
(560028.4495758876, 6362059.390493258),
(560028.4495758876, 6362058.890493258),
(560027.4495758876, 6362058.890493258),
(560027.4495758876, 6362058.390493258),
(560026.9495758876, 6362058.390493258),
(560026.9495758876, 6362057.890493258),
(560025.4495758876, 6362057.890493258),
(560025.4495758876, 6362057.390493258),
(560023.4495758876, 6362057.390493258),
(560023.4495758876, 6362060.390493258),
(560023.9495758876, 6362060.390493258),
(560023.9495758876, 6362061.890493258),
(560024.4495758876, 6362061.890493258),
(560024.4495758876, 6362063.390493258),
(560024.9495758876, 6362063.390493258),
(560024.9495758876, 6362064.390493258),
(560025.4495758876, 6362064.390493258),
(560025.4495758876, 6362065.390493258),
(560025.9495758876, 6362065.390493258),
(560025.9495758876, 6362065.890493258),
(560026.4495758876, 6362065.890493258),
(560026.4495758876, 6362066.890493258),
(560026.9495758876, 6362066.890493258),
(560026.9495758876, 6362068.390493258),
(560027.4495758876, 6362068.390493258),
(560027.4495758876, 6362068.890493258),
(560027.9495758876, 6362068.890493258),
(560027.9495758876, 6362069.390493258),
(560028.4495758876, 6362069.390493258),
(560028.4495758876, 6362069.890493258),
(560033.4495758876, 6362069.890493258),
(560033.4495758876, 6362070.390493258),
(560033.9495758876, 6362070.390493258),
(560033.9495758876, 6362070.890493258),
(560034.4495758876, 6362070.890493258),
(560034.4495758876, 6362071.390493258),
(560034.9495758876, 6362071.390493258),
(560034.9495758876, 6362071.890493258),
(560036.4495758876, 6362071.890493258)]

fig, axs = plt.subplots(2, 1, sharex = True, sharey = True)
a_points = np.array(points)
x = a_points[:, 0]
y = a_points[:, 1]
axs[0].plot(x,y)
center, phi, axes = find_ellipse(x, y)
print "center = ",  center
print "angle of rotation = ",  phi
print "axes = ", axes

axs[1].plot(x, y)
axs[1].scatter(center[0],center[1], color = 'red', s = 100)
axs[1].set_xlim(x.min(), x.max())
axs[1].set_ylim(y.min(), y.max())

plt.show()
``````