Program: gaufit2d
Purpose: Find position or major-, minor axis and position angle of
sources by fitting 2 dim. gauss function parameters to
data in a box..
Category: MODELS, FITTING
File: gaufit2d.c
Author: M.G.R. Vogelaar (GUI: J.P. Terlouw)
Description: CONTENTS:
-PURPOSE
-GRAPHICAL USER INTERFACE
-DESCRIPTION OF THE MODEL
-FIT PROCEDURE
-WARNING GRID CONVERSIONS
-LIST OF KEYWORDS/BUTTONS AND INPUT FIELDS
-APPENDIX I
-EXAMPLE OF OUTPUT IN LOG FILE
PURPOSE:
Program GAUFIT2d is a tool for finding the exact position of
the maximum in a source or to find major & minor axes and
position angle of a source by fitting a two dimensional gauss
function to the data. Also the volume under the 2d-gauss is
returned in units of the amplitude (see appendix 2).
Output is a list of fitted parameters on screen or in the
GIPSY Log file. The program works either in a coordinate system
of grids or in a system of seconds of arc.
If you press the button, an example run is shown.
GRAPHICAL USER INTERFACE:
The graphical interface is closely related to Hermes, i.e.
you are still able to run the program from within Hermes
and specify keywords on the command line. The GUI is a
nicer way to present keywords and an easier way to enter
values for these keywords. Keyword values entered on the
command line will have an effect on the status of the GUI
immediately. As an example start the program with something
like this:
gaufit2d inset=test f 0 box=-192 337 -182 346
Then the GUI starts with filled input fields for SET and BOX.
Note that for input fields the same input rules apply as for
keywords in Hermes
DESCRIPTION OF THE MODEL:
The function that this program uses to find model parameters
is:
F(x,y) = A * EXP[ -4*Ln(2)*( (xr/HW1)^2+(yr/HW2)^2 ) + Z0
where: xr = -x1 * sin(phi) - y1 * cos(phi)
yr = x1 * cos(phi) - y1 * sin(phi)
and: x1 = x - x0
y1 = y - y0
The parameters are fitted either in a grid system
or an system in which both axis units are arcsec.
The parameters are:
A, amplitude in units of the image data.
x0, center in X in grids.
y0, center in Y in grids.
HW1, first Full Width at Half Maximum in arcsec or grids
HW1, second Full Width at Half Maximum in arcsec or grids
phi, angle between X axis and first axis of ellipse,
angle always between -45 to 45 degrees.
Note that angles in the grid system can differ from
angles in the arcsec system if the length of a
grid in X is unequal to the length in Y.
Z0 Offset of the data (zero level) in units of the image
data.
FIT PROCEDURE:
The method used to fit these parameters is a least squares
method. This method needs initial estimates. The fit can be
very sensitive to these estimates, therefore you are able
to overrule the default estimates set by the program and
you are also able to fix one or more initial estimates to
decrease the degrees of freedom and to improve the results.
Estimates are calculated with first and second moments of
the data. The parameters of a 2-dim gauss function can be
expressed in these moments. For the calculations of
initial estimates, a cutoff value is introduced. This
value includes only data if the image value of a pixel
is more than 5% of the maximum value in a given box.
This maximum is in fact the average of the four pixels
with the highest value.
The least squares fit stops when successive iterations fail
to produce a decrement in reduced chi-squared less than
a certain value called the relative tolerance. Maximum
accuracy can be obtained by setting this value to 0 which
is the default in this program. If a fit needs more iterations
than the maximum number, you can set the tolerance to e.g.
0.05 (or increase the maximum number of iterations).
GAUFIT2D can do its calculations in either a (square)
grid system or a sky system where both units are seconds
of arc (from now on called the arcsec system).
In the arcsec system a position is calculated by multi-
plying the grid position with the grid spacing (which
can be different in the x- and y direction) and both
estimates- and fit routines do their calculations
in the arcsec system. If the header information from a
set allows grids converted to arcsec then the calculations
are default in the arcsec system. This can be overruled
by pressing the button (appears in: EXTRA KEYWORDS).
WARNING GRID CONVERSIONS:
The fit results (an ellipse with major, minor axis and
position angle) are plotted if the button is on.
The image data is always displayed in grids and therefore
the fit is also in grids. This can result
in confusing results because if the grid spacings in
both x and y are not equal, the position angle is
different in the arcsec and grid system.
LIST OF KEYWORDS/BUTTONS AND INPUT FIELDS:
Keywords: Each keyword corresponds to a button or input field!
EXAMPLE= (Button: )
Run an example. It creates a simple set called
"GAUFIT2Dexampleset". It fits parameters and shows the result
in GIDS.
INSET= (Input field: SET)
Input set and optionally subsets. The subset must be
two dimensional.
Model parameters are converted to arcsec if the input is a
spatial map, i.e. the first axis is a spatial longitude axis
and the second is a spatial latitude axis. If you want the
results in grids in any case use the button under
'Extra keywords'.
Maximum number of subsets is 1024. You can change the current
subset level by pressing one of the '<' or '>' buttons or
specify a new level in the 'subset' input field.
Examples: INSET=AURORA FREQ 10:20
INSET=AURORA FREQ
INSET=AURORA FREQ 8 10 15:18
BOX= (Input field: BOX)
Enter the limits in grid coordinates of
the box from which you want to extract data for the fit.
Examples: BOX=-20 -10 20 10
BOX=0 0 D 10 10 center at 0,0 size: 10x10
SUBSET= (Input field: SUBSET)
NOTE: This field can only be used AFTER you specified
subsets in the SET input field. Then you can change
the value in this field to force reading data from another
subset. Subsets can also be changed with the ">" and "<"
buttons to the right of the input field.
AMPLITUDE= (Input field: AMPLITUDE)
Estimate for the value of the central value in least
squares fit. Units are the same as the units of your data.
This parameter can be adjusted to be a better estimate
in your opinion or can be fixed with the /
button. A fixed value will always re-appear in the column
with fitted parameters.
MAJOR= (Input field: FWHM X)
Estimate for the value of the major axis in grids or in
arcsec. The value is in arcsec if a conversion from
grids to arcsec is possible and the button is
off. Otherwise the input must be in grids.
MINOR= (Input field: FWHM Y)
Estimate for the value of the minor axis in grids or
in arcsec.
CENTERX0= (Input field: CENTER X0)
Estimate in grids for the X-position of the maximum of
the source.
CENTERY0= (Input field: CENTER Y0)
Estimate in grids for the Y-position of the maximum of the
source.
ANGLE= (Input field: ANGLE)
Call the axes of an unrotated frame the X and Y axis,
and the main axes of a rotated ellipse X' and Y'. Then
ANGLE= is an estimate for the smallest angle (positive
or negative) between the X and X' (or Y and Y') axis.
The range is always from -45 deg to 45 deg.
For the arcsec system this angle is converted to a
real position angle, i.e. an angle between the major
axis and the direction of the North.
ZEROLEVEL= (Input field: ZERO LEVEL)
Estimate for the offset of the
two dimensional gauss. In some cases (e.g. antenna patterns)
this parameter can be set to zero and fixed.
All these parameter keywords also have prefixes F_ for the
corresponding fixed/free toggle used in the lsq.fit and
R_ to restore the previously calculated moments estimate.
QUIT= (Button: )
Terminate program.
FIT= (Button: )
Start a least squares fit using the estimates
that were defined before pressing this button.
NEXT= (Button: >)
Go to the next subset in the sequence.
Works only if more than one subset was entered.
PREV= (Button: <)
Go to the previous subset in the sequence.
LOG= (Button: )
Send all fit information to the Hermes Log and
Screen file. The button is only active if there is
a fit result.
Some buttons/input fields appear when you press the
button:
TOL= (Input field: )
Fractional tolerance for the chi square (chi2) default
set to 0.0.
Fitting is stopped when successive iterations fail
to produce a decrement in reduced chi-squared less
than TOL= The value cannot be less than a
certain minimum as set by the system. This means that
maximum accuracy can be obtained with TOL=0.0
LAB= (Input field: )
Value for mixing parameter default set to
0.001. Mixing parameter in the least squares fit function.
LAB= determines the initial weight of steepest descent
method relative to the Taylor method. LAB= should be
a small value (e.g. 0.001).
CLIP= (Input field: CLIP)
Value below which data is excluded in
the least squares fit. Default the clip is not set, meaning
that all data contributes to the fit. If you image data is
very noisy then a better fit could be made if you set a
clip near the value of the zero level that you expect.
MAXITS= (Input field: MAXITS)
Maximum number of iterations allowed in the least squares
fit routine. The default is 100.
GRIDS= (Button: )
If On: All calculations are in grids.
If off: If header information allows conversion from
grids to arcsec, then calculations are in arcsec.
Display:
VIEW= (Button: )
If necessary start GIDS. Display current
subset and ellipse representing the major,minor axis
with position angle.
ERASE= (Button: )
Erase current plotted ellipse in GIDS and restore image.
RESIDUE= (Button: )
Calculate difference between data and
model and display the result in GIDS. The button is only
active if a set is viewed and a fit is made.
The residue is stored in a set on disk. The name of the
set is "XXXXresidualXXXX". The set is removed after the
program has finished. Note that the program only
shows the residual in the current box.
Notes: APPENDIX 1:
A note about chi2. The quantity chi-squared is defined as:
chi2 = Sum[ (Y_i-Ymodel_i)^2/sigma_i^2 ]
and the reduced chi-squared is:
reduced-chi2 = chi2 / (N-n-1)
where N is the number of data pints in the fit and n the
number of parameters.
However, the values of sigma_i (sigma per pixel) is not
known and set to 1. Then the reduced chi2 is equal to
the variance of the fit (Bevington, Data Reduction and
Error Analysis for the Physical Sciences, p 188) and
chi2 is not a measure of the goodness of fit anymore.
So we prefer to list only the variance of the fit.
If you assume the fit is a good approximation of the model
then reduced-chi2 ~ 1 and because
reduced-chi2 = variance / sigma_av
you have an approximate value for sigma_av, the
weighted average of the individual variances.
APPENDIX 2:
F(x,y) = A * EXP[ -4*Ln(2)*( (xr/HW1)^2+(yr/HW2)^2 ) + Z0
The volume V under a 2d-gauss limited by the offset ZO
is given by: V = 2.A.HW1.HW2.PI.Erf(Inf).Erf(Inf)/Ln(256)
= 1.13309.A.HW1.HW2
This volume is expressed in units of the amplitude, therefore
it will be the same number while working in grids or
physical coordinates.
Example: EXAMPLE OF OUTPUT IN LOG FILE:
===================================== GAUFIT2D ================================
Date of fit : 03-Mar-1998 (17:06:38)
Name of set : [test]
Fit box : [-192 337 -182 346]
Iterations : 8
Variance of fit : 0.192932
Points in box : 110
# excluded by clip : 0
Number of blanks : 0
Volume under gauss : 568.882 in units of the amplitude
Coordinate system : Arcsec
Size of 1 grid : 4.32 x 5.38978 arcsec
Center X: -187.34 (grids) = 179.1631786 (DEGREE) = 11h56m 39.1s
Center Y: 341.41 (grids) = 54.1635537 (DEGREE) = 54d 9m48.79s
Major axis: 19.54'', Minor axis: 15.66'', Position angle: 155.3 (deg)
===============================================================================
Updates: Feb 26, 1998: VOG, Document created.
Jul 13, 2005: VOG, Repaired bug in QUIT= keyword handler
(avoid use of freefmatrix with null pointer)
Apr 14, 2009: VOG, Removed unused macro for NINT()
Apr 17, 2010: VOG, Extra fit information in LOG