>  Docs Center  >  Libraries  >  Markwardt  >  MPFITELLIPSE






  Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770
  UPDATED VERSIONs can be found on my WEB PAGE:


  Approximate fit to points forming an ellipse

Major Topics

  Curve and Surface Fitting

Calling Sequence

  parms = MPFITELLIPSE(X, Y, start_parms, [/TILT, WEIGHTS=wts, ...])


  MPFITELLIPSE fits a closed elliptical or circular curve to a two
  dimensional set of data points. The user specifies the X and Y
  positions of the points, and an optional set of weights. The
  ellipse may also be tilted at an arbitrary angle.
  IMPORTANT NOTE: this fitting program performs simple ellipse
  fitting. It will not work well for ellipse data with high
  eccentricity. More robust answers can usually be obtained with
  "orthogonal distance regression." (See FORTRAN package ODRPACK on
  netlib.org for more information).
  The best fitting ellipse parameters are returned from by
  MPFITELLIPSE as a vector, whose values are:
      P[0] Ellipse semi axis 1
      P[1] Ellipse semi axis 2 ( = P[0] if CIRCLE keyword set)
      P[2] Ellipse center - x value
      P[3] Ellipse center - y value
      P[4] Ellipse rotation angle (radians) if TILT keyword set
  If the TILT keyword is set, then the P[0] is meant to be the
  semi-major axis, and P[1] is the semi-minor axis, and P[4]
  represents the tilt of the semi-major axis with respect to the X
  axis. If the TILT keyword is not set, the P[0] and P[1] represent
  the ellipse semi-axes in the X and Y directions, respectively.
  The returned semi-axis lengths should always be positive.
  The user may specify an initial set of trial parameters, but by
  default MPFITELLIPSE will estimate the parameters automatically.
  Users should be aware that in the presence of large amounts of
  noise, namely when the measurement error becomes significant
  compared to the ellipse axis length, then the estimated parameters
  become unreliable. Generally speaking the computed axes will
  overestimate the true axes. For example when (SIGMA_R/R) becomes
  0.5, the radius of the ellipse is overestimated by about 40%.
  This unreliability is also pronounced if the ellipse has high
  eccentricity, as noted above.
  Users can weight their data as they see appropriate. However, the
  following prescription for the weighting may serve as a good
  starting point, and appeared to produce results comparable to the
  typical chi-squared value.
    WEIGHTS = 0.75/(SIGMA_X^2 + SIGMA_Y^2)
  where SIGMA_X and SIGMA_Y are the measurement error vectors in the
  X and Y directions respectively. However, this has not been
  robustly tested, and it should be pointed out that this weighting
  may only be appropriate for a set of points whose measurement
  errors are comparable. If a more robust estimation of the
  parameter values is needed, the so-called orthogonal distance
  regression package should be used (ODRPACK, available in FORTRAN
  at www.netlib.org).


  X - measured X positions of the points in the ellipse.
  Y - measured Y positions of the points in the ellipse.
  START_PARAMS - an array of starting values for the ellipse
                  parameters, as described above. This parameter is
                  optional; if not specified by the user, then the
                  ellipse parameters are estimated automatically from
                  the properties of the data.


  Returns the best fitting model ellipse parameters. Returned
  values are undefined if STATUS indicates an error condition.


  ** NOTE ** Additional keywords such as PARINFO, BESTNORM, and
              STATUS are accepted by MPFITELLIPSE but not documented
              here. Please see the documentation for MPFIT for the
              description of these advanced options.
  CIRCULAR - if set, then the curve is assumed to be a circle
              instead of ellipse. When set, the parameters P[0] and
              P[1] will be identical and the TILT keyword will have
              no effect.
  PERROR - upon return, the 1-sigma uncertainties of the returned
            ellipse parameter values. These values are only
            meaningful if the WEIGHTS keyword is specified properly.
            If the fit is unweighted (i.e. no errors were given, or
            the weights were uniformly set to unity), then PERROR
            will probably not represent the true parameter
            If STATUS indicates an error condition, then PERROR is
  QUIET - if set then diagnostic fitting messages are suppressed.
          Default: QUIET=1 (i.e., no diagnostics]
  STATUS - an integer status code is returned. All values greater
            than zero can represent success (however STATUS EQ 5 may
            indicate failure to converge). Please see MPFIT for
            the definitions of status codes.
  TILT - if set, then the major and minor axes of the ellipse
          are allowed to rotate with respect to the data axes.
          Parameter P[4] will be set to the clockwise rotation angle
          of the P[0] axis in radians, as measured from the +X axis.
          P[4] should be in the range 0 to !dpi.
  WEIGHTS - Array of weights to be used in calculating the
            chi-squared value. The chi-squared value is computed
            as follows:
                CHISQ = TOTAL( (Z-MYFUNCT(X,Y,P))^2 * ABS(WEIGHTS)^2 )
            Users may wish to follow the guidelines for WEIGHTS
            described above.


  ; Construct a set of points on an ellipse, with some noise
  ph0 = 2*!pi*randomu(seed,50)
  x = 50. + 32.*cos(ph0) + 4.0*randomn(seed, 50)
  y = -75. + 65.*sin(ph0) + 0.1*randomn(seed, 50)
  ; Compute weights function
  weights = 0.75/(4.0^2 + 0.1^2)
  ; Fit ellipse and plot result
  p = mpfitellipse(x, y)
  phi = dindgen(101)*2D*!dpi/100
  plot, x, y, psym=1
  oplot, p[2]+p[0]*cos(phi), p[3]+p[1]*sin(phi), color='ff'xl
  ; Fit ellipse and plot result - WITH TILT
  p = mpfitellipse(x, y, /tilt)
  phi = dindgen(101)*2D*!dpi/100
  ; New parameter P[4] gives tilt of ellipse w.r.t. coordinate axes
  ; We must rotate a standard ellipse to this new orientation
  xm = p[2] + p[0]*cos(phi)*cos(p[4]) + p[1]*sin(phi)*sin(p[4])
  ym = p[3] - p[0]*cos(phi)*sin(p[4]) + p[1]*sin(phi)*cos(p[4])
  plot, x, y, psym=1
  oplot, xm, ym, color='ff'xl


  MINPACK-1, Jorge More', available from netlib (www.netlib.org).
  "Optimization Software Guide," Jorge More' and Stephen Wright,
    SIAM, *Frontiers in Applied Mathematics*, Number 14.

