;+
; 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.
;
;   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
;
;   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%.
;
;   Users can weight their data as they see appropriate.  However, the
;   following prescription for the weighting appears to be the correct
;   one, and produces values 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, it should be pointed
;   out that this weighting is only 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.
;
; 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.
;
;   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.  
;
;   QUIET - if set then diagnostic fitting messages are suppressed.
;           Default: QUIET=1 (i.e., no diagnostics)
;
;   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.
;
;   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.
;
;   WEIGHTS - Array of weights to be used in calculating the
;             chi-squared value.  If WEIGHTS is specified then the ERR
;             parameter is ignored.  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)
;   plot, x, y, psym=1
;   phi = dindgen(101)*2D*!dpi/100
;   oplot, p(2)+p(0)*cos(phi), p(3)+p(1)*sin(phi)
;
; 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
;
;  $Id: mpfitellipse.pro,v 1.8 2003/01/24 04:04:30 craigm Exp $
;-
; Copyright (C) 1997-2000,2002,2003, 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
  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

  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(p(0)*p(0)*wc)*u2 $
      else           u2 = sqrt(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

  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]]

  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

  ;; 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

  return, result
end