title: Fitting an ellipse to a set of data points in python
body: <p>I have a 2D points (x,y), and I want to fit the ellipse using this post </p>
  
  <p><a href="https://stackoverflow.com/questions/13635528/fit-a-ellipse-in-python-given-a-set-of-points-xi-xi-yi">fit a ellipse in Python given a set of points xi=(xi,yi)</a></p>
  
  <p>But my result is <code>axes =  [ 0.93209407 nan]</code> since in function <code>ellipse_axis_length</code> the <code>down2</code> is a minus number, so res2 is invalid, how to do with this? and if I want to draw an ellipse according to the dataset, and calculate the error between the data points and the ellipse, how could I do?</p>
  
  <p>and the code is like this:</p>
  
  <pre><code>import numpy as np
  import numpy.linalg as linalg
  import matplotlib.pyplot as plt
  
  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 =  linalg.eig(np.dot(linalg.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_angle_of_rotation( a ):
      b,c,d,f,g,a = a[1]/2, a[2], a[3]/2, a[4]/2, a[5], a[0]
      return 0.5*np.arctan(2*b/(a-c))
  
  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 find_ellipse(x, y):
      xmean = x.mean()
      ymean = y.mean()
      x = x - xmean
      y = y - ymean
      a = fitEllipse(x,y)
      center = ellipse_center(a)
      center[0] += xmean
      center[1] += ymean
      phi = ellipse_angle_of_rotation(a)
      axes = ellipse_axis_length(a)
      x += xmean
      y += ymean
      return center, phi, axes
  
  if __name__ == '__main__':
  
      points = [( 0 , 3),
          ( 1 , 2),
          ( 1 , 7),
          ( 2 , 2),
          ( 2 , 4),
          ( 2 , 5),
          ( 2 , 6),
          ( 2 ,14),
          ( 3 , 4),
          ( 4 , 4),
          ( 5 , 5),
          ( 5 ,14),
          ( 6 , 4),
          ( 7 , 3),
          ( 7 , 7),
          ( 8 ,10),
          ( 9 , 1),
          ( 9 , 8),
          ( 9 , 9),
          (10,  1),
          (10,  2),
          (10 ,12),
          (11 , 0),
          (11 , 7),
          (12 , 7),
          (12 ,11),
          (12 ,12),
          (13 , 6),
          (13 , 8),
          (13 ,12),
          (14 , 4),
          (14 , 5),
          (14 ,10),
          (14 ,13)]
  
      fig, axs = plt.subplots(2, 1, sharex = True, sharey = True)
      a_points = np.array(points)
      x = a_points[:, 0]
      y = a_points[:, 1]
      axs[0].plot(x,y)
      center, phi, axes = find_ellipse(x, y)
      print "center = ",  center
      print "angle of rotation = ",  phi
      print "axes = ", axes
  
      axs[1].plot(x, y)
      axs[1].scatter(center[0],center[1], color = 'red', s = 100)
      axs[1].set_xlim(x.min(), x.max())
      axs[1].set_ylim(y.min(), y.max())
  
      plt.show()
  </code></pre>