File: mpfit2dpeak.pro

package info (click to toggle)
mpfit 1.85+2017.01.03-4
  • links: PTS, VCS
  • area: main
  • in suites: bullseye, sid
  • size: 788 kB
  • sloc: python: 126; makefile: 16; sh: 13
file content (508 lines) | stat: -rw-r--r-- 19,441 bytes parent folder | download | duplicates (4)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
;+
; NAME:
;   MPFIT2DPEAK
;
; 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:
;   Fit a gaussian, lorentzian or Moffat model to data
;
; MAJOR TOPICS:
;   Curve and Surface Fitting
;
; CALLING SEQUENCE:
;   yfit = MPFIT2DPEAK(Z, A [, X, Y, /TILT ...] )
;
; DESCRIPTION:
;
;   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
;   MPFIT and MPFIT2DFUN).
;
;   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
;   unweighted.
;
; RESTRICTIONS:
;
;   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.
;
;   Note that if PARINFO is supplied, PARINFO(*).VALUES is ignored.
;   If you wish to supply starting values, use the ESTIMATES keyword.
;
;   MPFIT2DPEAK works in two steps.  First, it computes initial
;   ESTIMATES if none are provided, not using MPFIT.  Second, it uses
;   the initial ESTIMATES to fit a refined solution using MPFIT.  The
;   first step, initial estimates, is not required to match any
;   constraints supplied with the PARINFO keyword parameter.  Thus, if
;   you don't supply ESTIMATES but do supply PARINFO, it is possible
;   for MPFIT to fail with an error that parameters exceed their
;   PARINFO limits.  To avoid this situation, call MPFIT2DPEAK with
;   ESTIMATES explicitly.
;
; COMPATIBILITY
;
;  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.
;
;
; INPUTS:
;
;   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
;
; OUTPUTS:
;   A - Upon return, an array of best fit parameter values.  See the
;       table above for the meanings of each parameter element.
;
;
; RETURNS:
;
;   Returns the best fitting model function as a 2D array.
;
; KEYWORDS:
;
;   ** 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
;             DOF = N_ELEMENTS(DEVIATES) - NFREE
;         Note that this doesn't account for pegged parameters (see
;         NPEGGED).
;
;   ERROR - upon input, the measured 1-sigma uncertainties in the "Z"
;           values.  If no ERROR or WEIGHTS are given, then the fit is
;           unweighted.
;
;   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
;               ignored.)
;               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
;              automatically.
;
;   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
;            uncertainties.  
;
;            *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.
;
;
; EXAMPLE:
;
; ; 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)
;
; 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:
;
;   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
;   Add more documentation about the interaction of ESTIMATES and
;     PARINFO, CM, 2013-05-28
;
;  $Id: mpfit2dpeak.pro,v 1.11 2013/07/18 03:25:40 cmarkwar Exp $
;-
; Copyright (C) 1997-2000, 2002, 2003, 2005, 2006, 2007, 2008, 2013 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 mpfit, mpfitfun, mpfit2dpeak, mpfit2dpeak_gauss, $
  mpfit2dpeak_lorentz, mpfit2dpeak_moffat, mpfit2dpeak_u

; Compute the "u" value = (x/a)^2 + (y/b)^2 with optional rotation
function mpfit2dpeak_u, x, y, p, tilt=tilt, symmetric=sym
  COMPILE_OPT strictarr
  widx  = abs(p[2]) > 1e-20 & widy  = abs(p[3]) > 1e-20 
  if keyword_set(sym) then widy = widx
  xp    = x-p[4]            & yp    = y-p[5]
  theta = p[6]
  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

; Gaussian Function
function mpfit2dpeak_gauss, x, y, p, tilt=tilt, symmetric=sym, _extra=extra
  COMPILE_OPT strictarr
  sz = size(x)
  if sz[sz[0]+1] EQ 5 then smax = 26D else smax = 13.

  u = mpfit2dpeak_u(x, y, p, tilt=keyword_set(tilt), symmetric=keyword_set(sym))
  mask = u LT (smax^2)  ;; Prevents floating underflow
  return, p[0] + p[1] * mask * exp(-0.5 * u * mask)
end

