GAUSS2DFIT

The GAUSS2DFIT function fits a two-dimensional, elliptical Gaussian equation to rectilinearly gridded data. Given a set of gridded data, Z = F(x, y),

where:

And the elliptical function is:

The rotated coordinate system is defined as:

The parameters of the ellipse U are:

The rotation is optional, and can be forced to 0, making the major and minor axes of the ellipse parallel to the X and Y axes.

Procedure Used and Other Notes

The peak/valley is found by first smoothing Z and then finding the maximum or minimum, respectively. GAUSSFIT is then applied to the row and column running through the peak/valley to estimate the parameters of the Gaussian in X and Y. Finally, CURVEFIT is used to fit the 2D Gaussian to the data.

Be sure that the 2D array to be fit contains the entire peak/valley out to at least 5 to 8 half-widths, or the curve-fitter may not converge.

This is a computationally-intensive routine. The time required is roughly proportional to the number of elements in Z.

This routine is written in the IDL language. Its source code can be found in the file gauss2dfit.pro in the lib subdirectory of the IDL distribution.

Syntax

Result = GAUSS2DFIT( Z, A [, X, Y] [, FITA=vector] [, MASK=array] [, /NEGATIVE] [, /TILT] )

Return Value

The result is a two-dimensional array with the same dimensions as Z that contains the function F(x, y) evaluated at each X and Y value. The coefficients of the fit are returned in the argument A.

Arguments

Z

The dependent variable. Z should be a two-dimensional array with dimensions (Nx, Ny). Gridding in the array must be rectilinear.

A

A named variable in which the coefficients of the fit are returned. A is returned as a seven element vector the coefficients of the fitted function. The meanings of the seven elements in relation to the discussion above is:

X

An optional vector with Nx elements that contains the X values of Z (i.e., Xi is the X value for Zi,j. If this argument is omitted, a regular grid in X is assumed, and the X location of Zi,j = i.

Y

An optional vector with Ny elements that contains the Y values of Z (i.e., Yj is the Y value for Zi,j. If this argument is omitted, a regular grid in Y is assumed, and the Y location of Zi,j = j.

Keywords

FITA

Set this keyword to a vector with the same number of elements as A. The vector should contain a zero for each fixed parameter and non-zero values for the elements of A to fit. If you do not specify FITA, then all parameters are assumed to be non-fixed.

MASK

Set this keyword to a 2D array with Nx x Ny dimensions, containing a 1 for each good data point, and a 0 for each missing data point. By default all data points are included in the fit.

NEGATIVE

Set this keyword to indicate that the Gaussian to be fitted is a valley (such as an absorption line). By default, a peak is fit.

TILT

Set this keyword to allow the orientation of the major and minor axes of the ellipse to be unrestricted. The default is that the axes of the ellipse must be parallel to the X and Y axes. Therefore, in the default case, A[6] is always returned as 0.

Examples

This example creates a 2D gaussian, adds random noise, adds some bad data points, and then applies GAUSS2DFIT.

; Create X and Y arrays:

nx = 128L

ny = 100L

X = FINDGEN(nx) # REPLICATE(1.0, ny)

Y = REPLICATE(1.0, nx) # FINDGEN(ny)

 

; Define input function parameters:

aAxis = nx/6.

bAxis = ny/10.

h = 0.5*nx

k = 0.6*ny

tilt = 30*!PI/180

A = [ 5., 10., aAxis, bAxis, h, k, tilt]

 

; Create an ellipse:

xprime = (X - h)*cos(tilt) - (Y - k)*sin(tilt)

yprime = (X - h)*sin(tilt) + (Y - k)*cos(tilt)

U = (xprime/aAxis)^2 + (yprime/bAxis)^2

 

; Create gaussian Z with random noise:

Zideal = A[0] + A[1] * EXP(-U/2)

Z = Zideal + RANDOMN(seed, nx, ny)

B = [] ; clear out the variable

 

; Make about 20% of the points be "bad" data.

bad = WHERE(RANDOMU(1, nx, ny) gt 0.8)

Z[bad] = 999

 

; Create the mask of the bad data points.

mask = REPLICATE(1, nx, ny)

mask[bad] = 0

 

;***** Fit the function *****

yfit = GAUSS2DFIT(Z, B, /TILT, MASK=mask)

 

; Report results:

PRINT, 'Should be: ', STRING(A, FORMAT='(6f10.4)')

PRINT, 'Is: ', STRING(B, FORMAT='(6f10.4)')

 

; Create an array with our fitted results

xprime = (X - B[4])*cos(B[6]) - (Y - B[5])*sin(B[6])

yprime = (X - B[4])*sin(B[6]) + (Y - B[5])*cos(B[6])

Ufit = (xprime/B[2])^2 + (yprime/B[3])^2

Zfit = B[0] + B[1] * EXP(-Ufit/2)

 

; Plot the results. The black dots are missing data.

im = IMAGE(BYTSCL(Z, MAX=20), RGB_TABLE=40, LAYOUT=[1,1,1])

 

; Contour plot of the fit.

c = CONTOUR(Zfit, /OVERPLOT, C_THICK=[4], COLOR='red')


Version History

4.0.1

Introduced

8.0

Added FITA and MASK keywords

See Also

COMFIT, GAUSSFIT, POLY_FIT, REGRESS, SFIT, SVDFIT, “Curve and Surface Fitting” (Chapter 7, Using IDL) in the help/pdf directory of the IDL installation.