0

I am trying to fit an ellipse over my different sets of 2D points using matplotlib. I'm using the ellipse fitting functions from here, and here's the code.

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

However, when I plot the ellipse using matplotlib, sometimes the ellipse is fitted well, and sometimes I need to rotate the ellipse 90 degrees to have it fit (the blue ellipse is rotated 90, the red ellipse is no additional rotation) Here is my code.

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)

And here are screenshots of the plots:

x1, y1 plot

x2, y2 plot

As you can see, the non-rotated (red) ellipse fits the x1,y1 data well, but the rotated ellipse (blue) fits the x2,y2 data.

I'm confused if I'm missing something here, when do I need to rotate the ellipse by 90º and when do I not need to?

Improve this question
2
  • I see the problem, I have to check the details of the algorithm to figure if there is something "obvious"....probably not today Commented May 4, 2020 at 14:49
  • Is there a specific reason why you use two different user accounts? Commented May 5, 2020 at 7:35

1 Answer 1

Reset to default
0

Here is my guess, without checking the math 100%. Looking at the definition and the Lagrangian to be solved, everything looks fine and logic. I think the vector a is correct and such are the internal parameters a to f. However, the code mentions that the function to be minimized is independent of a scaling factor. So when calculating the axis angle, one might run into an sign issue. My guess is, that this takes place in the arctan() function. One solution might be to add a sign check in the calculation of the a vector, i.e. a -> -a if a[-1] < 0.

Even better and probably also working (I just tested 2 cases from OP) is to replace

if a > c: return np.arctan( 2 * b / ( a - c ) ) / 2 else: return np.pi / 2 + np.arctan( 2 * b / (a - c ) ) / 2

by 

if a > c: return np.arctan2( 2 * b, ( a - c ) ) / 2 else: return np.pi / 2 + np.arctan2( 2 * b, ( a - c) ) / 2

If it is clear that the argument of the arctan is coming from a division, arctan2 is the one to go for.

One additional remark: I think that code with multiple return statements in one function should be omitted. In this short functions it is still easy to handle, but I don't think of it as good practice.

Update

When contacting the author of the original fit code, he mentioned that arctan2 is in use in version he provides on GitHub

This code looks much cleaner and I recommend to use that instead of the snippets from the home page.

Additional thoughts

Actually, I think it can be done in a more consistent and easier to follow way, how the parameters of the ellipse are extracted from the a vector. so I wrote this code down

# -*- 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()

I think this code should not have the 90° problem either. If removing all comments, it is quite compact and ( from the math point of view) readable.

Note

I encountered a problem in inv( S ). This matrix became singular. One workaround might be: rotate all data by a small angle and rotate the calculated t back. Also subtract the angle from the calculated angle phi.

Improve this answer
Sign up to request clarification or add additional context in comments.

Comments

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.