Program: FINDGAUSS
Purpose: Fitting and subtraction of 2-D gaussians
Category: FITTING, MODELS, TABLES
File: findgauss.c
Author: M. Vogelaar
Keywords:
INSET= Give set, subsets:
Maximum number of subsets is 2048.
BOX= Give box in ..... [entire subset]
** CUTOFF= Give cutoff ratio for gaussian [1/100]
This is the ratio of the gaussian function values
at central point and a point x on the boundary.
OUTSET= Give output set (, subsets):
Output set and subset(s) for the result. The number of
output subsets is the same as the number of input sub-
sets.
SUBFILE= Use parameter file to subtract gaussians? Y/[N]
Instead of calculating gaussians, it is possible to
use parameters from a previous session to subtract
the fitted gaussians. The parameters are written to
an ASCII file and can be edited before the actual sub-
traction.
If SUBFILE=Y:
PARAMFILE= Name of file with stored parameters: [Stop program]
File must be created in a previous run of this program.
Only a subtraction on INSET= is performed.
If SUBFILE=N:
PARAMFILE= Name of file to store parameters: [No file]
Write gauss parameters to file so that these parameters
can be edited and used for subtraction.
If a given file name already exists, APPEND= must be
specified.
APPEND= File exists, ok to append? [Y]/N
A specified file name already exists.
You can append to this file with APPEND=Y. If APPEND=N
you will be prompted for another filename.
Pre specification of initial parameters in the fit:
AMPLITUDE= Amplitude (map units): [Calculated]
BEAM= Beam: major, minor: [Calculated]
If indicated, in seconds of arc.
X0Y0= Centre in grids: [Calculated]
ANGLE= Angle of major axis wrt pos. Y-axis (deg): [Calculated]
ZERO= Zero level (map units): [Calculated]
MASK= Mask (0=fixed,1=free): [1,1,1,1,1,1,1]
FITSIZE= Give size (x*y) of fit box: [Calculated]
Give a fixed size for the 'fitbox' instead of
determination of suitable sizes by the program.
SIGMA= Rms of noise in map: [None]
VALRANGE= Range of values in search area: [4*sigma, ->]
Or if SIGMA= was not specified:
Range of values in search area: [All data]
Data values outside given range cannot be candidates
for a gauss fit.
** TOLERANCE= Convergence criterion. [0.01]
Fitting stops if successive iterations fail to
produce a decrement in reduced chi-squared less than
TOLERANCE=
** LAB= Mixing parameter: [0.01]
Mixing parameter, LAB determines the initial
weight of steepest descent method relative to the Taylor
method. LAB should be a small value (i.e. 0.01).
** CLIPLOHI= Give clip levels: [Include all data]
Examples:
1) CLIPLOHI= include all data
2) CLIPLOHI=3 include all data >= 3
3) CLIPLOHI=3 5 include all data >=3 and <= 5
Subtract only gaussians with amplitudes and beams in
range AMPRANGE= and BEAMRANGE=:
AMPRANGE= Amplitude range to subtract: [All amplitudes]
BEAMRANGE= Beam range to subtract: [All values]
TABNAME= Name of the table to store table data: [GAUSS]
Description: The program FINDGAUSS fits 2 dimensional gaussians to
sources in an image. The fitting is used for example
to get rid of point sources in a map or to find positi-
on and flux of these point sources. The set that is
examined is INSET=. The defined subset must be two
dimensional. The number of subsets can be greater than
1 but cannot exceed 2048. The items corresponding with
the grid spacings (CDELTn) must be available from the
the header of the input set, otherwise the program will
abort. If it can find the grid spacing, and the units
of the spacing is in degrees (along both subset axes),
then input will be done in seconds of arc else it will
be done in units as found in the header.
With BOX= the area that will be searched for gaussians
can be made smaller. The default is the entire (sub)set.
OUTSET= defines the name of the set where the result is
stored. This result is always the subtracted INSET= and
therefore can be used as an inspection set. If you want
to write to an already existing set, you will be warned
with the message: Will overwrite data, okay ? [Y] and
at the OKAY= prompt.
The application tries to allocate memory for your image
and when it is not able to allocate enough memory, it
aborts. Then you can try again with a smaller box.
If the result of the subtractions is not exactly what
you wanted or expected because a couple of unwanted
fits were extracted, there is no need to rerun the
application completely. It is possible to write the
gauss parameters to an ASCII file with PARAMFILE=
This file can be edited. To use it for subtraction,
you must give SUBFILE=Y first. Then you are prompted
with PARAMFILE= Pressing carriage return will stop the
program. If a file name is given and the file can be
opened, the parameters in this file are used to define
the gaussians. The parameter file is made in a run
of FINDGAUSS with SUBFILE=N For each subset a file
name can be specified. For each fit, the file is
extended with one line with 12 items:
1) An integer 0 or 1 (0: not selected for subtraction)
(1: selected for subtraction)
2) Value of fitted amplitude in map units.
3) Major axis in (if possible) seconds of arc.
4) Minor axis in (if possible) seconds of arc.
5) Centre X coordinate in grids.
6) Centre Y coordinate in grids.
7) Angle between major axis and positive y-axis in
degrees (In most cases the pos. y-direction is
also the direction of the North).
8) Zero level in map units.
9) Lower X value of (fit)box where gauss is fitted.
10) Lower Y value.
11) Upper X value.
12) Upper Y value.
The fitting:
FINDGAUSS finds the maximum in a subset by examining
data with values above a certain level or between two
levels (VALRANGE=) If SIGMA= has a value, the default
for VALRANGE= is one (lower) clip equal to 4 times the
value in SIGMA= If SIGMA= was not specified, all data
in a subset is used, but this will result in attempts to
fit gaussians with an amplitude comparable to the noise,
so it is better to define a threshold. One value, say x
will result in a value range starting with x and
including all values greater than x. Two values, say x,y
indicate a range where the search area is determined by
all values between x and y (x, y values included).
As mentioned, the fit parameters are:
1) Amplitude,
2) Major axis,
3) Minor axis,
4) Centre X,
5) Centre Y,
6) Position angle of major axis,
7) Zero level.
In the default situation, all parameters are free in the
fit (i.e. after the fit they can have other values than
the pre specified or estimated values).
You can fix parameters in the fit with the keyword
MASK= The default is MASK=1 1 1 1 1 1 1 (All free). If
MASK=1 0 0 1 1 1 1 you fix the major and minor axis etc.
Fixed parameters must have an initial value. This can be
a value estimated by the program, or it can be a pre
specified value. You can give pre specified values for
fixed and free parameters (which will in both cases
overrule the calculated estimates). Pre specifications
are given with the keywords:
AMPLITUDE= amplitude in map units,
ANGLE= position angle in degrees (angle between
major axis and positive Y-axis counted
counter clockwise). You will be warned
if the pos. Y-axis does not align with the
north in your map.
BEAM= major, minor axes (in units as prompted).
X0Y0= centre X, Y in grids.
ZERO= zero level in map units.
The fitting can be influenced by two variables given in
TOLERANCE= and LAB= TOLERANCE= is the relative tolerance.
Fitting of the gaussian function stops when
successive iterations fail to produce a decrement in
reduced chi-squared less than TOLERANCE=. If its value
is less than the minimum tolerance possible, it will be
set to this value. This means that maximum accuracy can
be obtained by setting TOLERANCE=0.0.
LAB= is the mixing parameter, LAB= determines the initial
weight of steepest descent method relative to the Taylor
method. LAB= should be a small value (i.e. 0.01).
The data that is used in a fit can be limited with the
keyword CLIPLOHI= If you give no values (default) than
all data in the 'fitbox' is used in the fit. If you give
one value, then all data greater than or equal to that
value will be included and if you give two values, then
data with values greater than (or equal to) the first
clip and less than (or equal to) the second clip
are included in the fit.
The decision to subtract or not is limited by the key-
words AMPRANGE= and BEAMRANGE= Default is the subtraction
of all fitted gaussians. If you give one value for a
range keyword, then all amplitudes (or beams) with a
value greater than this value will be included in the
subtraction. If you give two values, then all
amplitudes (or beams) with values greater than (or equal
to) the first range value and less than (or equal to)
the second value will be included in the subtraction.
Finding a 'fitbox':
If the maximum in a subset is found, a box is created
by appending a 'border' around this pixel. This box
has size 3x3. Then the box is extended with a new
surrounding border. If the maximum in the new extension
exceeds the maximum in the previous extension, the
first box becomes the 'fitbox'. This is also true if
the mean in the new extension exceeds the mean in the
previous extension or if the decrease in the mean is
less than SIGMA= (i.e. if SIGMA= was specified).
The last criterion is that the mean becomes smaller
than the maximum times a cutoff factor given in
CUTOFF= The default is 1/100. You have to experiment
with the keywords SIGMA= and CUTOFF= to get an idea
how the growth of a fitbox is stopped.
Fixing parameters:
You have the possibility to combine the fixing of
parameters with the pre specification of them. Good
results can be obtained by pre specifying the
beam (BEAM=) and angle (ANGLE=) but with all parameters
free in the fit (MASK=1 1 1 1 1 1 1)
If, for instance, you fix beam and angle, then more
gaussians will be found, but after the subtraction you
will also find more residuals because not all sources
have the pre specified properties, but just gave a 'good'
fit.
Use of tables:
The results of this application are stored in a table.
The default name of the table is GAUSS and the columns
are AMPLITUDE, MINOR, MAJOR, X0grid, Y0grid, X0phys,
Y0phys (physical counterparts of X0grid, Y0grid), ANGLE,
ZEROLEV, AMPerr, MAJerr, MINerr, X0err, Y0err, ANGerr,
ZEROerr, XLO, YLO, XHI, YHI (4 box coordinates in grids).
The name of the table can be changed, the name of the
columns not. If one table name is used twice, the
contents of the last table will overwrite the first one.
Notes: .......
Example: Example of a default file:
AMPLITUDE= / Program calculates estimate
AMPRANGE= 0.04 12 / Use in subtraction only fits
with ampl. in this range
ANGLE= 0 / Pre specify angle, 0 degrees
wrt north (to east)
BEAM= 19 13 / Pre specify beam, major axis
axis 19 arcsec, minor 13 arcsec.
BEAMRANGE= 5 40 / Use in subtraction only fits
with beams in this range
BOX= / Operate on entire subset
CUTOFF= 1/100 / Limit 'fitbox' to boundaries
with mean > cutoff*maxampl.
PARAMFILE=storeparm.txt / Sore fit parameters in ASCII
file
FITSIZE= / Do not use a fixed size for
the fitbox.
INSET=AURORA freq 0 / Set and subset
MASK= 1 1 1 1 1 1 1 / All parameters free in fit
OUTSET=SUBAURORA / Name of subtracted data
SIGMA= 0.01 / Rms of noise
SUBFILE=N / Do not use an ASCII file
with gauss parameters for
subtraction
TABNAME= / Use default name for table
VALRANGE= / Sigma was specified so
default is VALRANGE=4*0.01
which implies that only
data > 0.04 can be candidates
for a fit
X0Y0= / Program calculates estimate
ZERO= / Program calculates estimate
Updates: Sep 22, 1992: VOG, Document created.
Oct 4, 1993: VOG, Repaired bug in preset of FWHM
Mar,25 1994: VOG, Replaced some 'tofchar's by str2char
in routine 'putintable'.