我试图使用matplotlib在不同的二维点集上拟合一个椭圆。我使用了来自于 此处，这是代码。
from matplotlib.patches import Ellipse
import matplotlib.pyplot as plt
import numpy
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_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
但是，当我用matplotlib绘制椭圆时，有时椭圆拟合得很好，有时需要将椭圆旋转90度才能拟合（蓝色椭圆是旋转90度，红色椭圆是没有额外旋转的）这是我的代码。
def plot_ellipse(x, y):
a = fitEllipse(x, y)
center = ellipse_center(a)
phi = ellipse_angle_of_rotation2(a)
axes = ellipse_axis_length(a)
a, b = axes
ell = Ellipse(center, 2*a, 2*b, phi*180 / np.pi, facecolor=(1,0,0,0.2), edgecolor=(0,0,0,0.5))
ell_rotated = Ellipse(center, 2*a, 2*b, phi*180 / np.pi + 90, facecolor=(0,0,1,0.2), edgecolor=(0,0,0,0.5))
fig, ax = plt.subplots()
ax.add_patch(ell)
ax.add_patch(ell_rotated)
plt.scatter(x, y)
plt.show()
x1 = np.array([238, 238, 238, 237, 237, 237, 237, 237, 236, 236, 236, 236, 237, 238,
239, 240, 240, 241, 242, 243, 243, 244, 245, 246, 247, 248, 249, 250,
251, 252, 253, 254, 255, 255, 256, 257, 258, 259, 260, 261, 262, 263,
264, 265, 266, 266, 267, 267, 268, 268, 269, 269, 270, 270, 271, 271,
271, 271, 271, 272, 272, 272, 272, 272, 273, 273, 273, 273, 273, 273,
274, 274, 274, 274, 274, 274, 275, 275, 275, 275, 275, 275, 275, 275,
275, 275, 274, 274, 274, 274, 274, 274, 274, 274, 274, 273, 273, 273,
272, 272, 272, 272, 271, 271, 271, 270, 270, 269, 268, 268, 267, 266,
266, 265, 265, 264, 263, 262, 261, 260, 259, 258, 257, 256, 256, 255,
254, 253, 252, 251, 250, 249, 248, 247, 246, 245, 245, 244, 243, 242,
241, 240, 239, 238, 237, 236, 235, 235, 235, 234, 234, 233, 233, 232,
232, 232, 231, 231, 230, 230, 229, 229, 229, 229, 229, 229, 229, 229,
228, 228, 228, 228, 228, 228, 228, 228, 228, 228, 228, 228, 228, 228,
228, 228, 228, 228, 228, 228, 228, 228, 228, 228, 228, 228, 228, 228,
228, 229, 229, 229, 229, 229, 229, 229, 229, 229, 230, 230, 230, 230,
231, 231, 231, 232, 232, 232, 232, 233, 233, 233, 234, 235, 236, 237,
238, 239, 240, 241, 242])
y1 = np.array([283, 284, 285, 286, 287, 288, 289, 290, 291, 292, 293, 293, 294, 295,
296, 297, 298, 299, 300, 301, 301, 301, 301, 301, 301, 301, 301, 301,
301, 301, 301, 301, 301, 301, 301, 300, 300, 299, 299, 298, 298, 297,
297, 296, 296, 296, 295, 294, 293, 292, 291, 290, 289, 288, 287, 286,
286, 285, 284, 283, 282, 281, 280, 279, 278, 277, 276, 276, 275, 274,
273, 272, 271, 270, 269, 268, 267, 266, 265, 264, 264, 263, 262, 261,
260, 259, 258, 257, 256, 255, 254, 253, 252, 252, 251, 250, 249, 248,
247, 246, 245, 244, 243, 242, 242, 241, 240, 239, 238, 237, 236, 235,
234, 233, 233, 232, 232, 231, 230, 230, 229, 228, 228, 227, 227, 227,
227, 226, 226, 226, 226, 226, 226, 225, 225, 225, 225, 226, 226, 227,
228, 229, 229, 230, 231, 231, 232, 232, 233, 234, 235, 236, 237, 238,
239, 240, 241, 242, 243, 244, 245, 246, 246, 247, 248, 249, 250, 251,
252, 253, 254, 255, 256, 257, 257, 258, 259, 260, 261, 262, 263, 264,
265, 266, 267, 268, 269, 270, 271, 272, 273, 273, 274, 275, 276, 277,
278, 279, 280, 281, 282, 283, 284, 285, 285, 286, 287, 288, 289, 290,
291, 292, 293, 294, 295, 296, 297, 298, 299, 299, 300, 301, 301, 302,
303, 304, 304, 305, 306])
x2 = np.array( [235, 236, 237, 238, 239, 240, 241, 242, 243, 243, 244, 245, 246, 247,
248, 249, 250, 251, 252, 253, 254, 255, 255, 256, 257, 258, 259, 260,
261, 262, 263, 264, 265, 266, 266, 267, 268, 269, 270, 270, 271, 272,
273, 274, 274, 274, 275, 275, 276, 276, 277, 277, 278, 278, 279, 279,
279, 279, 280, 280, 280, 280, 281, 281, 281, 281, 282, 282, 282, 282,
282, 282, 282, 282, 282, 281, 281, 281, 281, 281, 281, 281, 281, 281,
280, 280, 280, 280, 280, 279, 279, 279, 279, 278, 278, 277, 277, 276,
276, 275, 275, 274, 274, 274, 273, 272, 271, 270, 269, 268, 267, 266,
265, 264, 263, 263, 262, 261, 260, 259, 258, 257, 256, 255, 254, 253,
252, 252, 251, 250, 249, 248, 247, 246, 245, 244, 243, 242, 241, 240,
240, 239, 238, 237, 236, 235, 234, 233, 232, 232, 231, 231, 230, 230,
229, 229, 228, 228, 227, 227, 227, 227, 227, 227, 227, 227, 227, 227,
226, 226, 226, 226, 226, 226, 226, 226, 226, 226, 226, 226, 226, 227,
227, 227, 227, 227, 227, 227, 227, 228, 228, 228, 228, 228, 228, 229,
229, 230, 230, 231, 231, 232, 232, 233, 233, 234, 234, 235, 235, 235,
236, 237, 238, 239, 240, 241, 242, 243, 244])
y2 = np.array(
[279, 280, 281, 282, 283, 284, 285, 286, 287, 287, 287, 288, 288, 289,
289, 290, 290, 290, 291, 291, 292, 292, 292, 292, 291, 291, 290, 290,
289, 289, 288, 288, 287, 287, 287, 286, 285, 284, 283, 282, 281, 280,
279, 278, 278, 277, 276, 275, 274, 273, 272, 271, 270, 269, 268, 267,
267, 266, 265, 264, 263, 262, 261, 260, 259, 258, 257, 256, 255, 255,
254, 253, 252, 251, 250, 249, 248, 247, 246, 245, 244, 244, 243, 242,
241, 240, 239, 238, 237, 236, 235, 234, 234, 233, 232, 231, 230, 229,
228, 227, 226, 225, 224, 224, 223, 222, 222, 221, 220, 219, 218, 217,
217, 216, 215, 215, 215, 215, 215, 215, 215, 214, 214, 214, 214, 214,
214, 214, 214, 215, 215, 216, 216, 217, 217, 217, 218, 218, 219, 219,
219, 220, 221, 222, 223, 223, 224, 225, 226, 226, 227, 228, 229, 230,
231, 232, 233, 234, 235, 236, 236, 237, 238, 239, 240, 241, 242, 243,
244, 245, 246, 247, 248, 249, 250, 251, 251, 252, 253, 254, 255, 256,
257, 258, 259, 260, 261, 262, 263, 264, 265, 266, 267, 268, 268, 269,
270, 271, 272, 273, 274, 275, 276, 277, 278, 279, 280, 281, 282, 282,
283, 284, 284, 285, 286, 287, 287, 288, 289])
plot_ellipse(x1, y1)
plot_ellipse(x2, y2)
而这里是绘图的截图。
如你所见，不旋转的（红色）椭圆能很好地拟合x1,y1数据，但旋转后的椭圆（蓝色）能拟合x2,y2数据。
我很困惑，如果我错过了什么，什么时候我需要将椭圆旋转90º，什么时候不需要？
这是我的猜测，没有100%检查数学。从定义和要解决的拉格朗日来看，一切看起来都很好，也很符合逻辑。我认为向量 a 是正确的，这样的内部参数 a 到 f. 但是，代码中提到，要最小化的函数与缩放因子无关。所以在计算轴角时，可能会遇到一个符号问题。我的猜测是，这个问题发生在 arctan() 函数。一个解决办法是在计算的 a 向量，即 a -> -a 如果 a[-1] < 0.
更好的，可能也是工作（我只是测试了2个案例从OP）是取代
if a > c:
return np.arctan( 2 * b / ( a - c ) ) / 2
else:
return np.pi / 2 + np.arctan( 2 * b / (a - c ) ) / 2
由
if a > c:
return np.arctan2( 2 * b, ( a - c ) ) / 2
else:
return np.pi / 2 + np.arctan2( 2 * b, ( a - c) ) / 2
的论点，则可知。arctan 是来自一个师。arctan2 是要去的。
补充一点。我认为有多个 return 一个函数中的语句应该省略。在这种短函数中，还是很容易处理的，但我认为这不是好的做法。
更新
在联系原配合码的作者时，他提到了 arctan2 是在他提供的版本中使用的。GitHub
这段代码看起来更简洁，我建议使用这段代码来代替主页上的片段。
其他想法
其实，我觉得可以用更一致、更容易操作的方式，椭圆的参数是如何从 a 所以我写下了这段代码
# -*- coding: utf-8 -*-
import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse
import numpy as np
RAD = 180. / np.pi
DEGREE = 1. / RAD
def rot( a ):
"""
simple rotation matrix in 2D
"""
return np.array(
[ [ +np.cos( a ), -np.sin( a ) ],
[ +np.sin( a ), +np.cos( a ) ] ]
)
def fit_ellipse( x, y ):
"""
main fit from the original publication:
http://nicky.vanforeest.com/misc/fitEllipse/fitEllipse.html
"""
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 ] = +2
C[ 2, 0 ] = +2
C[ 1, 1 ] = -1
E, V = np.linalg.eig( np.dot( np.linalg.inv( S ), C ) )
n = np.argmax( np.abs( E ) )
a = V[ :, n ]
return a
def ell_parameters( a ):
"""
New function substituting the original 3 functions for
axis, centre and angle.
We start by noting that the linear term is due to an offset.
Getting rid of it is equivalent to find the offset.
Starting with the Eq.
xT A x + bT x + c = 0 and transforming x -> x - t
we get a new linear term. By demanding that this term vanishes
we get the Eq.
b = (AT + A ) t.
Hence, an easy way to write down how to get t
"""
A = np.array( [ [ a[0], a[1]/2. ], [ a[1]/2., a[2] ] ] )
b = np.array( [ a[3], a[4] ] )
t = np.dot( np.linalg.inv( np.transpose( A ) + A ), b )
"""
the transformation changes the constant term, which we need
for proper scaling
"""
c = a[5]
cnew = c - np.dot( t, b ) + np.dot( t, np.dot( A, t ) )
Anew = A / (-cnew)
# ~cnew = cnew / (-cnew) ### debug only
"""
now it is in the form xT A x - 1 = 0
and we know that A is a rotation of the matrix
( 1 / a² 0 )
B = ( )
( 0 1 / b² )
where a and b are the semi axes of the ellipse
it is hence A = ST B S
We note that rotation does not change the eigenvalues, which are
the diagonal elements of matrix B. Moreover, we note that
the matrix of eigenvectors rotates B into A
"""
E, V = np.linalg.eig( Anew )
"""
so we have
B = VT A V
and consequently
A = V B VT
where V is of a form as given by the function rot() from above
"""
# ~B = np.dot( np.transpose(V), np.dot( Anew, V ) ) ### debug only
phi = np.arccos( V[ 0, 0 ] )
"""
checking the sin for changes in sign to detect angles above 180°
"""
if V[ 0, 1 ] < 0:
phi = 2 * np.pi - phi
### cw vs ccw and periodicity of pi
phi = -phi % np.pi
return np.sqrt( 1. / E ), phi * RAD, -t
"""
That's it. One might put some additional work/thought in the 180°
and cw vs ccw thing, as it is a bit messy.
"""
"""
creating some test data
"""
xl = np.linspace(-3,2.5, 10)
yl = np.fromiter( (2.0 * np.sqrt( 1 - ( x / 3. )**2 ) for x in xl ), np.float )
xl = np.append(xl,-xl)
yl = np.append(yl,-yl)
R = rot( -103.01 * DEGREE ) ### check different angles
# ~R = rot( 153 * DEGREE ) results in singular matrix !!!...strange
xyrot = np.array( [ np.dot(R, [ x, y ] )for x, y in zip( xl, yl ) ] )
xl = xyrot[:,0] + 7
yl = xyrot[:,1] + 16.4
"""
fitting
"""
avec = fit_ellipse( xl, yl )
(a, b), phi, t = ell_parameters( avec )
ell = Ellipse(
t, 2 * a, 2 * b, phi,
facecolor=( 1, 0, 0, 0.2 ), edgecolor=( 0, 0, 0, 0.5 )
)
"""
plotting
"""
fig = plt.figure()
ax = fig.add_subplot( 1, 1, 1 )
ax.add_patch( ell )
ax.scatter( xl ,yl )
plt.show()
我想这段代码应该也不会有90°的问题。如果去掉所有注释，它是相当紧凑的，（从数学的角度来看）可读性很强。
注
我在 inv( S )中遇到了一个问题。这个矩阵变成了奇数。一个变通的方法是：将所有数据旋转一个小角度，然后将计算出的 t 回。同时从计算的角度中减去角度 phi.