如何从 SymPy 中提取椭圆的中心坐标、高度、宽度和 phi 以绘制拟合椭圆?
how to extract center coordinates, height, width and phi of an ellipse from SymPy to plot a fitted ellipse?
我一直在使用 lsq-ellipse
包,在其中我使用以下代码获取椭圆的坐标:
from ellipse import LsqEllipse
from matplotlib.patches import Ellipse
coords_D0 = np.array(coords_D0)
reg = LsqEllipse().fit(coords_D0)
center_D0, width_D0, height_D0, phi_D0 = reg.as_parameters()
print(f'center: {center_D0[0]:.3f}, {center_D0[1]:.3f}')
print(f'width: {width_D0:.3f}')
print(f'height: {height_D0:.3f}')
print(f'phi: {phi_D0:.3f}')
但是,我的 coords_D0
变量由三个坐标组成,这导致了以下错误:
ValueError: Received too few samplesGot 3 features, 5 or more required.
但是,在查看了一些包和在线之后,我发现 sympy
也可以做 Ellipse,我知道你可以提取 centre, vradius and hradius
来自 sympy
。但是,我想知道如何从 sympy 获取 width、height 和 phi,它是否与 lsq-ellipse
[=43= 相同? ] 要在matplotlib 的Ellipse 中使用的包?我使用 matplotlib 中 lsq-ellipse
包中的值来形成椭圆部分,它可以在以下代码行中找到:
代码:
ellipse_D0 = Ellipse(xy=center_D0, width=2*width_D0, height=2*height_D0, angle=np.rad2deg(phi_D0),edgecolor='b', fc='None', lw=2, label='Fit', zorder=2)
我的坐标如下:
coords_D0 =
-1.98976 -1.91574
-0.0157721 2.5438
2.00553 -0.628061
# another points
coords_D1 =
-0.195518 0.0273673
-0.655686 -1.45848
-0.447061 -0.168108
# another points
coords_D2 =
-2.28529 0.91896
-2.43207 0.446211
-2.23044 0.200087
附带问题:
有没有办法将椭圆拟合到这些坐标(一般情况下,3 个坐标或更多)?
假设 OP 是关于最小体积封闭椭圆,我建议以下解决方案。
#! /usr/bin/python3
# coding=utf-8
import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse
import numpy as np
from mymodules3.mvee import mvee
coords= list()
coords.append( np.array([
-1.98976, -1.91574,
-0.0157721, 2.5438,
2.00553, -0.628061
]).reshape(3,-1) )
coords.append( np.array([
-0.195518, 0.0273673,
-0.655686, -1.4584,8
-0.447061, -0.168108,
]).reshape(3,-1)
)
coords.append( np.array([
-2.28529, 0.91896,
-2.43207, 0.446211,
-2.23044, 0.200087
]).reshape(3,-1)
)
fig = plt.figure()
ax = fig.add_subplot( 1, 1, 1 )
for i, col in enumerate( ['k', 'g', 'm'] ):
sol = mvee( coords[i] )
e =Ellipse(
sol[0],
width=2 * sol[1][0],
height=2 * sol[1][1],
angle=sol[2] * 180/3.1415926,
color=col, alpha=0.5,
zorder = -1000 + i
)
ax.scatter( coords[i][:,0], coords[i][:,1], c=col, zorder=10 * i )
ax.add_artist( e )
plt.show()
提供
mvee
是基于类似问题的 SE answer。
"""
NMI : 2021-11-11
Minimum Volume Enclosing Ellipsoids, see e.g.
NIMA MOSHTAGH : MINIMUM VOLUME ENCLOSING ELLIPSOIDS
or
Linus Källberg : Minimum_Enclosing_Balls_and_Ellipsoids (Thesis)
"""
from warnings import warn
from numpy import pi
from numpy import sqrt
from numpy import arccos
from numpy import dot, outer
from numpy import diag, transpose
from numpy import append
from numpy import asarray
from numpy import ones
from numpy import argmax
from numpy.linalg import inv
from numpy.linalg import norm
from numpy.linalg import eig
def mvee( data, tolerance=1e-4, maxcnt=1000 ):
"""
param data: list of xy data points
param tolerance: termination condition for iterative approximation
param maxcnt: maximum number of iterations
type data: iterable of float
type tolerance: float
return: (offset, semiaxis, angle)
return type: ( (float, float), (float, float), float )
"""
locdata = asarray( data )
N = len( locdata )
if not locdata.shape == ( N, 2):
raise ValueError ( " data must be of shape( n, 2 )" )
if tolerance >= 1 or tolerance <= 0:
raise ValueError (" 0 < tolerance < 1 required")
if not isinstance( maxcnt, int ):
raise TypeError
if not maxcnt > 0:
raise ValueError
count = 1
err = 1
d = 2
d1 = d + 1
u = ones( N ) / N
P = transpose( locdata )
Q = append( P, ones( N ) ).reshape( 3, -1 )
while ( err > tolerance):
X = dot( Q, dot( diag( u ), transpose( Q ) ) )
M = diag(
dot(
transpose( Q ),
dot(
inv( X ),
Q
)
)
)
maximum = max( M )
j = argmax( M )
step_size = ( maximum - d1 ) / ( d1 * ( maximum - 1 ) )
new_u = ( 1 - step_size ) * u
new_u[ j ] += step_size
err = norm( new_u - u )
count = count + 1
u = new_u
if count > maxcnt:
warn(
"Process did not converge in {} steps".format(
count - 1
),
UserWarning
)
break
U = diag( u )
c = dot( P, u )
A = inv(
dot(
P,
dot( U, transpose( P ) )
) - outer( c, c )
) / d
E, V = eig( A )
phiopt = arccos( V[ 0, 0 ] )
if V[ 0, 1 ] < 0:
phiopt = 2 * pi - phiopt
### cw vs ccw and periodicity of pi
phiopt = -phiopt % pi
sol = ( c, sqrt( 1.0 / E ), phiopt)
return sol
我一直在使用 lsq-ellipse
包,在其中我使用以下代码获取椭圆的坐标:
from ellipse import LsqEllipse
from matplotlib.patches import Ellipse
coords_D0 = np.array(coords_D0)
reg = LsqEllipse().fit(coords_D0)
center_D0, width_D0, height_D0, phi_D0 = reg.as_parameters()
print(f'center: {center_D0[0]:.3f}, {center_D0[1]:.3f}')
print(f'width: {width_D0:.3f}')
print(f'height: {height_D0:.3f}')
print(f'phi: {phi_D0:.3f}')
但是,我的 coords_D0
变量由三个坐标组成,这导致了以下错误:
ValueError: Received too few samplesGot 3 features, 5 or more required.
但是,在查看了一些包和在线之后,我发现 sympy
也可以做 Ellipse,我知道你可以提取 centre, vradius and hradius
来自 sympy
。但是,我想知道如何从 sympy 获取 width、height 和 phi,它是否与 lsq-ellipse
[=43= 相同? ] 要在matplotlib 的Ellipse 中使用的包?我使用 matplotlib 中 lsq-ellipse
包中的值来形成椭圆部分,它可以在以下代码行中找到:
代码:
ellipse_D0 = Ellipse(xy=center_D0, width=2*width_D0, height=2*height_D0, angle=np.rad2deg(phi_D0),edgecolor='b', fc='None', lw=2, label='Fit', zorder=2)
我的坐标如下:
coords_D0 =
-1.98976 -1.91574
-0.0157721 2.5438
2.00553 -0.628061
# another points
coords_D1 =
-0.195518 0.0273673
-0.655686 -1.45848
-0.447061 -0.168108
# another points
coords_D2 =
-2.28529 0.91896
-2.43207 0.446211
-2.23044 0.200087
附带问题:
有没有办法将椭圆拟合到这些坐标(一般情况下,3 个坐标或更多)?
假设 OP 是关于最小体积封闭椭圆,我建议以下解决方案。
#! /usr/bin/python3
# coding=utf-8
import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse
import numpy as np
from mymodules3.mvee import mvee
coords= list()
coords.append( np.array([
-1.98976, -1.91574,
-0.0157721, 2.5438,
2.00553, -0.628061
]).reshape(3,-1) )
coords.append( np.array([
-0.195518, 0.0273673,
-0.655686, -1.4584,8
-0.447061, -0.168108,
]).reshape(3,-1)
)
coords.append( np.array([
-2.28529, 0.91896,
-2.43207, 0.446211,
-2.23044, 0.200087
]).reshape(3,-1)
)
fig = plt.figure()
ax = fig.add_subplot( 1, 1, 1 )
for i, col in enumerate( ['k', 'g', 'm'] ):
sol = mvee( coords[i] )
e =Ellipse(
sol[0],
width=2 * sol[1][0],
height=2 * sol[1][1],
angle=sol[2] * 180/3.1415926,
color=col, alpha=0.5,
zorder = -1000 + i
)
ax.scatter( coords[i][:,0], coords[i][:,1], c=col, zorder=10 * i )
ax.add_artist( e )
plt.show()
提供
mvee
是基于类似问题的 SE answer。
"""
NMI : 2021-11-11
Minimum Volume Enclosing Ellipsoids, see e.g.
NIMA MOSHTAGH : MINIMUM VOLUME ENCLOSING ELLIPSOIDS
or
Linus Källberg : Minimum_Enclosing_Balls_and_Ellipsoids (Thesis)
"""
from warnings import warn
from numpy import pi
from numpy import sqrt
from numpy import arccos
from numpy import dot, outer
from numpy import diag, transpose
from numpy import append
from numpy import asarray
from numpy import ones
from numpy import argmax
from numpy.linalg import inv
from numpy.linalg import norm
from numpy.linalg import eig
def mvee( data, tolerance=1e-4, maxcnt=1000 ):
"""
param data: list of xy data points
param tolerance: termination condition for iterative approximation
param maxcnt: maximum number of iterations
type data: iterable of float
type tolerance: float
return: (offset, semiaxis, angle)
return type: ( (float, float), (float, float), float )
"""
locdata = asarray( data )
N = len( locdata )
if not locdata.shape == ( N, 2):
raise ValueError ( " data must be of shape( n, 2 )" )
if tolerance >= 1 or tolerance <= 0:
raise ValueError (" 0 < tolerance < 1 required")
if not isinstance( maxcnt, int ):
raise TypeError
if not maxcnt > 0:
raise ValueError
count = 1
err = 1
d = 2
d1 = d + 1
u = ones( N ) / N
P = transpose( locdata )
Q = append( P, ones( N ) ).reshape( 3, -1 )
while ( err > tolerance):
X = dot( Q, dot( diag( u ), transpose( Q ) ) )
M = diag(
dot(
transpose( Q ),
dot(
inv( X ),
Q
)
)
)
maximum = max( M )
j = argmax( M )
step_size = ( maximum - d1 ) / ( d1 * ( maximum - 1 ) )
new_u = ( 1 - step_size ) * u
new_u[ j ] += step_size
err = norm( new_u - u )
count = count + 1
u = new_u
if count > maxcnt:
warn(
"Process did not converge in {} steps".format(
count - 1
),
UserWarning
)
break
U = diag( u )
c = dot( P, u )
A = inv(
dot(
P,
dot( U, transpose( P ) )
) - outer( c, c )
) / d
E, V = eig( A )
phiopt = arccos( V[ 0, 0 ] )
if V[ 0, 1 ] < 0:
phiopt = 2 * pi - phiopt
### cw vs ccw and periodicity of pi
phiopt = -phiopt % pi
sol = ( c, sqrt( 1.0 / E ), phiopt)
return sol