;+
; NAME:
;   MPFITELLIPSE
;
; AUTHOR:
;   Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770
;   craigm@lheamail.gsfc.nasa.gov
;   UPDATED VERSIONs can be found on my WEB PAGE: 
;      http://cow.physics.wisc.edu/~craigm/idl/idl.html
;
; PURPOSE:
;   Approximate fit to points forming an ellipse
;
; MAJOR TOPICS:
;   Curve and Surface Fitting
;
; CALLING SEQUENCE:
;   parms = MPFITELLIPSE(X, Y, start_parms, [/TILT, WEIGHTS=wts, ...])
;
; DESCRIPTION:
;
;   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).
;
; INPUTS:
;
;   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:
;
;   Returns the best fitting model ellipse parameters.  Returned
;   values are undefined if STATUS indicates an error condition.
;
; KEYWORDS:
;
;   ** 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
;            uncertainties.
;
;            If STATUS indicates an error condition, then PERROR is
;            undefined.
;
;   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.
;
;
; EXAMPLE:
;
; ; 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
;
; REFERENCES:
;
;   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
;
;  $Id: mpfitellipse.pro,v 1.14 2010/01/25 03:38:03 craigm Exp $
;-
; Copyright (C) 1997-2000,2002,2003,2007,2008,2009,2010 Craig Markwardt
; This software is provided as is without any warranty whatsoever.
; Permission to use, copy, modify, and distribute modified or
; unmodified copies is granted, provided this copyright and disclaimer
; are included unchanged.
;-


FORWARD_FUNCTION mpfitellipse_u, mpfitellipse_eval, mpfitellipse, mpfit

; Compute the "u" value = (x/a)^2 + (y/b)^2 with optional rotation
function mpfitellipse_u, x, y, p, tilt=tilt, circle=circle
  COMPILE_OPT strictarr
  widx  = abs(p[0]) > 1e-20 & widy  = abs(p[1]) > 1e-20 
  if keyword_set(circle) then widy  = widx
  xp    = x-p[2]            & yp    = y-p[3]
  theta = p[4]

  if keyword_set(tilt) AND theta NE 0 then begin
      c  = cos(theta) & s  = sin(theta)
      return, ( (xp * (c/widx) - yp * (s/widx))^2 + $
                (xp * (s/widy) + yp * (c/widy))^2 )
  endif else begin
      return, (xp/widx)^2 + (yp/widy)^2
  endelse

end

; This is the call-back function for MPFIT.  It evaluates the
; function, subtracts the data, and returns the residuals.
function mpfitellipse_eval, p, tilt=tilt, circle=circle, _EXTRA=extra

  COMPILE_OPT strictarr
  common mpfitellipse_common, xy, wc

  tilt = keyword_set(tilt) 
  circle = keyword_set(circle)
  u2 = mpfitellipse_u(xy[*,0], xy[*,1], p, tilt=tilt, circle=circle) - 1.

  if n_elements(wc) GT 0 then begin
      if circle then u2 = sqrt(abs(p[0]*p[0]*wc))*u2 $
      else           u2 = sqrt(abs(p[0]*p[1]*wc))*u2 
  endif

  return, u2
end

function mpfitellipse, x, y, p0, WEIGHTS=wts, $
                       BESTNORM=bestnorm, nfev=nfev, STATUS=status, $
                       tilt=tilt, circular=circle, $
                       circle=badcircle1, symmetric=badcircle2, $
                       parinfo=parinfo, query=query, $
                       covar=covar, perror=perror, niter=iter, $
                       quiet=quiet, ERRMSG=errmsg, _EXTRA=extra

  COMPILE_OPT strictarr
  status = 0L
  errmsg = ''

  ;; Detect MPFIT and crash if it was not found
  catch, catcherror
  if catcherror NE 0 then begin
      MPFIT_NOTFOUND:
      catch, /cancel
      message, 'ERROR: the required function MPFIT must be in your IDL path', /info
      return, !values.d_nan
  endif
  if mpfit(/query) NE 1 then goto, MPFIT_NOTFOUND
  catch, /cancel
  if keyword_set(query) then return, 1

  if n_params() EQ 0 then begin
      message, "USAGE: PARMS = MPFITELLIPSE(X, Y, START_PARAMS, ... )", $
        /info
      return, !values.d_nan
  endif
  nx = n_elements(x) & ny = n_elements(y)
  if (nx EQ 0) OR (ny EQ 0) OR (nx NE ny) then begin
      message, 'ERROR: X and Y must have the same number of elements', /info
      return, !values.d_nan
  endif

  if keyword_set(badcircle1) OR keyword_set(badcircle2) then $
    message, 'ERROR: do not use the CIRCLE or SYMMETRIC keywords.  ' +$
    'Use CIRCULAR instead.'

  p = make_array(5, value=x[0]*0)

  if n_elements(p0) GT 0 then begin
      p[0] = p0
      if keyword_set(circle) then p[1] = p[0]
  endif else begin
      mx = moment(x)
      my = moment(y)
      p[0] = [sqrt(mx[1]), sqrt(my[1]), mx[0], my[0], 0]
      if keyword_set(circle) then $
        p[0:1] = sqrt(mx[1]+my[1])
  endelse

  common mpfitellipse_common, xy, wc
  if n_elements(wts) GT 0 then begin
      wc = abs(wts)
  endif else begin
      wc = 0 & dummy = temporary(wc)
  endelse

  xy = [[x],[y]]

  nfev = 0L & dummy = temporary(nfev)
  covar = 0 & dummy = temporary(covar)
  perror = 0 & dummy = temporary(perror)
  status = 0
  result = mpfit('mpfitellipse_eval', p, $
                 parinfo=parinfo, STATUS=status, nfev=nfev, BESTNORM=bestnorm,$
                 covar=covar, perror=perror, niter=iter, $
                 functargs={circle:keyword_set(circle), tilt:keyword_set(tilt)},$
                 ERRMSG=errmsg, quiet=quiet, _EXTRA=extra)

  ;; Print error message if there is one.
  if NOT keyword_set(quiet) AND errmsg NE '' then $
    message, errmsg, /info
  ;; Return if there is an error condition
  if status LE 0 then return, result

  ;; Sanity check on resulting parameters
  if keyword_set(circle) then begin
      result[1] = result[0]
      perror[1] = perror[0]
  endif
  if NOT keyword_set(tilt) then begin
      result[4] = 0
      perror[4] = 0
  endif

  ;; Make sure the axis lengths are positive, and the semi-major axis
  ;; is listed first
  result[0:1] = abs(result[0:1])
  if abs(result[0]) LT abs(result[1]) AND keyword_set(tilt) then begin
      tmp = result[0] & result[0] = result[1] & result[1] = tmp
      tmp = perror[0] & perror[0] = perror[1] & perror[1] = tmp
      result[4] = result[4] - !dpi/2d
  endif

  if keyword_set(tilt) then begin
      ;; Put tilt in the range 0 to +Pi
      result[4] = result[4] - !dpi * floor(result[4]/!dpi)
  endif

  return, result
end