NAME

specinholucy -- Extract point source spectra from longslit spectrum image with inhomogeneous background spectrum using an iterative technique based on Richardson-Lucy restoration

USAGE

specinholucy inpima errima dqima poitab Rootab poi_back bakima bakerr bakdq psfima psfmeth psfb interpol posmeth icsum skernel niter epspoi epsbac ntrial seed dqlim verbose

PARAMETERS

inpima = "" [file name]
The name of the input longslit spectrum image.
errima = "" [file name]
The name of the input statistical error image corresponding to the longslit spectrum image. If this parameter is set to null, multiple restoration trials cannot performed and no output error image will be produced.
dqima = "" [file name]
The name of the input data quality image corresponding to the longslit spectrum image. If set to null, no output data quality image will be produced.
poitab = "" [table file name]
Name of the STSDAS table file giving the X and Y value of the position of the point source and the slope (pixel^-1) with respect to the X (dispersion) axis.
Rootab = "" [table file name]
The root name of the output tables for the extracted point source spectra. The tables are post-fixed by an underscore and the item number of the profile as it appears in the input point source position table. This is a four column table with the X channel number (or X value if a World Coordinate System is found in the header), integrated point source flux, statistical error on the flux and the data quality listed.
poi_back = "" [boolean]
Switch to control the form of the output image(s). If switched on the output image (and error image if specified) consists of the restored background with the restored point source spectra (FWHM corresponding to the PSF); if switched off, the background image(s) only contains the restored background.
bakima = "' [file name]
Name of the output image for the restored background or background + point sources.
bakerr = "" [file name]
Name of the output statistical error image for the restored background or background + point sources. If no statistical error image was input, then no output error file is produced.
bakdq = "" [file name]
Name of the ouput data quality image for the restored background or background + point sources. If no data quality image was input, then no output data quality file is produced.
psfima = "" [file name]
The name of the input Point Spread Function (PSF) longslit spectrum image corresponding to the input spectrum image. The X dimension of this image must be the same as the input image.
psfmeth = "" [string]
Method used to fit the position of the centre of the Point Spread Function in the spatial (Y) direction. The options are C (Centroid of the profile) or G (Gaussian fit).
psfb = "" [integer]
The subsampling of the supplied Point Spread Function image in the Y (cross-dispersion) direction relative to the data image.
interpol = "" [string]
Method used to interpolate the Point Spread Function when shifting it to the position of the point source prior to restoration. The options are:
nearest - nearest neighbor interpolation
linear - linear interpolation
poly3 - third order interior polynomial
poly5 - fifth order interior polynomial
spline3 - cubic spline
posmeth = "" [string]
Method used to determine the position of the centre of the point source in the Y (spatial) direction. The options are E (Exact position from the input position table) or C (by Cross-correlation with the Point Spread Function).
icsum = "" [integer]
Number of X (dispersion) channels to sum in order to form a profile for determining the position in the Y (cross-dispersion) direction by cross-correlation with the Point Spread Function.
skernel = "" [float]
The standard deviation of the one-dimensional Gaussian smoothing kernel in the Y (cross-dispersion) direction used to smooth the iterated estimate of the background.
niter = "" [integer]
The number of iterations of the Lucy restoration for the point sources and background at each wavelength increment
epspoi = "" [float]
If in successive calls of the iterative restoration of the point sources and background for each wavelength increment, the estimated fluxes of all the point sources, do not differ by more this fractional value, then the restoration is converged.
epsbac = "" [float]
If in successive calls of the iterative restoration for the point sources and background for each wavelength increment, the estimate of the background does not change by more than this fractional value, then the restoration of the background is converged.
ntrial = "" [integer]
Number of Monte-Carlo trials for estimating the errors on the extracted point source spectrum and fitted background. This parameter can only be applied if the statistical error frame is input.
seed = "" [integer]
Value of the seed for the random number generator for producing Gaussian distributed errors for Monte Carlo simulation of statistical errors in the extracted spectra and restored background.
dqlim = "" [integer]
Limiting value of data quality above which to ignore data in input image
verbose = "" [boolean]
Switch to control output of information to monitor the progress of the iterations. If switched on, information on the extracted point source and background and their fractional changes per iteration are given as the iterations progress. See Example 2 below for the form of the output. If switched off, no progress report is output.

DESCRIPTION

