如何从 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