fitprofs input
input
List of input images to be fit. The images may be one dimensional
spectra (one or more spectra per image) or long slit spectra. Other
types of nonspectral images may also be used and for two dimensional
images the fitting direction will be determined from either the keyword
DISPAXIS in the image header or the dispaxis parameter.
lines =
List of lines, columns, or apertures to be selected from the input image
format. The default empty list, "", selects all vectors in the images.
The syntax is a list of comma separated numbers or ranges, where a range
is a pair of hyphen separated numbers.
bands =
List of bands for 3D images. The empty list, "", selects all bands.
dispaxis = )_.dispaxis , nsum = )_.nsum
Parameters for defining vectors in 2D and 3D images. The
dispersion axis is 1 for line vectors, 2 for column vectors, and 3 for band
vectors. A DISPAXIS parameter in the image header has precedence over the
dispaxis parameter. The default values defer to the package
parameters of the same name.
The following are the fitting parameters.
region =
Region of the input vectors to be fit specified as a pair of space
separated numbers. The coordinates are defined in terms of the linear
image header coordinate parameters. For dispersion corrected spectra this
is usually wavelength in Angstroms and for other data it is usually pixels.
A fitting region must be specified.
positions =
File of initial or fixed profile positions and (optional) widths. The
format consists of lines with one or two whitespace separated numbers.
Comments and any additional columns are ignored. The positions and
widths are specified in the coordinate units of the image, usually
wavelength for dispersion corrected spectra and pixels otherwise. The
width is interpreted as a sigma for Gaussian profile fitting.
background =
Initial or fixed linear background defined by the region
endpoints. If not specified the pixel values nearest the fitting
region endpoints are used. Other values are specified by a pair of
space separated numbers. When fitting the background this parameter
defines the initial guess otherwise the background is fixed as
specified.
sigma = 5.
Default profile width or sigma. This value is used for the initial/fixed
width when a value is not specified in the position file.
function = gaussian
Profile function to be fit. Currently the only choice is "Gaussian" for
Gaussian profiles though additional functions may be added in future.
fitbackground = yes
Fit the background? If "yes" a linear background across the fitting region
will be fit simultaneously with the profiles. If "no" the background will
be fixed.
fitpositions = all
Position fitting option. This may be "fixed" to fix all positions at their
initial values, "single" to fit a single shift to the positions while
keeping their separations fixed, or "all" to independently fit all the
positions.
fitsigmas = all
Profile width fitting option. This may be "fixed" to fix all widths at their
initial values, "single" to fit a single scale factor to the initial
widths of the profiles, or "all" to independently fit all the widths.
The following parameters are used for error estimates in the 'd',
'k', and 'e' key measurements. See the ERROR ESTIMATES section for a
discussion of the error estimates.
sigma0 = INDEF, invgain = INDEF
The pixel sigmas are modeled by the formula:
sigma**2 = sigma0**2 + invgain * I
where I is the pixel value and "**2" means the square of the quantity. If either parameter is specified as INDEF or with a value less than zero then no sigma estimates are made and so no error estimates for the measured parameters is made.
The following parameters determine the output of the task.
components =
All profiles defined by the position file are simultaneously fit but only
a subset of the fitted profiles may be selected for output. A profile
or component is identified by the order number in the position file;
i.e. the first entry in the position file is 1, the second is 2, etc.
The components to be output are specified by a range list. The empty
list, "", selects all profiles.
verbose = yes
Print fitting results and record of output images created on the
standard output (normally the terminal).
The fitting informations is printed to the logfile so there is normally
no need to redirect this output. The output may be turned off when
the task is run as a background task.
logfile = logfile
Logfile for fitting results. If not specified the results will not be
logged.
plotfile = plotfile
File to contain plot output. The plots show the image vector with
overplots of the total fit, the individual components, and the residuals.
The plotfile may be examined and manipulated later with tools such as
gkimosaic.
output =
List of output images. If not specified then no output images are created.
If images are specified the list is matched with the input list.
option = fit
Image output option. The choices are "fit" to output the fitted image
vector which is the sum of the fitted profiles (without a background),
or "difference" to output the data with the profiles subtracted.
clobber = no, merge = no
Clobber or modify any existing output images? If clobbering is not
enabled a warning is printed and any existing output images are not
modified. If clobbering is enabled then either new images are created
if merge is "no" or the new fits are merged with the existing images.
Merging is meaningful when only a subset of the input is fit such
as selected lines or apertures.
Fitprofs fits one dimensional profile functions to image vectors and outputs the fitting parameters, plots, and model or residual image vectors. This is done noninteractively using a file of initial profile positions and widths. Interactive profile fitting may be done with the deblending option of splot or stsdas.fitting.ngaussfit.
The input consists of images in a variety of formats. These include all the spectral formats as well as standard images. For two dimensional images (or the first 2D plane of higher dimensional images) either the lines or columns may be fit with possible summing of adjacent lines or columns to increase the signal-to-noise. A subset of the image apertures, lines, or columns may be specified or all image vectors may be fit.
The fitting parameters consist of a fitting region, a list of initial positions and widths, initial background endpoints, the fitting function, and the parameters to be fit or constrained. The coordinates and units used for the positions and widths are those defined by the standard linear coordinate header parameters. For dispersion corrected spectra these are generally wavelengths in Angstroms and otherwise they are generally pixels. A fitting region must be specified by a pair of numbers. The background parameter may be left empty to select the pixel values at the endpoints of the fitting region for defining the initial linear background. The position list consists of one or two columns. The positions are the first column and are required while the second column gives the profile width which defaults to the sigma parameter if it is missing.
The task currently fits only Gaussian profiles with simultaneous linear background. Other functions may be added in the future. The profile parameters fit are the central position, the peak amplitude, and sigma. The fitting may be constrained in number of ways. The linear background may be fixed or simultaneously fit with the profiles. The profile positions may be fixed, the relative separations fixed but a single zero point shift fit, or all positions may be fit simultaneously. The profile sigmas may also be fixed, the relative ratios of the widths fixed while fitting a single scale factor, or all widths fit simultaneously. The profile amplitudes are always fit.
The fitting technique uses a nonlinear iterative Levenberg-Marquardt algorithm to reduce the Chi-square of the fit. The execution time increases rapidly with the number of profiles fit so there is an effective limit to the number of profiles that can be fit at once.
The output includes a number of formats. The fitted parameters are recorded in a logfile (if specified) and printed on the standard output (if the verbose flag is set). This output includes the date, image vector, fitting parameters used, and a table of fitted or derived quantities. The parameters included some quantities relevant to spectral lines but others apply to any image data. The quantities are the profile center, the background or continuum at the center of the profile, the integral or flux of the profile (which is negative for profiles below the background), the equivalent width, the profile peak amplitude or core value, and the Gaussian sigma and full width at half maximum.
Summary plots are recored in a plotfile (if specified). The plots show the data with the total fit, individual profiles, and residuals overplotted. The plotfile may be examined and printed using the task gkimosaic as well as other tasks which interpret GKI metacode.
The final output consists of images in the same format as the input. The images may be of the total fit (sum of profiles without background) or of the difference (residuals) of the data minus the model.
Error estimates may be computed for the fitted parameters. This requires a model for the pixel sigmas. Currently this model is based on a Poisson statistics model of the data. The model parameters are a constant Gaussian sigma and an "inverse gain" as specified by the parameters sigma0 and invgain. These parameters are used to compute the pixel value sigma from the following formula:
sigma**2 = sigma0**2 + invgain * I
where I is the pixel value and "**2" means the square of the quantity.
If either the constant sigma or the inverse gain are specified as INDEF or with values less than zero then no noise model is applied and no error estimates are computed. Note that for processed spectra this noise model will not generally be the same as the detector readout noise and gain. These parameters would need to be estimated in some way using the statistics of the spectrum. The use of an inverse gain rather than a direct gain was choosed to allow a value of zero for this parameters. This provides a model with constant uncertainties.
The error estimates are computed by Monte-Carlo simulation. The model is fit to the data (using the sigmas) and this model is used to describe the noise-free spectrum. One hundred simulations are created in which random Gaussian noise is added to the noise-free spectrum based on the pixel sigmas from the noise model. The model fitting is done for each simulation and the absolute deviation of each fitted parameter to model parameter is recorded. The error estimate for the each parameter is then the absolute deviation containing 68.3% of the parameter estimates. This corresponds to one sigma if the distribution of parameter estimates is Gaussian though this method does not assume this.
The Monte-Carlo technique automatically includes all effects of parameter correlations and does not depend on any approximations. However the computation of the errors does take a significant amount of time.
1. The following example creates an artificial spectrum and fits it. It requires the artdata and proto packages be loaded.
cl> mk1dspec test slope=1 temp=0 lines=testlines nl=20
cl> mknoise test rdnoise=10 poisson=yes
cl> fields testlines fields=1,3 > fitlines
cl> fitprofs test reg="4000 8000" pos=fitlines
# Jan 3 17:42 test - Ap 1:
# Nfit=20, bkgd=YES, positions=all, sigmas=all, RMS=177.073
# center cont flux eqw core sigma fwhm
6829.874 1362.309 -7159.78 5.256 -157.466 18.14 42.72
7963.207 1507.092 -9451.88 6.272 -376.275 10.02 23.6
5688.533 1216.502 -8024.02 6.596 -370.239 8.646 20.36
6832.274 1362.615 -16332.1 11.99 -626.033 10.41 24.51
7216.591 1411.712 -11661.5 8.261 -417.378 11.15 26.25
6708.711 1346.83 -5725.68 4.251 -222.683 10.26 24.16
6433.66 1311.692 -8449.34 6.442 -331.589 10.17 23.94
6130.577 1272.974 -7032.61 5.525 -228.346 12.29 28.94
4567.955 1073.349 -6480.78 6.038 -249.199 10.38 24.43
5657.427 1212.529 -9856.24 8.129 -311.541 12.62 29.72
4219.741 1028.864 -6138.24 5.966 -238.262 10.28 24.2
4546.741 1070.638 -2892.72 2.702 -139.746 8.258 19.45
4605.13 1078.098 -6454.01 5.986 -259.21 9.933 23.39
6966.801 1379.801 -13722.6 9.945 -572.377 9.565 22.52
4258.94 1033.872 -4701.89 4.548 -214.153 8.759 20.63
5952.413 1250.213 -9507.61 7.605 -318.475 11.91 28.05
4530.854 1068.609 -706.732 0.6614 -62.6292 4.502 10.6
7815.642 1488.241 -2176.27 1.462 -271.426 3.199 7.533
5311.103 1168.286 -11653.5 9.975 -481.713 9.651 22.73
5022.378 1131.401 -8689.48 7.68 -329.906 10.51 24.75
REVISIONS
FITPROFS V2.11
Error estimates based on a simple noise model are now computed.
FITPROFS V2.10
This task is new.
The following CPU times were obtained with a Sun Sparcstation I. The number of pixels in the fitting region and the number of lines fit were varied. The worst case of fitting all parameters and a background was considered as well as the constrained case of fitting line positions and a single width with fixed background.
Npixels Nprofs Fitbkg Fitpos Fitsig CPU(sec) 100 5 yes all all 1.9 100 10 yes all all 3.3 100 15 yes all all 5.6 100 20 yes all all 9.0 512 5 yes all all 4.7 512 10 yes all all 10.0 512 15 yes all all 17.6 512 20 yes all all 27.8 1000 5 yes all all 8.0 1000 10 yes all all 18.0 1000 15 yes all all 31.8 1000 20 yes all all 50.2 1000 25 yes all all 72.8 1000 30 yes all all 100.2 512 5 no all single 2.8 512 10 no all single 5.3 512 15 no all single 8.6 512 20 no all single 12.8
Crudely this implies CPU time goes as the 1.4 power of the number of profiles and the 0.75 power of the number of pixels.
splot, stsdas.fitting.ngaussfit,