Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
tree: a39f276fa2
Fetching contributors…

Cannot retrieve contributors at this time

279 lines (254 sloc) 10.003 kb
;+
; 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: //ASP/Personal/mudies/scatterBrain/mpfitellipse.pro#1 $
;-
; 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
Jump to Line
Something went wrong with that request. Please try again.