如何从 SymPy 中提取椭圆的中心坐标、高度、宽度和 phi 以绘制拟合椭圆？

Question

我一直在使用 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 也可以做椭圆，我知道你可以从 centre, vradius and hradius 中提取em>sympy。但是，我想知道如何从 sympy 中获取 width、height 和 phi，它是否与 Ellipse 中使用的 lsq-ellipse 包相同matplotlib？我使用 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 个坐标或更多）？

Answer 1

假设 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