AIPS HELP file for SAD in 31DEC09
As of Sat Nov 21 20:28:50 2009
SAD: Task to fit Gaussian models to an image by least-squares
INPUTS
INNAME Input image name (name)
INCLASS Input image name (class)
INSEQ 0.0 9999.0 Input image name (seq. #)
INDISK 0.0 9.0 Input image disk drive #
INVERS -1.0 46655.0 Input image MF ext. version
<= 0 => new
BLC 0.0 4096.0 Bottom Left corner of fit
TRC 0.0 4096.0 Top Right corner of fit
DORESID -1.0 2.0 >0 => 1 Catalog residual map
OUTNAME Output image name (name)
OUTCLASS Output image name (class)
OUTSEQ -1.0 9999.0 Output image name (seq. #)
OUTDISK 0.0 9.0 Output image disk drive #
NGAUSS 0.0 40000.0 Max. Number of components
CPARM 0.0 9999.0 Flux cutoff levels in
descending order
ICUT 0.0 9999.0 Retry level
BWSMEAR -0.1 0.1 Bandwidth smearing corr.
SORT Sort Order of output ' '=RA
DOCRT -3.0 132.0 >0 => print fits on CRT
> 72 => terminal width
=0 => print suppressed
<0 => line printer
FITOUT
Disk file to save fit info
OUTVERS -1.0 46655.0 CC file version # for
solutions: -1 => none, 0 new
DOALL -1.0 1.0 >0 => fit multiple peaks
DOWIDTH -1.0 2.0 >0 => fit widths
GAIN 0.0 1.0 Amp-dependent part of retry
and warning levels
DPARM Reject components with
(1) peak < D(1)
(2) flux < D(2)
(3) rms > D(3) + GAIN * flux
(in quadrature)
(4) widths > D(4) cells
(5) peaks > D(5) cells
outside interior area
(6) peaks > D(6) cells
outside image
(7) total residual flux in
the fitting box
[no default]
(9) Use IN2NAME etc to set
image of rms
>= 2 => CPARM is in S/N
(see HELP)
(10) Expand fit box pixels
IN2NAME RMS image name (name)
IN2CLASS RMS image name (class)
IN2SEQ 0.0 9999.0 RMS image name (seq. #)
IN2DISK 0.0 9.0 RMS image disk drive #
PRTLEV 0.0 1023.0 Debug print level
PBPARM Primary beam parameters:
(1) level to believe - <= 0
means do not apply a primary
beam (2) > 0 use (3)-(7)
HELP SECTION
SAD
Type: Task
Use: SAD (Search And Destroy) is a task to find a number of potential
sources in a (sub)image, fit Gaussian components to these, and
optionally print the results, store them in a CC file, and create
a residual file. The searching may be restarted with different
parameters to remove additional sources. Only one plane of a cube
is fit at one time. However, the MF file keeps track of the plane
number and may contain data for more than one plane. With
additional tasks (to be written yet), one should be able to use
this feature to do spectral fitting of the objects found spatially
by SAD.
The display produced by SAD shows the results of the fitting,
estimates of the uncertainties in the answers (based on the actual
noise and theory rather than the quality of the fitting), and
attempts to deconvolve the Clean beam from the fit sources. When
the residual image near a source exceeds a limit controlled by the
cutoff (CPARM) and GAIN parameters, an H, L, or S will be shown
near the fit peak value. H is for a high point in the residual,
L for a rather negative point, and S for a large rms. If the
width of the display allows, the max residual is displayed with
the deconvolution results and an * is placed near excessive
residuals.
Adverbs:
INNAME......First image name (name). Standard defaults.
INCLASS.....First image name (class). Standard defaults.
INSEQ.......First image name (seq. #). 0 => highest.
INDISK......Disk drive # for the first image. 0 => any.
INVERS......Version number of a MF table extension file attached to
the input image and copied to the ouput residual image (if
any). <= 0 => create a new one, > 0 => use a pre-existing
one and subtract the components in it before searching for
more.
BLC.........Bottom left corner of area of image to fit.
TRC.........Top right corner of area of image to fit.
DORESID.....>0 -> Create a residual map of subimage
OUTNAME.....Residual map name. Standard defaults.
OUTCLASS....Residual map class. Standard defaults.
OUTSEQ......Residual map seq. #. 0 => highest unique.
OUTDISK.....Residual map disk no. 0 => highest with room.
NGAUSS......The maximum number of components to search for
0->10. Maximum number is 40000.
CPARM.......Search for potential sources down to these levels in
units of the input image or in units of S/N (when
DPARM(9) >= 2.0 and IN2NAME et al are specified). No
defaults - list in descending order, A value of 0 => no
more levels. Thus CPARM = 8, 4, 0.666, 0 => 3 rounds
of searching. It has been found useful to fit the
higest levels first, then look somewhat lower, and then
even lower. Low level signal around a strong source can
sometimes confuse the fitting. CPARM all 0 -> 3,0 in
S/N or 3*rms,0 in image units.
ICUT........If the peak residual in an island, after fitting a single
Gaussian exceeds ICUT in image or S/N units (as CPARM)
adjusted in quadrature by GAIN*flux, then the program
will attempt a 2-Gaussian fit to the island. 0 =>
CPARM(i) for each search level i.
BWSMEAR.....If not 0, the Clean beam will be smeared by a Gaussian
in the radial direction of FWHM proportional to the
radius from the pointing position times abs(BWSMEAR).
Set it roughly to the channel bandwidth divided by the
center frequency. The Clean beam is used as an initial
estimate of the source size and is used in the
deconvolution attempt to find the true size from the fit
size. The peak intensity printed will be corrected for
this effect; set BWSMEAR < 0 to suppress this
correction.
SORT........Sort order of output.
' ' -> 'R'-> Right Ascension
'D' -> Declination
'X' -> X pixel (~RA if no rotation)
'Y' -> Y pixel (~Dec if no rotation)
'S' -> Flux
DOCRT.......False (< 0) use the line printer if FITOUT = ' '
else write named FITOUT file only.
-2 causes the line feed characters in FITOUT to
be suppressed
-3 causes the documentation lines to be suppressed
Null (= 0) print suppressed.
True (> 0) use the terminal interactively. The task will
use the actual terminal width as a display limit
unless 72 < DOCRT < width. In that case, the display
limit will be DOCRT characters.
FITOUT......Disk file name in which to save the fit (line printer)
output. ' ' => use scratch and print immediately for
interactive jobs - batch jobs use FITOUT =
'PRTFIL:BATCHjjj.nnn' (jjj = job #, nnn = user #). When
FITOUT is not blank, multiple outputs are concatenated
and the file is not actually printed.
OUTVERS.....= 0 => Create a new CC file on input image.
< 0 => do not write any CC file
> 0 => append components to pre-existing CC file
SAD writes the deconvolved widths (0 on failure of
deconvolution) into the CC files, if the input image is a
Clean image. If DOWIDTH <= 0, then it writes point
components into the CC file.
DOALL.......>0 -> If an island has multiple (up to 4) peaks, fit
multiple gaussians. Otherwise fit 1 gaussian to each
island.
DOWIDTH.....<= 0 -> Fit all sources with Clean beam
>0 & <=1 -> Let fit find size of source starting with the
Clean beam
>1 -> Let program guess the sizes and then fit them
(this can be somewhat unstable)
GAIN........Retry, cutoff and warning levels are computed by
SQRT ( level**2 + (GAIN*flux)**2 )
where the level is the no-signal retry (ICUT), cutoff
(DPARM(2)) or warning (CPARM cutoff) level and flux is
the peak of the component in question. No default.
DPARM.......Components are fit within islands that are rectangular in
shape. It may be possible to get bad solutions, so the
program will reject components with:
(1) peak <= DPARM(1). No default; 0 is okay if you
mean it. Units are image units unless DPARM(9)>=2,
in which case they are in S/N.
(2) flux <= DPARM(2). No default; 0 okay
(3) rms in the box > sqrt (DPARM(3)**2 + (GAIN*flux)**2).
0 => 1000. Units are image units unless DPARM(9)>=2,
in which case they are in S/N.
(4) widths > DPARM(4) cells. 0 => 1000.
(5) peaks > DPARM(5) cells outside the island (fitting
area) but interior to the image. No default;
< 0 => inside the fitting area required.
(6) peaks > DPARM(6) cells outside image. No default.
This allows SAD to fit components legitimately on the
edge of the whole image while rejecting spurious ones
extrapolated to positions outside the fitting
rectangle but fully in the image.
(7) total residual flux in the fitting box >
SQRT (DPARM(7) + (GAIN*flux)**2). No default.
*** (9) > 0 => use IN2NAME, IN2CLASS, etc as an image of the
RMS so that it may vary with position. = 1 => use
the image only to evaluate the uncertainties in the
fit components. >= 2 => use the image also to
evaluate the noise at each pixel as data are read.
In that case CPARM, ICUT, DPARM(1), and DPARM(3) are
in S/N units.
*** (10) Expand (or contract) the box used to fit the
component.
IN2NAME.....RMS image name (name). NO DEFAULT.
The task RMSD may be used to prepare such an image.
Note that the RMS image must be at the same coordinates
more or less as the INNAME image but they do not have to
allign. The RMS is found by looking up the coordinate
of the coponent in the RMS image.
If DPARM(9) >= 2, the RMS image is ASSUMED to have the
same size and coordinates as the input image. Bad
things will happen, perhaps not in obvious ways, if this
is not true. In this case, the noise image is used
throughout rather than simply in evaluating the
uncertainties in the component parameters.
IN2CLASS....RMS image name (class). Standard defaults.
IN2SEQ......RMS image name (seq. #). 0 => highest.
IN2DISK.....RMS drive # for the first image. 0 => any.
PRTLEV......To find out which islands are rejected by the above for
what reasons, set PTRLEV to the sum of the desired
reasons: 1 peak, 2 flux, 4 rms, 8 widths, 16 X position
outside interior window, 32 Y position outside interior
window, 64 X position outside image, 128 Y position outside
image., and 256 for residual flux. 511 does all of them.
If you don't want to see the low peak ones, set PRTLEV
to an even number. Any that are rejected for low peak
will not be displayed even if they are also rejected for
a selected reason.
PBPARM......Primary beam parameters:
(1) Lowest beam value to believe: 0 -> do not do the
primary beam correction. This correction is done
to the printed parameters only. The beam value used
is max (PBPARM(1), that computed from (2)-(7)).
(2) > 0 => Use beam parameters from PBPARM(3)-PBPARM(7)
Otherwise use default parameters for the VLA (or
ATCA where appropriate)
(3-7)..For all wavelengths, the beam is described by the
function:
1.0 + X*PBPARM(3)/(10**3) + X*X*PBPARM(4)/(10**7) +
X*X*X*PBPARM(5)/(10**10) + X*X*X*X*PBPARM(6)/(10**13)
X*X*X*X*X*PBPARM(7)/(10**16)
where X is (distance from the pointing position in arc
minutes times the frequency in GHz)**2. For the VLA,
these parms are, by default, given by Perley's fits:
0.0738 GHz -0.897 2.71 -0.242
0.3275 -0.935 3.23 -0.378
1.465 -1.343 6.579 -1.186
4.885 -1.372 6.940 -1.309
8.435 -1.306 6.253 -1.100
14.965 -1.305 6.155 -1.030
22.485 -1.417 7.332 -1.352
43.315 -1.321 6.185 -0.983
For the ATCA, these are by default:
1.5 GHz -1.049 4.238 -0.8473 0.09073 -5.004E-3
2.35 -0.9942 3.932 -0.7772 0.08239 -4.429E-3
5.5 -1.075 4.651 -1.035 0.12274 -6.125E-3
8.6 -0.9778 3.875 -0.8068 0.09414 -5.841E-3
20.5 -0.9579 3.228 -0.3807 0.0 0.0
See explain for details
EXPLAIN SECTION
SAD: Task to find potential sources and fit Gaussian models to an image
by least squares.
DOCUMENTORS: W. Jaffe, Leiden
RELATED PROGRAMS : SLFIT, IMEAN, MAXFIT, IMFIT, JMFIT
PURPOSE
SAD attempts to find all sources in a subimage whose peak is brighter
than a given level. It searches the subimage specified by BLC and TRC
for all points above this level and merges such points in contiguous
"islands". For each island, initial estimates of the strength and size
are generated. Then the gaussian fitting algorithm used in JMFIT is
called to determine the least square guassian fit. If DOALL is < 0.,
only one gaussian is fit per island, with initial estimates generated
from 2nd moments. If DOALL is TRUE and multiple peaks (above a cutoff
= CPARM(i)) are found within the island, then up to 4 gaussians are
fit. If the peak residual of a single Gaussian fit exceeds ICUT then
a more complex 2-Gaussian fit ai attempted. The results of the fit
are printed on the screen or the line printer and written into a
ModelFit (MF) extension table. If DORESID is TRUE, the gaussians are
subtracted from the input subimage, and a residual imageis cataloged,
using OUTDISK, OUTNAME etc. If DOWIDTH(1) is negative, a gaussian with
fixed size = CLEAN beam will be fit.
WARNING!!!
The list generated by SAD is not COMPLETE in any statistical sense. You
will be warned about islands where the fitting algorithm failed. Note
that sizes etc. of sources near the noise are very unreliable. For
something approaching completeness Walter Jaffe recommends:
Start with a well Cleaned map
Set DOWIDTH = FALSE and only fit point sources.
Do a series of search on progressively more tapered maps and
compare results to determine flxes of extended sources.
Make residual maps and inspect them to see where SAD blew it.
Suggestion from Eric Greisen:
DOWIDTH false always produces poor results even for only slightly
extended objects. For decent S/N, the width fitting is quite reliable.
COMMENTS ABOUT SOME PARAMETERS
BLC, TRC:
These specify the area to be searched for sources. If an
residual file is requested, it covers the entire input image.
NGAUSS:
The maximum number of sources to find. Currently limited
to 40000.
CPARM:
Search for source peaks down to this level. If specified
as 0., SAD will determine the real (i.e. source-free) R.M.S. in
the full image area and take CPARM(1) to be 3.0 times this value.
If you make some CPARM(i) too small, you'll get lots of nonsense.
The task uses CPARM(1), then restarts using CPARM(2), etc so long as
the CPARM(i) > 0. This lets you extract the strongest points first
without low-level interference and then go after weaker things.
DOCRT:
DOCRT = 1. will cause a list of fluxes, positions, and
extension parameters, with their estimated errors, to be
displayed on your terminal.
DOCRT = -1. will cause the list to be printed on the line
printer.
DOCRT = 0. causes listing to be supressed.
DORESID
Catalog the post-fit residual map.
DOWIDTH
If you think most sources are unresolved, use DOWIDTH(1)
= -1.0. This will fit simple point sources to each peak (and do
very badly on those few that are slightly resolved).
COMMENTS ON THE USE OF SAD
FLUX DENSITY DETERMINATION:
When attempting to obtain the flux density of a well-resolved
source, the task IMEAN, which integrates the map values in a specified
rectangle, is often more accurate than fitting the source with several
Gaussian components and summing the integrated flux densities. Of
course for a well resolved source a gaussian fit is only a crude
approximation in any case. TVSTAT allows you to use an irregular area
for this determination.
PEAK FLUX DENSITY DETERMINATION:
The verb MAXFIT, a simple fitting of the peak of a component with a
second degree interpolation, is much faster than JMFIT and useful to
obtain the approximate peak and position of a component.
ERRORS OF PARAMETERS:
The error estimates should be regarded as tentative, but more
relaible than those from formal standard error estimates. The latter
are meaningless since the objects are normally well fit using only 1-4
Clean beam areas. This means that the number of parameters fit exceeds
the number of truly independent samples (neglecting the samples of 0
which were part of the selection of the island but not part of the
formal fit of course). The errors are produced using an estimate of the
true (signal free) rms R. It is taken from the image header keyword
ACTNOISE, if present and positive, or by fitting the histogram of the
image with a Gaussian. Note that this fit excludes both blanked and
pure zero pixels, thereby allowing blanked pixels to be REMAGed to
zero when needed. Using, John Ball's article in Method in
Computational Physics, 14, 177, 1975, the noise in the peak of the fit
Gaussian is estimated as:
Delta(P) = R * (Clean-beam area) / (fit area)
the noise in the widths as:
Delta(W) = Delta(p) / P * W
the noise in the center positions as:
Delta(X) = Delta(W) / 2
The noise in the total flux is found by the usual combinatorial formulae
and the noise in the position angle is estimated from Condon's article
Delta(PA) = Delta(p) / P * sqrt (2 WxWn / (WxWx - WnWn))
where Wx is the maximum width and Wn the minimum width.
The estimates from Ball have been replace by more complicated
formulae from Condon. The following remarks now apply. An estimate of
the each error is determined from theory based on the actual rms of the
image (neglecting signal portions) or the rms given in header parameter
ACTNOISE (if present and positive). AIPS verb ACTNOISE may be used to
set this header parameter. See the explain file for SAD for details.
Theory gives expressions for the errors in two limiting cases: point
source (the beam area > 0.9*the fitted gaussian area) and expanded
source (the beam area < 0.1*the fitted gaussian area). The formulae are
taken from J. Condon paper 'Errors in Elliptical Gausian Fit', AA, 1996.
The intermediate case is handled by interpolation between the two limit
cases. It is not clear if this interpolation provides a good estimation
of the error in general relation of beam and fitted gausian sizes and
position angle.
DECONVOLUTION:
When fitting to a Clean map, SAD deconvolves the Clean beam from
the fitted component size. The nominal deconvolution is obtained by
deconvolving the fit from the Clean beam. A value of 0.0 means that the
source is smaller than the Clean beam in some dimension. The minimum
and maximum values are obtained by deconvolving the source beam
parameters with all combinations of 0.7 * error and listing the extreme
values.
ON THE MESSAGE OUTPUT FORMAT
For readability and compactness the message formats are chosen to
fit in the width of the "printer" (real printer or CRT screen width
given by DOCRT). The UNITS of the printed flux are scaled so that the
weakest found flux lies in the range 1.0 to 999.0. The UNIT will be Jy
or milliJY or whatever (indicated in the header to the printout). If
another source has a flux greater than 999.999 in this UNIT, the flux
will be divided by 1000., before printing, and an asterisk "*" will be
printed before the flux to indicate this. A flux greater than
999,999.999 will be divided by 1.e6 and be prefixed by a "#".
If non-point sources are fit (DOWIDTH > 0), and a CLEAN beam size
can be found in the header, the derived source sizes and are
orientations are deconvolved, with the results printed in a second
listing. The numbers printed are the deconvolution of the best fit
major axis, minor axis and position angle. (arcsec, arcsec, degrees),
followed by the one-sigma lower limits on these quantities and the one
sigma upper limits.
ON THE MEANING OF BMAX, BMIN, AND BPA
Formally the position angle of anything, such as the major axis of
a gaussian source, is defined as the angle of the thing East from North.
But the direction North only coincides with a column of the pixel array
when you are on a meridian through the phase center. Thus in general
two different positional angles can be defined, either relative to local
North, or relative to the local Y-direction of the pixel matrix.
Similar small differences can occur in the definition of an astronomical
major axis versus a pixel-based major axis. Note for example that the
formal astronomical position angle of the point spread function of an
unresolved source changes with position in the field, due to the
curvature of the North/South coordinate system.
Sooooooooo, following the convention in JMFIT, the extention parameters
bpa, bmin, and bmax printed in the message file are the astronomical
values, relative to local North etc., while the parameters put into the
CC file are the pixel based values. This latter so that operations like
UVSUB will work correctly. This has the consequence, for example, that
for point sources the message file values differ from the values for the
BEAM, whereas the CC file values should agree with the BEAM. The MF
file contains everything and the kitchen sink, so all parameters are
present there.