This task applies a point-source and background spectral decomposition, using a restoration method developed by Lucy and Walsh, to each column of a 2D spectrum. The restoration at each spectral element is independent and the technique can be applied to a 2D spectrum where the background spectrum is spatially varying along the slit; this is referred to as an inhomogeneous 2D spectrum. Each spectral element is individually restored as a sum of point sources, at the specified positions, and a background, which is distinguished from the point sources by a Gaussian smoothing kernel. The iteration procedure restores the point sources and background of all the spectral increments each iteration. The restoration technique is Richardson-Lucy (viz. convergence to Maximum Likelihood).
The X direction of the input image is the dispersion direction and Y the spatial (i.e. cross-dispersion) direction. The point sources are specified only by their Y position. The Point Spread Function (PSF) for each column is taken from a PSF spectrum image (i.e. dispersion along the X axis and spatial direction along the Y axis) which must have the same number of X-channels as the input image but the 2nd (spatial, Y) dimension can be smaller than the input data image.
The point sources are indicated by the three columns of an STSDAS table. Only a single value is required to specify the Y position of the spectrum of a point source in the 2-D spectrum image, but since spectra with a tilt can be accommodated then both X and Y values are required. The three columns are:
Column 1    X coordinate of the point source spectrum (pixels)
Column 2    Y coordinate of the point source spectrum (pixels)
Column 3    Tilt of the spectrum with respect to the X axis (pixel^-1)
The X and Y positions use the IRAF convention that the bottom left corner pixel is (1,1) and its centre is at (1.0,1.0). The positions of the point sources are not restricted to pixel centres.
If the positions of the point sources are exactly specified, the Y pixel offsets with respect to the PSF are used in shifting the PSF to the point source position. If the positions of the point sources are not exactly known then they can be computed by the task by cross-correlating the point source spectrum and the object spectrum and fitting a Gaussian to the cross-correlation peak. The centre of the PSF is computed either by simple centroid or by a Gaussian fit. For symmetrical profiles these two methods give identical results.
Before the iterations commence, the columns of the PSF spectrum image are extracted, shifted according to the offsets of all the point sources, and an array of shifted PSF spectra is produced. The PSF can be oversampled (in the spatial dimension) with respect to the input data image, which is useful when the PSF of the data is undersampled.
The initial values for the iterations use a constant value for the flux in all the point sources as a function of wavelength and a flat image for the background.
At each iteration, the shifted versions of the PSF are scaled by the current estimates for the fluxes of the point sources and added together with the estimate of the background to create the current estimated image. For each X-section, this estimate is iterated until the maximum number of iterations or until the fractional change in the flux of all restored point sources is less than ESPOI and the fractional change in the 2D background is less than EPSBAC. Thus full convergence is deemed to occur when the two convergence criteria are met: i.e.
a) the fractional change in the flux of each point source per iteration is < EPSPOI
b) the fractional change in the flux of the continuum channels per iteration is < EPSBCK
Both these convergence criteria apply to all spectral channels combined. An example of the output showing the progress of the iterations is given below.
For the point sources, the output tables list, for each column of the input image, the X value of the column (wavelength if CRVAL1, CD1_1 were present in the header of the input image), the flux in the point source and the statistical error (if an input statistical error image was present and the number of Monte-Carlo trials was >1).
If the number of trials is set greater than 1 and a statistical error spectrum image is input, then an estimate of the corresponding statistical error on the restored longslit spectrum image can be obtained. Multiple trials are performed on the spectrum formed from the input image and the error value randomly chosen from a Gaussian distribution using a random number seed. The output error is the rms on the mean of the multiple trials.
Points in the input image which have flagged data quality (i.e. a data quality value >= dqlim) are not used in the restoration for iteration to the solution. The output restored data quality image is a copy of the input data quality image.

EXAMPLES

1. To extract a single point source from an image with no associated statistical error or data quality files:
       inpima = "testsp1.fits"  Name of input flux file
       errima = ""              Name of input statistical error file
        dqima = ""              Name of input data quality file
       poitab = "testpo1.tab"   Name of table of X,Y pixel centres and slope of point source spectra

                                OUTPUT tables and images
       Rootab = "test1t"        Root name for output point source extracted spectrum tables
     poi_back = yes             Output fitted point sources+background (Y) or background only (N)
       bakima = "testout1.fits" Name of output fitted point source+background or background file
       bakerr = ""              Name of output fitted point source+background or background statistical error file
        bakdq = ""              Name of output fitted point source+background or background data quality file

                                POINT SPREAD FUNCTION Spectrum and Position
       psfima = "psf1.fits"     Name of input Point Spread Function spectrum file
      psfmeth = "c"             Method for fitting PSF peak (C=Centroid; G=Gaussian)
         psfb = 1               Subsampling factor (Y-direction) for PSF image
     interpol = "spline3"       Interpolation method used for shifting PSFs

                                POSITIONS of Point Source Spectra
      posmeth = "c"             Method for fitting input point source peak (E=exact; c=Cross-correlation)
        icsum = 1               Number of columns of input data to sum for cross-correlation analysis

                                BACKGROUND CHANNEL parameter
      skernel = 4.              Sigma for smoothing kernel for floating prior

                                ITERATIONS and Stopping Criteria
        niter = 200             Number of iterations for point source + background restoration
       epspoi = 0.001           Stopping criterion for fractional change in point sources per set of iterations
       epsbac = 0.001           Stopping criterion for fractional change in background per set of iterations

       ntrial = 1               No. of Monte Carlo trials for error estimation 
        (seed = 7)              Seed for random number generator
       (dqlim = 100)            Limiting data quality value above which to ignore data
     (verbose = no)             Report parameters of restoration during execution
