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

Posted by Gianni on Stack Overflow See other posts from Stack Overflow or by Gianni
Published on 2012-11-29T21:59:20Z Indexed on 2012/11/29 23:04 UTC
Read the original article Hit count: 341

Filed under:
|
|

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

© Stack Overflow or respective owner

Related posts about python

Related posts about geometry