SCULIB_FIT_2D_GAUSSIAN

←Prev
SURF Programming Interface
Next→
TOC ↑

SCULIB_FIT_2D_GAUSSIAN

fit a 2D Gaussian to an image

Description:

This routine fits a 2D Gaussian to an image of a source. It takes as input data a map made of the source with pixels on a square grid.

If status is good on entry, the routine will calculate the 0th ( $U_{0}$ ) and 1st order ( $U_{x}$ , $U_{y}$ ) moments of the image on the map. Map pixels with bad quality are ignored:

The 0th order is calculated as

U_{0} = \sum_{i j} f_{i j},

(27)

where i,j are the pixel indices. The 1st order moment in X is calculated as

U_{x} = \sum_{i j} x f_{i j},

(28)

where x is the x position of pixel i,j; and the 1st order moment in Y is calculated as

U_{y} = \sum_{i j} y f_{i j},

(29)

where y is the y position of pixel i,j.

If the 0th order moment of the image was 0, i.e. there is no image on the map, then a warning message will be output and all image parameters set to 0 or bad quality. Otherwise, the x,y centre of the image is calculated from:-

X_{}

cen = U_x

U_{0}

(30)

Y_{}

cen = U_y

U_{0}

(31)

With this information the 2nd order moments can be calculated about the centre of the image:-

U_{x x} = \sum_{i j} {(x - X_{c e n})}^{2} f_{i j}

(32)

U_{y y} = \sum_{i j} {(y - Y_{c e n})}^{2} f_{i j}

(33)

If the image had a 2D Gaussian profile,

f = A exp (- a_{x}^{2} x^{2}) exp (- a_{y}^{2} y^{2}),

(34)

then it can be shown that for x:-

a_{x}^{2} = \frac{U_{0}}{2 U_{x x}}

(35)

and the half-width at half-maximum is calculated from:-

X_{}

HWHM =

∘ ln (2)ax(36)

Lastly, the volume under the image and the variance on it are calculated. These will be the sum of all map pixels and their variances that lie within a radius of 2 $*$ max(HWHM_x, HWHM_y) from X_CENTRE, Y_CENTRE. If the HWHM in either axis is zero then a warning message will be output and PEAK_QUAL set bad. Otherwise, PEAK will be set to volume / ( $π$ $*$ HWHM_x $*$ HWHM_y).

Invocation

CALL SCULIB_FIT_2D_GAUSSIAN (IDIM, JDIM, MAP_DATA, MAP_VARIANCE, MAP_QUALITY, X, Y, PEAK, PEAK_VAR, PEAK_QUAL, X_CENTRE, Y_CENTRE, X_HWHM, Y_HWHM, STATUS)

Arguments

IDIM = INTEGER (Given)

‘I’ dimension of map

JDIM = INTEGER (Given)

‘J’ dimension of map

MAP_DATA (IDIM, JDIM) = REAL (Given)

map data values

MAP_VARIANCE (IDIM, JDIM) = REAL (Given)

map variance

MAP_QUALITY (IDIM, JDIM) = REAL (Given)

map quality

X (IDIM) = REAL (Given)

x axis of map

Y (JDIM) = REAL (Given)

y axis of map

PEAK = REAL (Returned)

peak height of Gaussian

PEAK_VAR = REAL (Returned)

variance on PEAK

PEAK_QUAL = INTEGER (Returned)

quality of PEAK

X_CENTRE = REAL (Returned)

x of image centre

Y_CENTRE = REAL (Returned)

y of image centre

X_HWHM = REAL (Returned)

half-width half-max of image in x-axis

Y_HWHM = REAL (Returned)

half-width half-max of image in y-axis

STATUS = INTEGER (Given and returned)

global status

Copyright

←Prev
SURF Programming Interface
Next→
TOC ↑