0

I have been working with lsq-ellipse package where I get the coordinates of ellipse with the following code below:

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}') 

However, my coords_D0 variable consists of three coordinates which caused the following error:

ValueError: Received too few samplesGot 3 features, 5 or more required. 

But, after looking into some packages and online, I found that sympy also can do Ellipse and I understand that you can extract the centre, vradius and hradius from sympy. But, I would like to know how to get the width, height and phi from sympy and will it be the same as the lsq-ellipse package to be used in Ellipse of matplotlib? I use the values from lsq-ellipse package in matplotlib to form the ellipse part and it can be found in the following code line:

Code:

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) 

My coordinates are the following:

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 

Side Question:

Is there a way to fit an ellipse to these coordinates (in general, 3 coordinates or more)?

16
  • You use fit. Two crossing ellipses can have 4 points in common, so there is no unique solution giving less than 5. Did you read the docs? Commented Nov 9, 2021 at 7:53
  • @mikuszefski which package are you hinting? Do you mean fit in the Lsq package I get the error for requiring more features more than 5. Commented Nov 10, 2021 at 8:15
  • Hi, so if I got it right , the LsqEllipse().fit() method takes a list of (x, y) to fit an ellipse. To give a proper unique result, you need 5 or more values. So naturally you get an error. Commented Nov 10, 2021 at 16:19
  • @mikuszefski I understand. I have been looking online and found these links (stackoverflow.com/questions/39693869/…) (stackoverflow.com/questions/47873759/…) (stackoverflow.com/questions/52818206/…) (stackoverflow.com/questions/39693869/…), however, I can’t get proper centre point, width, height and phi using my coordinates. Commented Nov 11, 2021 at 18:58
  • So we are clear that you cannot fit an ellipse providing 3 points. Which coordinates are we talking about then? I think I put an easy to follow code here Commented Nov 12, 2021 at 6:21

1 Answer 1

2

Assuming that the OP is about the Minimum Volume Enclosing Ellipse, I'd suggest the following solution.

#! /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() 

providing

MVEE Solution

The mvee is is based on an SE answer on a similar question.

""" 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 
Sign up to request clarification or add additional context in comments.

3 Comments

Thank you for the answer. I looked into the code, link and I managed to get it to work. I just have a question, about the tolerance and maxcnt. What parameters I could use to fit the ellipse to the data points?
@WDpad159 As you can see, my image is using standard settings. maxcnt shouldn't be the problem as you expect convergence and the loop should terminate before this is reached. I did not go into fully understanding the theory so I do not have a good feeling for what this tolerance means in details. Sure, it is the change of length in the u vector. which somehow encodes the eigenvalues of A. As this is a gradient method, it is testing how much the gradient is changing. If you have a good assumption about the shape of the extremum that gives you a hint how far you are away...
... if you have now clue, it might be a very shallow extremum and even with a very small gradient you can be miles off. tldr; I don't know. Play with it.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.