Установите Ellipse поверх 2D-данных, используя Python и matplotlib

Question

Я пытаюсь наложить эллипс на мои разные наборы 2D точек, используя 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

Как видите, не повернутый (красный) Эллипс хорошо подходит для данных x1, y1, но повернутый эллипс (синий) подходит для данных x2, y2.

Я запутался, если что-то здесь упустил, когда мне нужно повернуть эллипс на 90º а когда мне не нужно?

Answer 1

Вот мое предположение, без проверки математики на 100%. Глядя на определение и лагранжиан, который будет решен, все выглядит хорошо и логично c. Я думаю, что вектор 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 - это то, что go для.

Еще одно замечание: я думаю, что код с несколькими операторами 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.