Modification History

  Ported from MPFIT2DPEAK, 17 Dec 2000, CM
  More documentation, 11 Jan 2001, CM
  Example corrected, 18 Nov 2001, CM
  Change CIRCLE keyword to the correct CIRCULAR keyword, 13 Sep
      2002, CM
  Add error messages for SYMMETRIC and CIRCLE, 08 Nov 2002, CM
  Found small error in computation of _EVAL (when CIRCULAR) was set;
      sanity check when CIRCULAR is set, 21 Jan 2003, CM
  Convert to IDL 5 array syntax (!), 16 Jul 2006, CM
  Move STRICTARR compile option inside each function/procedure, 9
    Oct 2006
  Add disclaimer about the suitability of this program for fitting
    ellipses, 17 Sep 2007, CM
  Clarify documentation of TILT angle; make sure output contains
    semi-major axis first, followed by semi-minor; make sure that
    semi-axes are always positive (and can handle negative inputs)
      17 Sep 2007, CM
  Output tilt angle is now in range 0 to !DPI, 20 Sep 2007, CM
  Some documentation clarifications, including to remove reference
    to the "ERR" keyword, which does not exist, 17 Jan 2008, CM
  Swapping of P[0] and P[1] only occurs if /TILT is set, 06 Nov
    2009, CM
  Document an example of how to plot a tilted ellipse, 09 Nov 2009, CM
  Check for MPFIT error conditions and return immediately, 23 Jan 2010, CM

© 2022 Harris Geospatial Solutions, Inc. |  Legal
My Account    |    Store    |    Contact Us