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

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]

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 enter image description here

plot of 2D point and center of fit ellipse enter image description here

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)

Would you please post a figure to show how the result looks like? — Tengis
– Tengis, Commented Nov 29, 2012 at 22:05
@Eric: the 2 graph is the 2D Points and the center of the ellipse [ 560415.53298363+0.j 6368878.84576771+0.j]. The function adds +0.j at x and y — Gianni Spear
– Gianni Spear, Commented Nov 29, 2012 at 22:48

unutbu · Accepted Answer · 2012-11-29 23:13:44Z

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

enter image description here

@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?
@Gianni: I don't know of a better algorithm, but I did define a find_ellipse function which packages away all the (xmean, ymean) shifting, so the code looks nicer.
in 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)
for (i, (x, y)) in enumerate(points): PointArray2D32f[0, i] = (x, y) # Fits ellipse to current contour. (center, size, angle) = cv.FitEllipse2(PointArray2D32f)
@Gianni: Yes, I think the two methods are producing similar results, but expressed differently. cv.FitEllipse2 reports the same center as ellipse_center. What cv.FitEllipse2 reports as size is twice what ellipse_axis_length reports. I think this is because cv.FitEllipse2 is reporting the length of the diameter, while ellipse_axis_length is reporting the radius. Finally, cv.FitEllipse2 seems to be reporting the angle in degress, while ellipse_angle_of_rotation is given in radians: 180-(0.708/pi*180) is roughly 139 degrees which is within spitting distance of 130 degrees.

Collectives™ on Stack Overflow

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

1 Answer 1

6 Comments

Linked

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

6 Comments

Linked

Related