AIPS HELP file for JMFIT in 31DEC24
As of Thu Oct 10 7:45:41 2024
JMFIT: 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 #
BLC 0.0 4096.0 Bottom Left corner of fit
TRC 0.0 4096.0 Top Right corner of fit
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 4.0 Number of components
CTYPE 0.0 5.0 Model types; one for each
component
0->1=Gaussian
2=Zero level; 3=Zero+slope
4=Zero+slope+curve;
5=see HELP JMFIT
Guess of model parameters
GMAX Peak of component (JY)
0-> Use maximum value
GPOS -50.0 16434.0 (X,Y) position (pixels)
0-> Use position of max
GWIDTH -180.0 180.0 (BMAJ, BMIN, PA) of comp.
(pixels,pixels,deg)
0->Use clean beam
@ Fit model parameters
FMAX @ Peak of component (JY)
@ 0-> Use maximum value
FPOS @ -50.0 16434.0 (X,Y) position (pixels)
@ 0-> Use position of max
FWIDTH @ -180.0 180.0 (BMAJ, BMIN, PA) of comp.
@ (pixels,pixels,deg)
@ 0->Use clean beam
FSHIFT @ RA/Dec shift to put comp 1
@ on nearest pixel
DOMAX $ Solve for GMAX? >0 -> yes
$ returns the error
DOPOS $ Solve for GPOS? >0 -> yes
$ returns the error
DOWIDTH $ Solve for GWIDTH? >0-> yes
$ returns the error
BWSMEAR 0.0 0.1 Bandwidth smearing corr.
RADIUS Radius for finding RMS
NITER 0.0 4000.0 Maximum # of iterations
0->20*NGAUSS
Solve for model parameters?
DOPRINT -4.0 132.0 <=0 -> Print maps and
solutions on line printer
>0 -> print on terminal
=0 -> maps not printed
NDIG Number digits in printed
maps (3 - 6, else fit)
FITOUT
Disk file to save fit info
DOOUTPUT -1.0 2.0 >0 -> 1 Catalog residual map
OFFSET -1.0 1.0 Cutoff level. 0-> None
DOMODEL -1.0 1.0 > 0 => put solutions in a CC
file with input image
OUTVERS -1.0 MF table version number
-1 => none, 0 => new
STVERS -1.0 STars output file version
-1 => none, 0 => new
PRTLEV Level of messages desired:
0 normal, 1 some extra,
2 iteration results
PBPARM Primary beam parameters:
(1) level to believe - <= 0
means do not apply a primary
beam (2) > 0 use (3)-(7)
EFACTOR 0.0 20.0 Scale width sigmas in the
deconvolution tests
HELP SECTION
JMFIT
Type: Task
Use: JMFIT is a task to fit up to four Gaussian components to a portion
of an image. It can also solve for a constant, linear, or
quadratic two-dimensional ``baseline'' surface. The program
estimates the error in the fits using the image rms and theory.
The rms is determined from the image header keyword ACTNOISE (if
present and positive) or from a robust determination of the rms
using the full (or a circular portion) image plane. Note that
pixels which are exactly zero are not used in this fit, allowing
blanked pixels to be REMAGed to zero if you wish.
The answers are returned in the parameters FMAX, FPOS, FWIDTH
and the uncertainties are returned in DOMAX, DOPOS, and DOWIDTH.
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.
BLC.........Bottom left corner of area of image to fit.
TRC.........Top right corner of area of image to fit.
Maximum area is 40000 pixels (200x200)
If the data are a cube (BLC(3) > 1), then the plane number
if displayed. For "standard" spectral-line cubes (XYV)
this will be the spectral channel. NOTE: there is no
requirement that the data be in standard coordinates or
transposition - although am elliptical Gaussian in
RA-velocity is a curious concept.
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 number of components to use in the fitting.
0->1. Maximum number is four.
CTYPE.......Each component type.
0=>1. Two-dimensional elliptical Gaussian
2=Solve for zero level
3=Solve for zero level and slope
4=Solve for zero, slope and curvature
5=Set the six baseline parameters as desired.
See EXPLAIN for use of GMAX,GPOS and GWIDTH
GMAX........The peak value guess for each component.
0=>Use value with largest absolute value in the
BLC,TRC window
GPOS........The position (X,Y) guess for components. The values are
in pixels in the order (X1,Y1,X2,Y2,X3,Y3,X4,Y4).
0=>Use pixel location with largest absolute value.
GWIDTH......The major axis, minor axis and position of major axis
guess for components. The values are pixels with
degrees for position angle and the order is
(MJ1,MN1,PA1,MJ2,MN2,PA2,...etc) 0->Use clean beam, if
available; otherwise it will be taken as a circular
gaussian of diameter 2.
********* The following are output adverbs:
FMAX........The peak value fit for each component is returned.
FPOS........The position (X,Y) fit for components is returned. The
values are in pixels in the order
(X1,Y1,X2,Y2,X3,Y3,X4,Y4).
FWIDTH......The major axis, minor axis and position of major axis
fit for components are returned. The values are pixels
with degrees for position angle and the order is
(MJ1,MN1,PA1,MJ2,MN2,PA2,...etc)
FSHIFT......The values that one could put in RASHIFT and DECSHIFT
to move the first component to the nearest integer
pixel. These values include any previous shifts and are
corrected for rotation.
********* The following are input/output adverbs:
DOMAX.......Flags for GMAX: if > 0 fit this parameter, else hold
fixed. Returned value is uncertainty in the fit
parameter.
DOPOS.......Flags for GPOS: if > 0 fit this parameter, else hold
fixed. Returned value is uncertainty in the fit
parameter.
DOWIDTH.....Flags for GWIDTH: if > 0 fit this parameter, else hold
fixed. Returned value is uncertainty in the fit
parameter.
********* The following are input adverbs:
BWSMEAR.....If > 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 BWSMEAR. Set it roughly
to the channel bandwidth divided by the center frequency.
RADIUS......If = 0, the rms used for error estimates is taken from
the image header (keyword ACTNOISE) or found by fitting
the full image plane. If RADIUS > 0, then the rms is
found by fitting only those pixels within RADIUS of the
center of the BLC-TRC box. If RADIUS < 0, then
abs(RADIUS) is used as the rms. Pixels which are
exactly zero are not used in the fitting for rms. A
robust method is used and, if it fails within RADIUS,
the method will be applied to the full image plane.
NITER.......The maximum number of iterations to use in the
fitting. 0-> NGAUSS * 20.
DOPRINT.....<0 -> Plot map, model and residual map and list fit
information on the line printer
When FITOUT is not blank, DOPRINT=-2 suppresses
the page-feed character on page headers and
DOPRINT=-3 suppresses page headers and most other
header information.
DOPRINT=-4 produces output to FITOUT with one line
per source component in a special CfA format.
=0 -> List fit info in message file only
>0 -> Plot map, model, and residual map and list fit
information on the terminal
NDIG........Number of digits in printed maps. If NDIG = 3-6 and is
less than or equal the number that will fit, then NDIG
is used.
FITOUT......Disk file name in which to save the 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.
DOOUTPUT....>0 -> Catalog residual map with fitted components
and write them in a CC file attached to the output
image. If DOOUTPUT > 1.5, the components written
to the CC file are not deconvolved from the beam.
OFFSET......0-> Include all points in fitting area. Otherwise
disregard all points less than OFFSET*MAX, where
MAX is largest value in fitting window. If MAX
is less than zero, disregard all point greater
than OFFSET*MAX.
DOMODEL.....If true (> 0), put the deconvolved solutions in a
new CC file attached to the input image.
OUTVERS.....The results are written into an MF (Model Fit) table file
with version OUTVERS unless OUTVERS is set < 0. If
OUTVERS points at a pre-existing table, the results are
appended to the file. OUTVERS = 0 always means to make a
new MF table.
STVERS......The results may also be written to a STars table
attached to the image. -1 => do not do this, 0 => write
a new one, > 0 => add to the specified version.
PRTLEV......Print level: 0 => normal, 1 => some from fitting,
2 => all from fitting (to diagnose problems)
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.
See EXPLAIN for details and defaults
EFACTOR.....The range of possible deconvolved widths is found by
trying the deconvolution with Major, Minor, and PA each
-FACTOR*sigma, 0, +FACTOR*sigma from the fit value. The
highest value found in the 27 tests and the lowest value
found are reported as well as the values found at the
nominal (fit) values. 0 -> 1.3.
EXPLAIN SECTION
JMFIT: Task to fit Gaussian models to an image by least squares.
DOCUMENTORS: E.B.Fomalont NRAO/VLA, F. Schwab NRAO, CV
RELATED PROGRAMS : SLFIT,IMEAN,MAXFIT,IMFIT,SAD
PURPOSE
JMFIT fits up to four Gaussian-shaped components to a selected
part of an image. One of the components can be a baseline function
with a zero level, slope and curvature term. JMFIT is most commonly
used to derive the position, peak and integrated intensity and angular
size of a source which is not too extended. An initial guess for the
parameters, some of which are picked as defaults, must be supplied
before running the task. Solution and error estimates are generated
and the residual image after the fit can be printed on the
line-printer. An arbitrary selection of parameters may be held
constant in the solution. Default initial guesses chosen by the task
may suffice in the case of a single-component model; otherwise,
sufficiently good guesses for at least the positions must be provided
by the user.
The fitting algorithm is one which is based on a algorithm of
Davidon. Occasionally, the solution will converge on an obviously
unacceptable fit. If this occurs, then try a better initial guess.
When fitting several components to a blobby source, the fitted
parameters may be absurd. Careful selection of fixed parameters then
is necessary.
COMMENTS ABOUT SOME PARAMETERS
BLC, TRC:
The fitting area should be chosen as small as possible; and
several disconnected components should be fit separately. The
fitting area is limited to an area of 10000 pixels. If the data are a
cube (BLC(3) > 1), then the plane number if displayed. For "standard"
spectral-line cubes (XYV) this will be the spectral channel.
NGAUSS:
The number of components to fit. The maximum number is four, and
0->1.
CTYPE:
The component types, placed in a scalar array of length 4.
0->1 Elliptical Gaussian component.
2 Zero level.
3 Zero level and slope.
4 Zero level, slope and curvature.
5 Insert baseline parameters as follows:
GMAX = zero level
GPOS = slope (intensity per pixel),
Orientation (deg, measured from N through E)
GMAX = Curvature (intensity per square pixel),
Ellipticity of curvature (-1 to 1),
Orientation (deg, measured from N through E)
DOMAX, DOPOS and DOMAX are used to hold parameters fixed
for Types 1 and 5 only.
GMAX:
The initial guess of the model intensity may be supplied. The
units of GMAX are the same as those in the map. The default of 0 will
place the most extreme value in the fitting area (negative or
positive) in GMAX for the first component. Any subsequent components
with 0 default will be given the value of 0.1 times the extreme value.
GPOS:
The initial guess of the model position. The location must be
given by a pair of pixel coordinates. The default of 0 will insert
the location of the extreme value for all components. Note that GPOS
has the meaning of slope and orientation for a baseline component.
GWIDTH:
The defaults for the component widths are generally reasonable:
either the clean beam size or a circular beam of two pixels FWHP. To
try to get the minimization algorithm off to a faster start for
circular Gaussian models, the task will introduce a slight ellipticity
before beginning the fitting. This is not done if either axis is to
be held constant.
NITER:
The maximum number of iterations, NITER, is defaulted to 20 times
NGAUSS if the input value is 0. Whenever convergence is not achieved
in NITER iterations, JMFIT so informs the user and ceases execution.
This is to discourage the user from accepting a poor solution. When
convergence is not achieved, the user should increase NITER if it
appears that the program may have been converging on an acceptable
solution, but doing so slowly. Otherwise, one should alter the
initial guesses, try fitting to data within a somewhat smaller or
larger window, or alter the value of OFFSET.
DOWIDTH:
In fitting complicated sources, it is common to hold some of the
component diameters fixed.
DOPRINT:
Set DOPRINT = -1 in most cases. This produces an automatic hard
copy of the solutions and a digital map of the input image, the
initial guess and the post-fit residuals. These maps are useful for
assessing the quality of the fit. When DOPRINT = -1, the progress of
the (iterative) minimization algorithm is reported on the line-printer
output. Otherwise, it is displayed on the terminal.
DOOUTPUT:
Catalog the post-fit residual map. If all parameters are held
fixed, no fitting is done and a residual map is generated Only the
fitted area is cataloged. A CC file is also written with this output
image listing the components. If DOUTPUT > 1.5, the components
written are not deconvolved from the beam. If 0 < DOOUTPUT <= 1.5,
the components are deconvolved if possible.
OFFSET:
This adverb permits the exclusion of low valued points when doing
the fit. If the extremum value (MAX) in the fitted area is positive,
then all points less than OFFSET*MAX are ignored in the fit. If the
extremum value in the fitted area is negative, then all points greater
than OFFSET*MAX are ignored in the fit. If OFFSET = 0, then all
points are used.
PRIMARY BEAM CORRECTION
JMFIT corrects an image for the primary beam attenuation of
the antennas. The function used to model the primary beam for normal
VLA frequencies
F(x) = 1.0
+ parm(3) * 10E-3 * x
+ parm(4) * 10E-7 * x*x
+ parm(5) * 10E-10 * x*x*x
+ parm(6) * 10E-13 * x*x*x*x
+ parm(7) * 10E-16 * x*x*x*x*x
where x is proportional to the square of the distance from the
pointing position in units of [arcmin * freq (GHz)]**2, and F(x)
is the multiplicative factor to divide into the image intensity at the
distance parameter x. For other antennas, the user may read
in appropraite constants in PBPARM(3) through PBPARM(7). The
flag, PBPARM(2) must be set to a positive number to invoke this
option and PBPARM(3) must not be zero.
This correction scales with frequency and has a cutoff
beyond which the map values are set to an undefined pixel value GIVEN
in PBPARM(1). At the VLA frequencies the default cutoff is
1.485 GHz 29.8 arcmin
4.885 GHz 9.13 arcmin
15 GHz 2.95 arcmin
22.5 GHz 1.97 arcmin
and occurs at a primary beam sensitivity of 2.3 percent of the value at
the beam center. Corrections factors < 1 are forced to be 1.
The estimated error of the algorithm is about 0.02 in (1/F(x))
and thus leads to very large errors for x>1500, or at areas
outside of the primary response of 20 percent. The cutoff level
may be specified with DPARM(1).
Default values of PBPARM for the VLA are 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
For the Karl G Jansky VLA ("EVLA"), the defaults are frequency
dependent. If the observing frequency is between two tabulated
frequencies, then the beam is computed for each of the tabulated
frequencies and then interpolated to the observing frequency. The
values used are far too numerous to give here, see EVLA Memo 195,
"Jansky Very Large Array Primary Beam Characteristics" by Rick Perley,
revision dated June 2016. Obtain it from
http://library.nrao.edu/evla.shtml
COMMENTS ON THE USE OF JMFIT
SETTING UP THE PARAMETERS
For some simple cases the defaults set by JMFIT are good enough
starting values. Some examples are as follows: (Always insert the
appropriate input map and always set BLC and TRC to the smallest area
needed for the fit. The verb TVWINDOW can be used to set the window.)
1. Fit to one Gaussian component.
Nothing to specify except flags.
2. Fit to one Gaussian and zero level
NGAUSS=2;CTYPE=1,2;GO
Additional flags can be specified
3. Fit to one Gaussian, zero-level and slope
NGAUSS=2;CTYPE=1,3;GO
Additional flags can be specified
4. Fit to two Gaussians and a zero level
NGAUSS=3;CTYPE=1,1,2;GO
Depending on the source complexity it may be
important to set some of the fitting flags
5. Fit one Gaussian with a zero level, slope and
curvature only in E/W direction
NGAUSS=2;CTYPE=1,5;DOMAX=1,1;DOPOS=1,1,1,-1
DOWIDTH=1,1,1,1,-1,-1;GPOS=0,0,0,90;
GWIDTH=0,0,0,0,0,90
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.
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:
An estimate of the each error is determined from theory based on
the actual rms (R) of the image (neglecting signal portions) or the
rms given in header parameter ACTNOISE (if present and positive).
AIPS task IMEAN or verb ACTNOISE may be used to set this header
parameter. 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. The formulae now used are:
M = 1 (Clean beam area > 0.9 fit area)
M = sqrt (8 * ARbeam/ARimag) (Clean beam area < 0.1 fit area)
M = sqrt (0.8 + 0.25*(ARbeam/ARimag-0.1)) (else)
Then
Delta(P) = M * R
Delta(W) = W * Delta(P) / P
Delta(PA) = sqrt(2) * (Smay+Smin)/(Smay^2+Smin^2) * Delta(P) / P
Delta(X) = sqrt[(Delta(Smay)*sin(PA))^2 + (Delta(Smin)*cos(PA))^2]
/ sqrt (8 * ln (2.0))
Delta(Y) = sqrt[(Delta(Smay)*cos(PA))^2 + (Delta(Smin)*sin(PA))^2]
/ sqrt (8 * ln (2.0))
Delta(F) = Delta(P) * ARimag/ARbeam * sqrt(1+2*ARbeam/ARimag)
where P is fit peak, W = Smaj or Smin are fit widths, X and Y are
positions, PA is position angle, ARimag = Smaj*Smin, ARbeam=Bmaj*Bmin.
DECONVOLUTION:
When fitting to a clean map, JMFIT deconvolves the Clean beam from
the fitted component size. The nominal deconvolution is obtained by
deconvolving the fit from the Clean beam (corrected for bandwidth
smearing). A value of 0.0 means that the source is smaller than the
corrected Clean beam in some dimension. The minimum and maximum
values are obtained by deconvolving the source beam parameters with
all 27 combinations of EFACTOR * (-1, 0, 1) * uncertainties in the
major axis, minor axis, and position angle. The extremas in these
parameters over all 27 tries are listed. The default EFACTOR (1.3)
appears to work well with the considerations below.
An estimate is given concerning whether the source should be
viewed as resolved or unresolved. The task assumes that the
component is probably unresolved if:
(a) the deconvolution of the fit answers has the major axis 0
OR
(b) the fit total flux minus the error in the total flux is
less than the peak AND the minimum deconvolved major axis
is 0.
The task assumes the conponent is probably resolved if
(c) the total fit flux minus the error in the total fit flux is
greater than peak flux AND the minimum deconvolved major axis
is greater than 0.
The task is undecided about resolution if
(d) the deconvolution of the fit answers has the major axis greater
than zero.
AND
(e) the fit total flux minus the error in the total flux is
less than the peak BUT the minimum deconvolved major axis
is greater than 0.
OR
(f) the total fit flux minus the error in the total fit flux is
greater than peak flux BUT the minimum deconvolved major axis
is 0.
Note: the total flux and its error are corrected for primary beam and
the peak is corrected for primary beam and bandwidth smearing (if such
corrections are requested) in the tests described above. If the
component is unresolved, the best estimate of its total flux is its
peak brightness. In cases where the task is uncertain, use caution in
deciding if the component is resolved. Noise seems preferentially to
make sources appear resolved when they are not. Note too that a
"resolved" source may be clearly unresolved along the minor axis.
OUTPUT FROM THE MINIMIZATION SUBROUTINE
With iterative algorithms, it is important to monitor the
progress of the iteration. For this reason JMFIT, at each iteration,
prints out the value of the chi-squared function (the sum of the
squares of the residuals) and the Euclidean norm of the gradient of
the chi-squared function (it is trying to zero the gradient).
Chi-squared is labeled F in the printout, and the gradient norm is
labeled simply GRADIENT. The residuals are measured in the same units
as the units of the input data, but the units are not printed out by
JMFIT.
Also, at each iteration, the current estimates of the solution
parameters are printed out in their natural order (i.e., amplitude,
x-position, y-position, major axis fwhm, minor axis fwhm, and position
angle for Gaussian components), component-by-component. (Fixed
parameters are not included.) The units are input data amplitude
units for amplitudes, pixel coordinates for positions, pixels for
widths, and radians for angles. For the usual orientation of the RA
and DEC axes, pi/4 is added to the position angles in these printouts.
PRTLEV is provided to suppress some or all of this output. When
DOPRINT = -1 this information is printed on the line-printer.
Otherwise, it is printed on the terminal. Thus, if the user is
annoyed by the quantity of output which is written to the terminal, he
can set DOPRINT = -1 to have it instead written on the line-printer
(and vice-versa, of course). If PRTLEV < 2, the output is greatly
reduced.