fit a ellipse in Python given a set of points xi=(xi,yi)
10,152
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()
Author by
Gianni Spear
Updated on June 26, 2022Comments
-
Gianni Spear almost 2 years
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]thanks in advance for help
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) -
Gianni Spear over 11 years@untubu thanks. By the way, maybe i need to find a buld-in function to fit an ellipse fro a given 2D points. I am looking with google, but do you know some function?
-
unutbu over 11 years@Gianni: I don't know of a better algorithm, but I did define a
find_ellipsefunction which packages away all the(xmean, ymean)shifting, so the code looks nicer. -
Gianni Spear over 11 yearsin openCV there is cv.FitEllipse2: this is a short example (= with points i got a message error) PointArray2D32f = cv.CreateMat(1, len(points), cv.CV_32FC2)
-
Gianni Spear over 11 yearsfor (i, (x, y)) in enumerate(points): PointArray2D32f[0, i] = (x, y) # Fits ellipse to current contour. (center, size, angle) = cv.FitEllipse2(PointArray2D32f)
-
Gianni Spear over 11 yearsplease see the end of my first post (do you think can be correct?)
-
unutbu over 11 years@Gianni: Yes, I think the two methods are producing similar results, but expressed differently.
cv.FitEllipse2reports the same center asellipse_center. Whatcv.FitEllipse2reports assizeis twice whatellipse_axis_lengthreports. I think this is becausecv.FitEllipse2is reporting the length of the diameter, whileellipse_axis_lengthis reporting the radius. Finally,cv.FitEllipse2seems to be reporting the angle in degress, whileellipse_angle_of_rotationis given in radians:180-(0.708/pi*180)is roughly 139 degrees which is within spitting distance of 130 degrees.