The input PSF spectrum has the same sampling as the input data. The input table lists the X,Y position and slope of the spectrum of the point source:
128.000    138.000      0.000
The PSF centre is determined by centroid fitting and the position of the point source spectrum in each column will be determined by cross-correlation with the PSF. Thus the value of Y in the input position table does not need to be exact (a value within a few pixels is adequate). The input has high signal-to-noise and only 1 column is to be summed before determining its position by cross-correlation. Two hundred iterations for the point-source and background will be made unless convergence is achieved. Convergence is achieved if the fractional change in the flux of the point source is < 0.001 and the fractional change in the background is less than 0.001 for each iteration of each X-section. The PSF will be shifted to the point source position using cubic spline interpolation. The smoothing kernel has a sigma of 4.0 and since no statistical errors are specified multiple trials are not performed. The output table is test1t_1.tab. No input data quality was specified so no data quality is output.
2. To extract three point sources, which are tilted, from an image and to apply multiple trials to estimate the errors:
cl>  specpoint inpima="testspec1.fits" errima="testspec1er.fits"
     dqima="testspec1dq.fits" poitab="demopos.tab" Rootab="test1t" 
     poi_back=yes bakima="testout1.fits" bakerr="testout1er.fits" 
     bakdq="testout1dq.fits" psfima="testpsf1.fits" psfmeth="g" 
     psfb=2 interpol="spline3" posmeth="e" icsum=10 skernel=8.0 
     niter=50 epspoi=5.0E-4 epsbac=5.0E-4 ntrial=100 seed=13 
     dqlim=10 verbose=yes 
The input PSF spectrum is subsampled by 2 with respect to the input data. The input table lists the X,Y position and slope of the spectrum of the point sources:
    128.000    120.500      0.000
    128.000    199.330      0.044
    128.000     88.890     -0.044
The PSF centre is to be determined by fitting a Gaussian and the position of the point source spectrum in each column is taken from the input position table. The number of channels to sum for cross-correlation is ignored. For each spectral element 50 iterations will be performed or until the fractional change in all the point sources is less than 0.0005 per iteration and the fractional change per iteration in the background is less than 0.0005. The PSF will be shifted to the point source position using a cubic spline. The smoothing kernel has a sigma of 8.0. 100 multiple trials are performed and output statistical error and data quality frames are written. Report on the progress of the iterations (for each column) will be written to the terminal.
!--------------------------------------------------------------------------------
!-- Starting iteration no.   39 
!-- Point sources and background spectra NOT converged
!-- Point source no.    1; % flux change    0.03808
!-- Point source no.    2; % flux change    0.02767
!-- Point source no.    3; % flux change    0.04788
!-- All point sources: Maximum % flux change    0.04788
!-- Background spectrum: % total flux change    0.26408
 Star #   Flux (#  64) / Back       Flux (# 128) / Back       Flux (# 192) / Back
    1     4877.8       105.59       5009.6       69.850       4822.2       101.27
    2     534.68       3.5033       500.60       10.941       534.53       2.5654
    3     945.39       5.7587       964.88       17.673       1032.5       3.8568
!-- Completed combined iteration no.   39
!--------------------------------------------------------------------------------
The lines show the results from the 39th iteration on the restoration of the three point sources and background at all X-channels to reach the fractional change per iteration in the restored point sources of 0.0005 and 0.0005 in the background for all wavelength increments. The convergence of all the point sources and the background for all spectral elements is reported. The fractional change in flux of each point source through this iteration is then listed in % and the maximum change. Then the fractional change in the background for all the wavelength increments is listed (in %). The following lines then list, for each point source spectrum, the flux and the flux in the mean background under the point source spectrum at three X-channels, selected at 0.25, 0.5 and 0.75 of the X range.

TIMING

The processing is dominated by convolutions and scales roughly as the number of pixels. The processing time is also proportional to the number of point source spectra to be extracted.

BUGS/FEATURES

1. The input longslit spectrum and PSF images must have a Y dimension which is a multiple of two (for the FFT's to succeed).
2. If there is more than one point source which has a small separation (less than order three times the FWHM of each point source) and their positions are determined through fitting of the cross-correlation peak ("posmeth=c"), then the determination of the positions may be in error. Either only one peak will be found or two peaks with incorrect positions depending on the signal-to-noise and the separation of the profiles in pixels.
3. If any X channel of a point source spectrum falls within 2 pixels of the (Y) top or bottom of the image, then that source is not considered for extraction. No output spectrum table will be produced for such spectra, but the table numbering will correspond to the original input table (i.e. if first spectrum falls at Y<2, then first output spectrum is named Rootab_2.tab).

REFERENCES

An introduction is provided in:
Walsh, J. R. "Long slit spectral extraction using restoration". ST-ECF Newsletter, 28, pp. 5-6, 2001.
and the technique is fully described in
Lucy, L. B., Walsh, J. R., 2002. "Iterative techniques for the decomposition of long-slit spectra". Submitted to AJ.

NAME

USAGE

PARAMETERS

DESCRIPTION

EXAMPLES

TIMING

BUGS/FEATURES

REFERENCES

SEE ALSO