; Lorentzian Function
function mpfit2dpeak_lorentz, x, y, p, tilt=tilt, symmetric=sym, _extra=extra
  COMPILE_OPT strictarr
  u = mpfit2dpeak_u(x, y, p, tilt=keyword_set(tilt), symmetric=keyword_set(sym))
  return, p[0] + p[1] / (u + 1)
end

; Moffat Function
function mpfit2dpeak_moffat, x, y, p, tilt=tilt, symmetric=sym, _extra=extra
  COMPILE_OPT strictarr
  u = mpfit2dpeak_u(x, y, p, tilt=keyword_set(tilt), symmetric=keyword_set(sym))
  return, p[0] + p[1] / (u + 1)^p[7]
end

function mpfit2dpeak, z, a, x, y, estimates=est, tilt=tilt, $
                      gaussian=gauss, lorentzian=lorentz, moffat=moffat, $
                      perror=perror, sigma=sigma, zerror=zerror, $
                      chisq=chisq, bestnorm=bestnorm, niter=iter, nfev=nfev, $
                      error=dz, weights=weights, measure_errors=dzm, $
                      nfree=nfree, dof=dof, $
                      negative=neg, parinfo=parinfo, $
                      circular=sym, circle=badcircle1, symmetric=badcircle2, $
                      errmsg=errmsg, status=status, $
                      query=query, quiet=quiet, _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 functions MPFIT and MPFIT2DFUN ' + $
        'must be in your IDL path', /info
      return, !values.d_nan
  endif
  if mpfit(/query)      NE 1 then goto, MPFIT_NOTFOUND
  if mpfit2dfun(/query) NE 1 then goto, MPFIT_NOTFOUND
  catch, /cancel
  if keyword_set(query) then return, 1

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

  ;; Reject too few data
  if n_elements(z) LT 8 then begin
      message, 'ERROR: array must have at least eight elements', /info
      return, !values.d_nan
  endif

  sz = size(z)
  if sz[0] LT 2 then begin
      message, 'ERROR: array must be 2-dimensional', /info
      return, !values.d_nan
  endif
  nx = sz[1]
  ny = sz[2]
  
  ;; Fill in the X and Y values if needed -- note clever promotion to
  ;;                                         double if needed
  if n_elements(x) EQ 0 then x = findgen(nx)*(z[0]*0+1)
  if n_elements(y) EQ 0 then y = findgen(ny)*(z[0]*0+1)
  if n_elements(x) LT nx then begin
      message, 'ERROR: X array was not large enough', /info
      return, !values.d_nan
  endif
  if n_elements(y) LT ny then begin
      message, 'ERROR: Y array was not large enough', /info
      return, !values.d_nan
  endif

  ;; Make 2D arrays of X and Y values -- if the user hasn't done it
  if n_elements(x) NE n_elements(z) then   xx = x[*] # (y[*]*0 + 1)   else xx = x
  if n_elements(y) NE n_elements(z) then   yy = (x[*]*0 + 1) # y[*]   else yy = y

  ;; Compute the weighting factors to use
  if (n_elements(dz) EQ 0 AND n_elements(weights) EQ 0 AND $
      n_elements(dzm) EQ 0) then begin
      weights = z*0+1        ;; Unweighted by default
  endif else if n_elements(dz) GT 0 then begin
      weights = dz * 0   ;; Avoid division by zero
      wh = where(dz NE 0, ct)
      if ct GT 0 then weights[wh] = 1./dz[wh]^2
  endif else if n_elements(dzm) GT 0 then begin
      weights = dzm * 0   ;; Avoid division by zero
      wh = where(dzm NE 0, ct)
      if ct GT 0 then weights[wh] = 1./dzm[wh]^2
  endif

  if n_elements(est) EQ 0 then begin
      ;; Here is the secret - the width is estimated based on the volume
      ;; above/below the average.  Thus, as the signal becomes more
      ;; noisy the width automatically broadens as it should.

      maxx = max(x, min=minx) & maxy = max(y, min=miny)
      maxz = max(z, whmax)    & minz = min(z, whmin)
      nx = n_elements(x)
      dx = 0.5 * [x[1]-x[0], x[2:*] - x, x[nx-1] - x[nx-2]]
      ny = n_elements(y)
      dy = 0.5 * [y[1]-y[0], y[2:*] - y, y[ny-1] - y[ny-2]]
      ;; Compute cell areas
      da = dx # dy

      totvol = total(da*z)               ;; Total volume under curve
      av = totvol/(total(dx)*total(dy))  ;; Average height

      ;; Compute the spread in values above and below average... we
      ;; take the narrowest one as the one with the peak
      wh   = where(z GE av, ct1)
      sdx1 = total(xx[wh]^2)/ct1 - (total(xx[wh])/ct1)^2
      sdy1 = total(yy[wh]^2)/ct1 - (total(yy[wh])/ct1)^2
      wh   = where(z LE av, ct2)
      sdx2 = total(xx[wh]^2)/ct2 - (total(xx[wh])/ct2)^2
      sdy2 = total(yy[wh]^2)/ct2 - (total(yy[wh])/ct2)^2
      wh   = 0  ;; conserve memory

      if keyword_set(pos) then goto, POS_PEAK
      if keyword_set(neg) then goto, NEG_PEAK

      ;; Compute volume above/below average
      if (sdx1 LT sdx2 AND sdy1 LT sdy2) then begin 
          ;; Positive peak
          POS_PEAK:
          centx = xx[whmax]
          centy = yy[whmax]
          peakz = maxz - av
      endif else if (sdx1 GT sdx2 AND sdy1 GT sdy2) then begin 
          ;; Negative peak
          NEG_PEAK:
          centx = xx[whmin]
          centy = yy[whmin]
          peakz = minz - av
      endif else begin
          ;; Ambiguous case
          if abs(maxz - av) GT (minz - av) then goto, POS_PEAK $
          else                                  goto, NEG_PEAK
      endelse
      peakvol = totvol - total(da*(z<av))
      width   = sqrt(peakvol / (6*abs(peakz)))

      est = [av, peakz, width, width, centx, centy, 0, 1]
      guess = 1
  endif

  ;; Check the number of parameter estimates
  if n_elements(quiet) EQ 0 then quiet=1

  np = 7
  ;; Parameter checking for individual function types
  if keyword_set(moffat) then begin               ;; MOFFAT
      fun = 'mpfit2dpeak_moffat'
      np = 8
  endif else if keyword_set(lorentz) then begin  ;; LORENTZIAN
      fun = 'mpfit2dpeak_lorentz'
  endif else begin                               ;; GAUSSIAN
      fun = 'mpfit2dpeak_gauss'
  endelse
  if n_elements(est) LT np then begin
      message, 'ERROR: parameter ESTIMATES must have at least '+strtrim(np,2)+$
        ' elements', /info
      return, !values.d_nan
  endif
  p0 = replicate(est[0]*0, np > n_elements(est))
  p0[0] = est

  ;; Function call
  fargs = {tilt: keyword_set(tilt), symmetric: keyword_set(sym)}
  a = mpfit2dfun(fun, xx, yy, z, 0, p0[0:np-1], weights=weights, $
                 bestnorm=bestnorm, nfev=nfev, status=status, $
                 parinfo=parinfo, perror=perror, niter=iter, yfit=yfit, $
                 quiet=quiet, errmsg=errmsg, nfree=nfree, dof=dof, $
                 functargs=fargs, _EXTRA=extra)

  ;; Print error message if there is one.
  if NOT keyword_set(quiet) AND errmsg NE '' then $
    message, errmsg, /info

  ;; Make sure the parameters are sane
  if status NE 0 then begin
      ;; Width is positive
      a[2] = abs(a[2])
      a[3] = abs(a[3])
      if keyword_set(sym) then a[3] = a[2]
      ;; Make sure that theta is in the range [0,pi]
      a[6] = ((a[6] MOD !dpi) + 2*!dpi) MOD !dpi
      a = a[0:np-1]

      if n_elements(perror)   GT 0 then sigma = perror
      if n_elements(bestnorm) GT 0 then chisq = bestnorm

      if n_elements(yfit) EQ nx*ny then begin
          yfit = reform(yfit, nx, ny, /overwrite)
      endif

      zerror = a[0]*0
      if n_elements(dof) GT 0 AND dof[0] GT 0 then begin
          zerror[0] = sqrt( total( (z-yfit)^2 ) / dof[0] )
      endif

      return, yfit
  endif

  return, !values.d_nan
end