在Python中给定一组点xi =（xi，yi） [英] fit a ellipse in Python given a set of points xi=(xi,yi)
本文介绍了在Python中给定一组点xi =（xi，yi）的处理方法，对大家解决问题具有一定的参考价值，需要的朋友们下面随着小编来一起学习吧！
问题描述
当我运行这些函数时，最终结果看起来很奇怪，因为中心和轴线长度与2D点不成比例
center = [560415.53298363 + 0.j 6368878.84576771 + 0.j]
旋转角度=（-0.0528033467597-5.55111512313e-17j）
轴= [ 0.00000000-557.21553487j 6817.76933256 + 0.j]
提前致谢帮助
从numpy.linalg导入numpy作为np
导入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（DT，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]
返回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 * ba * c
x0 =（c * db * f）/ num
y0 =（a * fb * d）/ num
返回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 /（ac））
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 * （b + b * b * 2 * b * d * fa * c * g） /（（ac）*（ac））） - （c + a））
down2 =（b * ba * c）*（（ac）* np.sqrt（1 + 4 * b * b / （ac）*（ac））） - （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]从pylab导入的
*
plot（x， y）
show（）
a = fitEllipse（x，y）
center = ellipse_center（a）
phi = ellipse_angle_of_rotation（a）
axes = ellipse_axis_length（a）
printcenter =，center
printrotation of angle =，phi
打印axes =，axes
from pylab import *
plot（x，y）
plot（center [0：1]，center [1：]，color ='red'）
show（）
每个顶点都是xi，y，i点
二维点和拟合中心椭圆图
我得到以下结果：
import cv
PointArray2D32f = cv.CreateMat（1，len（po （$）
for（i，（x，y））枚举（点）：
PointArray2D32f [0，i] =（x，y）
＃适合椭圆到当前轮廓。
（center，size，angle）= cv.FitEllipse2（PointArray2D32f）
（center，size，angle）
（（560030.625,6362066.5），（10.480490684509277,17.20206642150879）， 144.34889221191406）
解决方案 x 和
y 的值非常大，所以fitEllipse 返回假冒结果值之间的差异。例如，如果您尝试打印出特征值 E ，您会看到 数组（[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]）
它们几乎都是零！显然这里有一些数字不准确。
您可以通过将数据的平均值移近零来解决问题，以便数值更正常化，并且数字之间的差异变得更加显着。
x = a_points [:, 0]
y = a_points [:, 1]
xmean = x.mean（）
ymean = y.mean（）
x = x-xmean
y = y-ymean
然后，您可以成功找到中心，phi和坐标轴，然后将中心重新移回（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
$ b $ 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 DT，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]
返回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 * ba * c
x0 = c * db * f）/ num
y0 =（a * fb * d）/ num
返回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] 0]
return 0.5 * np.arctan（2 * b /（ac））
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 * fa * c * g）
down1 =（b * ba * c）*（（ca）* np.sqrt（1 + 4 * b * b /（（ac）*（ac））） - （c + a））
down2 =（b * ba * c）*（（ac）* np.sqrt（1 + 4 * b （上/下2）$ b（（ac）*（ac））） - （c + a））
res1 = $ b return np.array（[res1，res2]）
def find_ellipse（x，y）：
xmean = x.mean（）
ymean = y.mean（ ）
x - = xm ean
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
返回中心，phi，轴
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.390
（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.390
（560023.9495758876,6362063.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， （560027.4495758876,6362069.890493258），
（560027.9495758876,6362068.890493258），
（560027.9495758876,6362069.390493258），
（560028.4495758876,6362069.390493258），
（560028.4495758876,6362069.890493258），
（560033.4495758876， 6362069.890
（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）
printcenter =，center
printrotation of rotation =，phi
printaxes =，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（）
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)
解决方案
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()
这篇关于在Python中给定一组点xi =（xi，yi）的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
查看全文
