Я пытаюсь подогнать эллипс к набору точек, описываемых координатами x и y.
Я нашел подробное объяснение, как это сделать, здесь http://nicky.vanforeest.com/misc/fitEllipse/fitEllipse.html и опробовал код, но, похоже, он не работает.
Он правильно находит центр, но угол и оси совершенно неправильны, как вы можете видеть на этом изображении: https://i.imgur.com/VLEeNKQ.png
Красные точки - мои точки данных, а синяя точка -эллипс взят из полученных параметров.Теперь, данные не являются идеальным эллипсом, но они не подходят.Я хотел бы приблизить соответствие к фактическим данным.
Вот код, о котором идет речь.
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])
def ellipse_angle_of_rotation2( a ):
b,c,d,f,g,a = a[1]/2, a[2], a[3]/2, a[4]/2, a[5], a[0]
if b == 0:
if a > c:
return 0
else:
return np.pi/2
else:
if a > c:
return np.arctan(2*b/(a-c))/2
else:
return np.pi/2 + np.arctan(2*b/(a-c))/2
А вот мой полный набор данных.У кого-нибудь есть идеи, почему он не устанавливается правильно?
# --------------------------------------------------------------------------
x = np.array([ 5727.53135, 7147.62235, 10330.93573, 8711.17228, 7630.40262,
4777.24983, 4828.27655, 9449.94416, 5203.81323, 6299.44811,
6494.21906])
y = np.array([67157.77567 , 66568.50068 , 55922.56257 , 54887.47348 ,
65150.14064 , 66529.91705 , 65934.25548 , 55351.57612 ,
63123.5103 , 67181.141725, 56321.36025 ])
# -----------------------------------------------------------------------------
a = fitEllipse(x,y)
center = ellipse_center(a)
#phi = ellipse_angle_of_rotation(a)
phi = ellipse_angle_of_rotation2(a)
axes = ellipse_axis_length(a)
# get the individual axes
a, b = axes
from matplotlib.patches import Ellipse
import matplotlib.pyplot as plt
ell = Ellipse(center, a, b, phi)
fig, ax = plt.subplots(subplot_kw={'aspect': 'equal'})
ax.add_artist(ell)
ell.set_clip_box(ax.bbox)
ax.set_xlim(0, 100000)
ax.set_ylim(0, 100000)
plt.show()
scat = plt.scatter(x, y, c = "r")