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]

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)

1 Answer: 

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