The GAUSS2DFIT function fits a two-dimensional, elliptical Gaussian equation to rectilinearly gridded data.
And the elliptical function is:
The parameters of the ellipse U are:
The rotated coordinate system is defined as:
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.
Coefficients of the computed fit are returned in argument A .
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 computationally-intensive routine takes approximately 4 seconds for a 128 by 128-element array on a Sun SPARC LX. 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.
The dependent variable. Z should be a two-dimensional array with dimensions ( N _{ x} , N _{ y} ). Gridding in the array must be rectilinear.
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:
This example creates a 2D gaussian, adds random noise and then applies GAUSS2DFIT.
nx = 128 & ny = 100 ; Define array dimensions.
A = [ 5., 10., nx/6., ny/10., nx/2., .6*ny]
;
Define input function parameters.
X = FINDGEN(nx) # REPLICATE(1.0, ny)
Y = REPLICATE(1.0, nx) # FINDGEN(ny)
;
Create X and Y arrays.
U = ((X-A[4])/A[2])^2 + ((Y-A[5])/A[3])^2
;
Create an ellipse.
Z = A[0] + A[1] * EXP(-U/2) ; Create gaussian Z.
Z = Z + RANDOMN(seed, nx, ny) ; Add random noise, SD = 1.
yfit = GAUSS2DFIT(Z, B) ; Fit the function, no rotation.
PRINT, 'Should be: ', STRING(A, FORMAT='(6f10.4)')
PRINT, 'Is: ', STRING(B(0:5), FORMAT='(6f10.4)')
;
Report results.