text: import numpy
  import matplotlib.pyplot as plt
  from matplotlib.patches import Ellipse
  
  """Demonstration of least-squares fitting of ellipses
  
      __author__ = "Ben Hammel, Nick Sullivan-Molina"
      __credits__ = ["Ben Hammel", "Nick Sullivan-Molina"]
      __maintainer__ = "Ben Hammel"
      __email__ = "bdhammel@gmail.com"
      __status__ = "Development"
  
      Requirements 
      ------------
      Python 2.X or 3.X
      numpy
      matplotlib
  
      References
      ----------
      (*) Halir, R., Flusser, J.: 'Numerically Stable Direct Least Squares 
          Fitting of Ellipses'
      (**) http://mathworld.wolfram.com/Ellipse.html
      (***) White, A. McHale, B. 'Faraday rotation data analysis with least-squares 
          elliptical fitting'
  """
  
  class LSqEllipse:
  
      def fit(self, data):
          """Lest Squares fitting algorithm 
  
          Theory taken from (*)
          Solving equation Sa=lCa. with a = |a b c d f g> and a1 = |a b c> 
              a2 = |d f g>
  
          Args
          ----
          data (list:list:float): list of two lists containing the x and y data of the
              ellipse. of the form [[x1, x2, ..., xi],[y1, y2, ..., yi]]
  
          Returns
          ------
          coef (list): list of the coefficients describing an ellipse
             [a,b,c,d,f,g] corresponding to ax**2+2bxy+cy**2+2dx+2fy+g
          """
          # 弄一个栈吧 然后for循环 这样能实现continue 当cond<0时候能返回
          x, y = numpy.asarray(data, dtype=float)
          # stack=[]
          # stack.append([x,y])
          # while stack:
          # x,y=stack.pop()
          #Quadratic part of design matrix [eqn. 15] from (*)
          # reference:https://blog.csdn.net/liucc09/article/details/82763635
          D1 = numpy.mat(numpy.vstack([x**2, x*y, y**2])).T
          #Linear part of design matrix [eqn. 16] from (*)
          D2 = numpy.mat(numpy.vstack([x, y, numpy.ones(len(x))])).T
  
          #forming scatter matrix [eqn. 17] from (*)
          S1 = D1.T*D1
          S2 = D1.T*D2
          S3 = D2.T*D2
  
          #Constraint matrix [eqn. 18]
          C1 = numpy.mat('0. 0. 2.; 0. -1. 0.; 2. 0. 0.')
  
          #Reduced scatter matrix [eqn. 29]
          M=C1.I*(S1-S2*S3.I*S2.T)
  
          #M*|a b c >=l|a b c >. Find eigenvalues(特征值) and eigenvectors(特征向量) from this equation [eqn. 28]
          eval, evec = numpy.linalg.eig(M)
  
          # eigenvector must meet constraint 4ac - b^2 to be valid.
  
          cond = 4*numpy.multiply(evec[0, :], evec[2, :]) - numpy.power(evec[1, :], 2)
          # print(cond)
          # if cond[0,0]<0:
          #     break
  
              # cond[0,0] = abs(cond[0,0]) # 强行拟合
          a1 = evec[:, numpy.nonzero(cond.A > 0)[1]]
  
          #|d f g> = -S3^(-1)*S2^(T)*|a b c> [eqn. 24]
          a2 = -S3.I*S2.T*a1
  
          # eigenvectors |a b c d f g>
          self.coef = numpy.vstack([a1, a2])
          self._save_parameters()
  
      def _save_parameters(self):
          """finds the important parameters of the fitted ellipse
  
          Theory taken form http://mathworld.wolfram
  
          Args
          -----
          coef (list): list of the coefficients describing an ellipse
             [a,b,c,d,f,g] corresponding to ax**2+2bxy+cy**2+2dx+2fy+g
  
          Returns
          _______
          center (List): of the form [x0, y0]
          width (float): major axis
          height (float): minor axis
          phi (float): rotation of major axis form the x-axis in radians
          """
  
          #eigenvectors are the coefficients of an ellipse in general form
          #a*x^2 + 2*b*x*y + c*y^2 + 2*d*x + 2*f*y + g = 0 [eqn. 15) from (**) or (***)
          a = self.coef[0,0]
          # print(a)
          b = self.coef[1,0]/2.
          c = self.coef[2,0]
          d = self.coef[3,0]/2.
          f = self.coef[4,0]/2.
          g = self.coef[5,0]
  
          #finding center of ellipse [eqn.19 and 20] from (**)
          x0 = (c*d-b*f)/(b**2.-a*c)
          y0 = (a*f-b*d)/(b**2.-a*c)
  
          #Find the semi-axes lengths [eqn. 21 and 22] from (**)
          numerator = 2*(a*f*f+c*d*d+g*b*b-2*b*d*f-a*c*g)
          # numerator = abs(numerator)
          denominator1 = (b*b-a*c)*( (c-a)*numpy.sqrt(1+4*b*b/((a-c)*(a-c)))-(c+a))
          # denominator1 = abs(denominator1)
          denominator2 = (b*b-a*c)*( (a-c)*numpy.sqrt(1+4*b*b/((a-c)*(a-c)))-(c+a))
          # denominator2 = abs(denominator2)
          # TODO:有些denominator1 denominator2 numerator小于零
          # if denominator1>0 and denominator2>0 and numerator>0:
          width = numpy.sqrt(numerator/denominator1)
          height = numpy.sqrt(numerator/denominator2)
          # angle of counterclockwise rotation of major-axis of ellipse to x-axis [eqn. 23] from (**)
          # or [eqn. 26] from (***).
          phi = .5*numpy.arctan((2.*b)/(a-c))
  
          self._center = [x0, y0]
          self._width = width
          self._height = height
          self._phi = phi
          # print([x0,y0])
      @property
      def center(self):
          return self._center
  
      @property
      def width(self):
          return self._width
  
      @property
      def height(self):
          return self._height
  
      @property
      def phi(self):
          """angle of counterclockwise rotation of major-axis of ellipse to x-axis
          [eqn. 23] from (**)
          """
          return self._phi
  
      def parameters(self):
          return self.center, self.width, self.height, self.phi
  
  
  def make_test_ellipse(center=[1,1], width=1, height=.6, phi=3.14/5):
      """Generate Elliptical data with noise
      
      Args
      ----
      center (list:float): (<x_location>, <y_location>)
      width (float): semimajor axis. Horizontal dimension of the ellipse (**)
      height (float): semiminor axis. Vertical dimension of the ellipse (**)
      phi (float:radians): tilt of the ellipse, the angle the semimajor axis
          makes with the x-axis 
  
      Returns
      -------
      data (list:list:float): list of two lists containing the x and y data of the
          ellipse. of the form [[x1, x2, ..., xi],[y1, y2, ..., yi]]
      """
      t = numpy.linspace(0, 2*numpy.pi, 1000)
      x_noise, y_noise = numpy.random.rand(2, len(t))
      
      ellipse_x = center[0] + width*numpy.cos(t)*numpy.cos(phi)-height*numpy.sin(t)*numpy.sin(phi) + x_noise/2.
      ellipse_y = center[1] + width*numpy.cos(t)*numpy.sin(phi)+height*numpy.sin(t)*numpy.cos(phi) + y_noise/2.
  
      return [ellipse_x, ellipse_y]
  
  