>  Docs Center  >  Libraries  >  Markwardt  >  MPFIT2DPEAK






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


  Fit a gaussian, lorentzian or Moffat model to data

Major Topics

  Curve and Surface Fitting

Calling Sequence

  yfit = MPFIT2DPEAK(Z, A [, X, Y, /TILT ...] )


  MPFIT2DPEAK fits a gaussian, lorentzian or Moffat model using the
  non-linear least squares fitter MPFIT. MPFIT2DPEAK is meant to be
  a drop-in replacement for IDL's GAUSS2DFIT function (and requires
  The choice of the fitting function is determined by the keywords
  GAUSSIAN, LORENTZIAN and MOFFAT. By default the gaussian model
  function is used. [ The Moffat function is a modified Lorentzian
  with variable power law index. ] The two-dimensional peak has
  independent semimajor and semiminor axes, with an optional
  rotation term activated by setting the TILT keyword. The baseline
  is assumed to be a constant.
      GAUSSIAN A[0] + A[1]*exp(-0.5*u)
      LORENTZIAN A[0] + A[1]/(u + 1)
      MOFFAT A[0] + A[1]/(u + 1)^A[7]
      u = ( (x-A[4])/A[2] )^2 + ( (y-A[5])/A[3] )^2
        where x and y are cartesian coordinates in rotated
        coordinate system if TILT keyword is set.
  The returned parameter array elements have the following meanings:
      A[0] Constant baseline level
      A[1] Peak value
      A[2] Peak half-width (x) -- gaussian sigma or half-width at half-max
      A[3] Peak half-width (y) -- gaussian sigma or half-width at half-max
      A[4] Peak centroid (x)
      A[5] Peak centroid (y)
      A[6] Rotation angle (radians) if TILT keyword set
      A[7] Moffat power law index if MOFFAT keyword set
  By default the initial starting values for the parameters A are
  estimated from the data. However, explicit starting values can be
  supplied using the ESTIMATES keyword. Also, error or weighting
  values can optionally be provided; otherwise the fit is


  If no starting parameter ESTIMATES are provided, then MPFIT2DPEAK
  attempts to estimate them from the data. This is not a perfect
  science; however, the author believes that the technique
  implemented here is more robust than the one used in IDL's
  GAUSS2DFIT. The author has tested cases of strong peaks, noisy
  peaks and broad peaks, all with success.


  This function is designed to work with IDL 5.0 or greater.
  Because TIED parameters rely on the EXECUTE() function, they cannot
  be used with the free version of the IDL Virtual Machine.


  Z - Two dimensional array of "measured" dependent variable values.
      Z should be of the same type and dimension as (X # Y).
  X - Optional vector of x positions for a single row of Z.
          X[i] should provide the x position of Z[i,*]
      Default: X values are integer increments from 0 to NX-1
  Y - Optional vector of y positions for a single column of Z.
          Y[j] should provide the y position of Z[*,j]
      Default: Y values are integer increments from 0 to NY-1


  A - Upon return, an array of best fit parameter values. See the
      table above for the meanings of each parameter element.


  Returns the best fitting model function as a 2D array.


  ** NOTE ** Additional keywords such as PARINFO, BESTNORM, and
              STATUS are accepted by MPFIT2DPEAK but not documented
              here. Please see the documentation for MPFIT for the
              description of these advanced options.
  CHISQ - the value of the summed squared residuals for the
          returned parameter values.
  CIRCULAR - if set, then the peak profile is assumed to be
              azimuthally symmetric. When set, the parameters A[2)
              and A[3) will be identical and the TILT keyword will
              have no effect.
  DOF - number of degrees of freedom, computed as
        Note that this doesn't account for pegged parameters (see
  ERROR - upon input, the measured 1-sigma uncertainties in the "Z"
          values. If no ERROR or WEIGHTS are given, then the fit is
  ESTIMATES - Array of starting values for each parameter of the
              model. If ESTIMATES is not set, then the starting
              values are estimated from the data directly, before
              fitting. (This also means that PARINFO.VALUES is
              Default: not set - parameter values are estimated from data.
  GAUSSIAN - if set, fit a gaussian model function. The Default.
  LORENTZIAN - if set, fit a lorentzian model function.
  MOFFAT - if set, fit a Moffat model function.
  MEASURE_ERRORS - synonym for ERRORS, for consistency with built-in
                    IDL fitting routines.
  NEGATIVE - if set, and ESTIMATES is not provided, then MPFIT2DPEAK
              will assume that a negative peak is present -- a
              valley. Specifying this keyword is not normally
              required, since MPFIT2DPEAK can determine this
  NFREE - the number of free parameters in the fit. This includes
          parameters which are not FIXED and not TIED, but it does
          include parameters which are pegged at LIMITS.
  PERROR - upon return, the 1-sigma uncertainties of the parameter
            values A. These values are only meaningful if the ERRORS
            or WEIGHTS keywords are 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* you can assume that the true reduced chi-squared
            value is unity -- meaning that the fit is implicitly
            assumed to be of good quality -- then the estimated
            parameter uncertainties can be computed by scaling PERROR
            by the measured chi-squared value.
              DOF = N_ELEMENTS(Z) - N_ELEMENTS(A) ; deg of freedom
              PCERROR = PERROR * SQRT(BESTNORM / DOF) ; scaled uncertainties
  QUIET - if set then diagnostic fitting messages are suppressed.
          Default: QUIET=1 (i.e., no diagnostics)
  SIGMA - synonym for PERROR (1-sigma parameter uncertainties), for
          compatibility with GAUSSFIT. Do not confuse this with the
          Gaussian "sigma" width parameter.
  TILT - if set, then the major and minor axes of the peak profile
          are allowed to rotate with respect to the image axes.
          Parameter A[6] will be set to the clockwise rotation angle
          of the A[2] 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) )
            Here are common values of WEIGHTS:
                1D/ERR^2 - Normal weighting (ERR is the measurement error)
                1D/Y - Poisson weighting (counting statistics)
                1D - Unweighted
            The ERROR keyword takes precedence over any WEIGHTS
            keyword values. If no ERROR or WEIGHTS are given, then
            the fit is unweighted.


  ; Construct a sample gaussian surface in range [-5,5] centered at [2,-3]
  x = findgen(100)*0.1 - 5. & y = x
  xx = x # (y*0 + 1)
  yy = (x*0 + 1) # y
  rr = sqrt((xx-2.)^2 + (yy+3.)^2)
  ; Gaussian surface with sigma=0.5, peak value of 3, noise with sigma=0.2
  z = 3.*exp(-(rr/0.5)^2) + randomn(seed,100,100)*.2
  ; Fit gaussian parameters A
  zfit = mpfit2dpeak(z, a, x, y)


  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

  New algorithm for estimating starting values, CM, 31 Oct 1999
  Documented, 02 Nov 1999
  Small documentation fixes, 02 Nov 1999
  Documented PERROR for unweighted fits, 03 Nov 1999, CM
  Copying permission terms have been liberalized, 26 Mar 2000, CM
  Small cosmetic changes, 21 Sep 2000, CM
  Corrected bug introduced by cosmetic changes, 11 Oct 2000, CM :-)
  Added POSITIVE keyword, 17 Nov 2000, CM
  Removed TILT in common, in favor of FUNCTARGS approach, 23 Nov
    2000, CM
  Added SYMMETRIC keyword, documentation for TILT, and an example,
    24 Nov 2000, CM
  Changed SYMMETRIC to CIRCULAR, 17 Dec 2000, CM
  Really change SYMMETRIC to CIRCULAR!, 13 Sep 2002, CM
  Add error messages for SYMMETRIC and CIRCLE, 08 Nov 2002, CM
  Make more consistent with comparable IDL routines, 30 Jun 2003, CM
  Defend against users supplying strangely dimensioned X and Y, 29
    Jun 2005, CM
  Convert to IDL 5 array syntax (!), 16 Jul 2006, CM
  Move STRICTARR compile option inside each function/procedure, 9 Oct 2006
  Add COMPATIBILITY section, CM, 13 Dec 2007
  Clarify documentation regarding what happens when ESTIMATES is not
    set, CM, 14 Dec 2008

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