diff options
Diffstat (limited to 'noao/twodspec/multispec/doc/MSalgo.ms')
| -rw-r--r-- | noao/twodspec/multispec/doc/MSalgo.ms | 1032 |
1 files changed, 1032 insertions, 0 deletions
diff --git a/noao/twodspec/multispec/doc/MSalgo.ms b/noao/twodspec/multispec/doc/MSalgo.ms new file mode 100644 index 00000000..0c64e2b3 --- /dev/null +++ b/noao/twodspec/multispec/doc/MSalgo.ms @@ -0,0 +1,1032 @@ +.de FX +.nf +.ps -2 +.ss 25 +.cs R 25 +.. +.de EX +.ps +2 +.ss +.cs R +.fi +.. +.EQ +delim $$ +.EN +.RP +.TL +Algorithms for the Multi-Spectra Extraction Package +.AU +Francisco Valdes +.K2 +.TU +.AB +The algorithms for the Multi-Spectra Extraction Package (\fBmultispec\fR) +in the Image Reduction and Analysis Facility (\fBIRAF\fR) is described. +The basic aspects of the general two dimensional aperture spectra extraction +problem are first discussed. +The specific algorithms for extraction of multi-aperture plate and +Echelle digital data are presented. Results of the authors experiments +with this type of data are included. +The detailed specification of the package is given in a second document, +\fIDetailed Specifications for the Multi-Spectra Extraction Package\fB. +.AE +.NH +Introduction +.PP +There are an increasing number of astronomical instruments which produce +multiple spectra on a two dimensional detector. +The basic concept is to use one dimension for wavelength, +the dispersion dimension, and the other, the cross dimension, for +packing additional information during a single exposure. +For example, the cross dimension can be different objects or +different spectral orders. The classic multi-spectra instrument is +the Echelle spectrograph. New instruments are the aperture plate and +Medusa spectrographs. +.PP +There is an additional aspect of the multi-spectra format; namely, +the individual spectra can contain spatial data. An example of +this would be multiple slit spectra in which each slit spectra contains +sky signal and object signal. In the following +discussion we limit the spectra to be simple aperture spectra in +which we only desire to sum the intensities at each wavelength to form +a one dimensional spectrum. +.PP +The analysis of multi-spectra aperture data consists of two steps; the +separation and extraction into individual aperture spectra +and the calibration and measurement of the spectra. These steps can +either be incorporated into one analysis package or two separate +packages. There are advantages to the first approach since some +aspects of the individual spectra are directly related by the physical +geometry of the multi-spectra format. However, because packages for +the analysis of individual spectra exist we begin by dividing the +reduction into separate extraction and analysis tasks. It is +important to realize, however, that the existing analysis tools are not well +suited to reducing the larger number of spectra and treating sets of +spectra together. +.PP +The latter part of this paper describes the algorithms for the +extraction of two types of data; the multi-aperture plate (MAP) +and Echelle used with digital detectors. However, +it is important to keep the more general problem in mind +and the remainder of this introduction considers the different +conceptual aspects of the multi-spectra extraction task. +Table 1 lists many of the general properties of multi-spectra aperture data. +The other two columns give possible alternatives that each property may take. + +.TS +center box; +c s s +c s s +c c s += = = +c || c | c. +Table 1: Aspects of Multi-Spectral Data + +Property Alternatives +detector digital photographic +alignment aligned skewed +blending blended unblended +aperture holes slits +spectral resolution low high +.TE + +.PP +The detector determines what kind of calibration procedures are +required to produce intensity from the measurements. +A digital detector requires sensitivity calibrations on all scales. +This is the "flat field" problem. There are also corrections for +bias and dark current. Photographic detectors require +intensity calibration. Data which are not aligned with the natural +dimensions of the digital image require extra procedures. Two types +of non-alignment are a rotation of the dispersion dimension relative +to the pixel dimension and a "wiggling" or "snaking" of the dispersion +dimension. Blending refers to the degree of contamination along the +cross dimension. Blended data requires extra effort to correct for +the overlap between different spectra and to determine the background. +The aperture defines the extent of the spectra in the cross dimension. +The two most relevant choices are holes and slits. In some +instruments, like the Echelle, the size of the aperture can be varied +at the time of the observations. Finally, the spectral resolution is +important in conjunction with digital detectors. If the resolution is +high then quartz flat field calibrations are relatively easy because +the spectral +signature need not be considered. Otherwise, the flat field problem +is more difficult because the gain variations of the detector +must be separated from the natural spectral intensity variation of the +quartz. +.PP +There is always some confusion of terms when talking about multi-spectra +data; in particular, the terms x, y, line, and band. +The image pixel dimensions are refered to as x and y. We assume +for the moment that the alignment of the multi-spectra format is such +that x corresponds to the cross dimension and y to the dispersion +dimension. If this is not the case a rotation or interpolation +program can be used to approximately orient the data in this way. +A line is the set of intensity values as a function of x at constant y. +In other words, a line is a cut across the dispersion dimension. +A band is the average of more than one line. +The image residing on disk will generally be organized +such that x varies more rapidly and a line of the image is easily +obtained. In this form a display of the image will have the spectra +running vertically. In the Cyber extraction package the data is +organized with x corresponding to the dispersion dimension. +.NH +Multi-Spectra Image Formats +.PP +The remainder of this paper will refer to two specfic and very +different multi-spectra formats; the Kitt Peak Multi-Aperture Plate +System and the Kitt Peak Echelle Spectrograph. +.NH 2 +Kitt Peak Multi-Aperture Plate System +.PP +The reduction of data from multi-aperture plate observations is the +driving force for the development of a multi-spectra extraction +package. This application turns out to have most of the worse aspects +of the properties listed in Table 1. The multi-aperture plate spectrograph uses +digital dectectors with low resolution, the spectra are blended and +change alignment along the pixel dimension. Furthermore, the camera +has a variable point-spread function and focus, +suffers from flexture problems, has a different illumination for +the quartz than object exposures, and unexplained background level +variations (in the CRYOCAM). There are two detectors which have been +used with the multi-aperture plate system, the Cryogenic Camera +(CRYOCAM) and the High Gain Video Spectrometer (HGVS). +.NH 2 +Echelle +.PP +As some of the algorithms were developed the Echelle data was brought +to my attention. It is considerably simpler than the MAP data because +it is unblended and of high spectral resolution. +Furthermore, the quartz exposure +can be made wider than the object exposures and flexture is not a +problem. The principle problem in this data was the +prevalence of cosmic rays. It pointed to the need to maintain generality +in dealing with both the MAP data and other types of +multi-spectra data which have different profiles, may or may not be +merged, and may or may not have different widths in quartz and object. +Dealing with the cosmic ray problem lead to a very effective solution +usable in both the Echelle and multi-aperture plate data. +.NH +User Level Extraction Logic +.PP +The user should generally only be concerned with the logical steps of +extracting the individual spectra from the multi-spectra image. This +means that apart from specifying the detector system and the format +he should not deal with details of the detector and the format. +In the paper, +\fIDetailed Specifications for the Multi-Spectra Extraction Package\fB, +the \fBIRAF\fR extraction package design and program specifications +are described. +.NH +Flat Fields +.PP +There are two types of flat field situations depending on the spectral +resolution. When the resolution is high then the spectral signature of +the continum calibration source, a quartz exposure, will be unimportant +and variations in the signal will be due to detector sensitivity variations. +In this case the quartz frame, or average of several frames, is the flat +field and division of the object frames by the quartz frame is all that +is required. However, a special +image division program is desirable to handle the region of low or absent +signal between the the spectra. This is described in section 4.2. +.PP +In the alternate case of lower resolution the quartz spectral signature is +larger than the detector response variations. A flat +field in which the intrinsic quartz spectrum is removed is produced by +assuming that the true value of a pixel is given by the smoothed average +of the pixels near that point in position and wavelength and taking +the ratio of the data value to the smoothed value. +This requires a special smoothing program described in section 4.1. +After the flat field is generated then the same image division +program used for the Echelle data is applied. +The image division and smoothing programs are general image operators and +not specific to the Multi-Spectra Extraction Package. +.NH 2 +MULTISPEC_FLAT +.PP +The multi-aperture plate data varies in both dimensions. Thus, any averaging +to smooth the image must take this variation into account. In the Cyber +a flat field for the multi-aperture plate data smooths across the dispersion +by modeling the spectra. This is a difficult task to do accurately because +the true shape of the spectra is not known and the counts vary greatly +and rapidly in this dimension. This approach has the further difficulty +that it is not possible to average several quartz exposures directly. +.PP +The alternate approach to modeling is statistical averaging. +Averaging across the dispersion requires very high order polynomials +because of the rapid variations; +the spectra are typically spaced about 8 pixels apart and there are on +the order of 50 spectra. On the other hand, variations along the dispersion +are much slower even if the spectra are slightly skewed; a bad case would +have two peaks in 800 pixels along the y dimension. This kind +of variation is tractable with relatively simple averaging polynomials +and is the one adopted for the multi-aperture plate data. +.PP +The flat fields are produced by a quadratic moving average along the +y direction. This means that the region centered at a given pixel +is fitted by a least-squares quadratic polynomial and the value of the +polynomial at that point is the appropriate statistical average. +The width of the moving average is an adjustable parameter. +At the edges of the frame where it is not possible to center a region of +the specified width about the desired pixel the polynomial fit is used to +extrapolate the average value to the edge. +.PP +Because the quadratic fit will +be influenced by bad pixels an attempt is made to detect and smooth over +the bad pixels. This is accomplished by comparing the smoothed values to +the observed values and ignoring pixels with a value of + +.EQ (1) + chi = | observed - smoothed | / sqrt smoothed +.EN + +greater than a specified value. Then the smoothing is recalculated and tested +again for bad pixels. This iteration continues until no further pixels +are rejected. +.PP +Following the smoothing the flat field is produced by ratioing the raw quartz +to the smoothed quartz. Pixels of low signal (specified by the +parameter \fIconstant\fR ) +are treated by the equation + +.EQ + r = (data + (constant - smoothed) ) / constant . +.EN + +The resultant flat field image is then divided into the object frames in +the manner described in the next section. +.PP +Experience with data from the Cryogenic Camera has proved very good. +The flat field which is produced can be examined on a display. It +shows fringing at red wavelengths and is not too strongly affected +by bad pixels. Some further effort, however, could go into smoothing +over the bad pixels. +.PP +The smoothing operation on data from the Cryogenic Camera actually +consists of four steps. The quartz exposures are first averaged. +The average quartz is rotated so that the dispersion +direction is the most rapidly varying or x dimension. Then the +smoothing is performed along x followed by another rotation to return +the flat field image to its original orientation. The reason for the +rotations is that they can be done quickly and efficiently whereas +smoothing along the y dimension is very slow and coding an efficient +version is much more complicated than doing a single line at a time. +.NH 2 +FLAT_DIVIDE +.PP +The Echelle data has quartz frames which can be used directly as flat fields. +One just has to divide the object frames by the quartz or average of several +quartz. However, in image division consideration has to be given the +problem of division by zero or very small numbers. In direct imaging this +may not be much of a problem but in multi-spectra data the region between +the spectra and near the edges of the spectra will have very low counts. +Another aspect of image division for making flat field corrections is the +scaling of the result. The flat field integer image data must be large +to give accurate relative response values. However, one wants to divide +an object frame by values near unity. +This section describes a special image division operator allowing the user +to specify how to handle these cases. +.PP +The parameters are a \fIdivision threshold\fR +(default of zero) and a \fIthreshold violation value\fR. Values of the +denominator above the \fIthreshold\fR are treated separatedly from those +below the \fIthreshold\fR. The denominator image is scaled to have an +average of one for pixels above the \fIthreshold\fR. The pixel by pixel +division is then performed for those points for which the denominator +is above the \fIthreshold\fR. Pixels for which the denominator is below the +\fIthreshold\fR are set to the \fIthreshold violation value\fR in the resultant +image if the \fIviolation value\fR is specified. If the value is not +specified then the numerator value is taken as the resultant value. +The divisions can be done in place or the result put into a new image file. +.PP +For the multi-spectra situation where the object spectra have a +smaller width than the quartz, as in the Echelle, one can either +set the \fIthreshold +violation value\fR to zero or not set it at all resulting in either +exactly zero or background values between the spectra while still flattening +the spectra. This allows looking at the flattened spectra without the +annoying "grass" between the spectra caused by dividing by small +values. +.NH +Extraction +.NH 2 +MULTIAP_EXTRACT +.PP +The extraction of spectra from multi-aperture plate images consists of +a series of steps. The steps are executed from a script. +The command + +.FX +ms> multiap_extract "ap165.*", "", 165, 50 +.EX + +will take the flattened images, ap165.*, from aperture plate 165 with 50 +spectra and automatically locate the spectra, model the profiles, and +extract the one dimensional spectra. The script consists of +a number of steps as described below. +.PP +\fBFind_spectra\fR (section 6) initializes the \fBmultispec\fR data file +and does a peak search to determine the initial positions of the +spectra. +\fBFind_bckgrnd\fR fits a polynomial of order 1 (or more) for the pixels which +are not near the spectra as defined by \fBfind_spectra\fR. +.PP +The spectra are then modeled in bands of 32 lines by the model profiles +described in section 8.1. The first \fBmodel_fit\fR uses three Gaussian +parameters for +each spectra measuring the peak intensity, peak position, and width. +The second \fBmodel_fit\fR adds a fourth parameter to modify the wings of the +profile. +.PP +The \fBmodel_extract\fR program extracts the spectra line by line and also +detects and removes cosmic ray events which do not fit the model +profiles (see section 9). +In outline, the extraction of blended spectral data uses the +model profiles to determine the fraction of light +from each of the neighboring spectra at the pixel in question. The +appropriate fraction of the +.ul +observed +pixel intensity (minus the background) is +assigned to the luminosities of the spectra. There are two versions +of the \fBmodel_extract\fR extraction. The first simultaneously fits the +peak intensity of all the spectra and the second uses the +data value at the peak of each spectra to normalize the model. The +first method is slow and accurate and the second is fast and approximate. +Because the models are used in extraction only to define the relative +contributions of neighboring spectra to the total observed pixel luminosity +the speed of the approximate method far outweighs the need for +accuracy. However, cleaning the spectra of cosmic rays is a different +matter and is discussed further in section 12. +.PP +After the extraction the spectra are correlated with the aperture plate +description using \fBap_plate\fR (see Section 10) to determine the +relative wavelength offsets and assign identification information to +the spectra. +.PP +For successive frames it is not necessary to resort to the initial +steps of finding the spectra and fitting from scratch. The \fBcopy_params\fR +routine makes a new copy of the fitting database. Small shifts in positions +of the spectra and the peak intensities are determined by doing a two +parameter fit for the peak and position using the previously determined +shape parameters. +Changes in the shape of the spectra are then determined by the three +and four parmater fits. Because the solution is likely to be close to +the previously determined shape the transfering of one solution from a +previously solved image is faster than starting from scratch. +Note that the shapes as well as the positions and peak intensities +of all frames including the object exposures are allowed to change. +.PP +The spectra are then extracted from the image by \fBmodel_extract\fR and the +process repeats for the succeeding images. +.PP +One useful feature is the ability to specify the bands or lines to be +modeled or extracted. +This feature is useful for diagnosising the programs quickly. +The default is all bands or lines. +.NH 2 +ECHELLE_EXTRACT +.PP +The extraction of the unblended Echelle spectra is performed +begins in a similar way with \fBfind_spectra\fR and \fBfind_bckgrnd\fR. +The extraction and cleaning, however, uses \fBstrip_extract\fR which +adds up the instrumental counts for each unblended spectra at each +wavelength to get the total luminosity. +.NH +FIND_SPECTRA -- Finding the Spectra +.PP +The first step in the extraction and processing of multi-spectra data is +to locate the spectra. This can be done interactively by +the user but it is far preferable to automate the process; +particularly since multi-spectra data can have a large number of +spectra and frames. The approach is to find the peaks in a line, or +average of lines, sort the peaks found in some manner, such as by +strength, and select the expected number of peaks from the top of the +list. +Beyond this simple outline are several algorithmic details such as how +to define and locate valid peaks and how to sort the list of peaks. +Peak finding is a general problem and a subroutine for peak finding is +described below. The \fBfind_spectra\fR program provides an +interface between the \fBmultispec\fR data file and the +general peak finding algorithm. +.PP +The \fBpeaks\fR function takes arrays of x (position) and y (value) points +and the number of +points in the arrays and returns the number of peaks found. It also +returns the estimated positions of the peaks in the x array and the +extimated peak values in the y array in order of peak likelihood. +There is one user parameter, the smoothing \fIwidth\fR. +The choice of the \fIwidth\fR parameter is dicatated by how closely and how +wide the peaks are to be expected. +The algorithm takes a region of \fIwidth\fR points +centered on each x point and fits a quadratic; + +.EQ +y sub fit = a + b x + c x sup 2~. +.EN + +A peak is defined +when the slopes, $b sub 1$ and $b sub 2$, of two neighboring points +$x sub 1$ and $x sub 2$ change +sign from positive to negative and the curvatures, $c sub 1$ and $c +sub 2$, are less than -0.001 for both points. +The quadratic polynomials define two estimated peak positions + +.EQ +x sub 1 sub peak = x sub 1 - b sub 1 / (2 * c sub 1 ),~~ +x sub 2 sub peak = x sub 2 - b sub 2 / (2 * c sub 2 )~. +.EN + +The offsets are then normalized to give a linear interpolation +fraction +$f = ( x sub 1 sub peak - x sub 1 ) / ( x sub 2 sub peak - x sub 1 sub +peak )$ in the interval between the two points. +The estimated position of the peak is then + +.EQ +x sub peak = f * ( x sub 1 - x sub 2 ) +.EN + +and the estimated peak value is the average value of the two quadratic +polynomials at $x sub peak$. The curvature at the peak is +estimated by $c sub peak = c sub 1 + f * (c sub 1 - c sub 2 )$. +Finally, the peaks are sorted by the magnitude of the peak curvature. +.PP +This peak finding algorithm works quite well. I have also used it to +automatically locate peaks in the extracted one dimensional spectra +and then do peak correlations between spectra to find a relative +wavelength solution. Some such use of this program may be implemented +in either future additions to the Multi-Spectra Extraction Package or +the Spectral Reduction Package. +.PP +In \fBfind_spectra\fR the number of spectra to be found is specified by +the user. The user should have previously looked at an image +on a display or done a profile plot across the +dispersion to count the observed spectra. +Additional parameters specify the columns in which the spectra +are to be found and the minimum separation and width of the spectra. +The column specification allows the elimination of problems with defective +areas of the detector such as the LED in the Cryogenic Camera. The minimum +width and separation provide algorithmic constraints to the spectra finding +procedure. +.PP +The peaks are found at two or more points in the +multi-spectra image for a band of 32 lines using a +\fBpeaks\fR \fIwidth\fR parameter of 5. After the peaks are found +at a number of bands in the image a linear fit is made to determine any small +slope of the spectra relative to the columns. +The reason for specifying only a few bands is that the process of +finding the peaks is moderately slow and only two bands are needed for +determining the initial position angle of the spectra to the y +dimension. Furthermore, some bands do not give a satisfactory result +because of extraneous data such as the LED in the CRYOCAM or bad +focus. Another possibility is that a spectrum may go off the edge +and in order to "find" the spectrum for partial extraction bands that +include the on image part of the spectrum must be specified. +.NH +FIND_BCKGRND -- Background +.PP +The background on a multi-spectra image is the result of very broad +scattering as opposed to the narrower scattering which produces +distinguishable wings on individual spectra. +Modeling of the background in a Cryogenic Camera multi-aperture plate +image shows that the background is well explained by a broad +scattering function. +It is not reasonable, however, to model the scattering to this detail +in actual extractions. +Instead a smooth polynomial is fitted to the pixels not covered by +spectra. The order of the polynomial is a specified parameter. +For the CRYOCAM MAP data a quadratic is appropriate. +.PP +The algorithm is the same for all multi-spectra data except for the +choice of parameters. First, the location of the spectra must be +determined. This is done by the \fBfind_spectra\fR program. There +are two principle parameters, a buffer distance and the order of the +fitting polynomial. Each line, or average of several lines, is fitted +by least-squares for the points lying farther than the buffer +distance from any spectra. If there are no points which completely +stradle the spectra, i.e. located on each side of the spectra, then +the order of the fitting polynomial is ignored and a constant, or +first order polynomial, is determined. +A hidden parameter specifying the columns allowed for searching for +background points is available so that bad parts of the image can be +ignored. +.PP +A difference in philosophy with the Cyber programs is that the +determined background is not subtracted from the image data. It is +instead kept in the database for the image. Generally, it is better to +modify the basic image data as little as possible. This is the approach +used in the Multi-Spectra Extraction Package. +.NH +Spectra Profiles +.NH 2 +MODEL_FIT -- Models for Multi-Spectra Images +.PP +The object of modeling is to separate blended spectra for extraction +and to remove artifacts, such as cosmic rays, which do not fit +the model. The models should have the minimum number of parameters +which give residuals approaching the detector statistics, they +should incorporate constraints based on the physics of the +detector/camera system, and the models must be ammenable to a +statistical fitting algorithm which is stable. +There are a large number of possibilities. +.PP +An important point to bear in mind during the following discussion is +the necessary accuracy of the model fitting. In the design proposed +here the model fitting is not used for determining the smooth quartz. +Use of a model for making a flat field would require a very accurate +model and using an average quartz is not possible. However, for +extraction the model is used only to indicate the +relative fraction of light for each spectrum when the spectra are +blended. The cleaning application is more critical but not nearly so +much as in the flat field modeling. Thus, though we do a good job of +model fitting (better the the Cyber modeling) some additional features +such as smoothing along the spectra are not included. +Also, though some improvement can be gained by the additional shape +parameters in the fit, they are not necessary for the required purpose +and can be left out leading to a faster extraction. +.PP +During the course of my investigation I tried more than one hundred +models and combinations of constraints. Some general results of this +study follow. +The model which I find gives the best results has six parameters not +counting the background. The model is defined by the following +equations where x is the cross dimension. + +.EQ (1) +I = I sub 0 exp (-s * ( DELTA x sup 2 )) +.EN +.EQ +DELTA x = (x - x sub 0 ) +.EN +.EQ +s = s sub 0 + s sub 1 y sup 2 + s sub 2 y sup 3 +.EN +.EQ +y = DELTA x / sqrt { DELTA x sup 2 + x sub 1 sup 2 } +.EN + +The model consists of a intensity scale parameter, $I sub 0$, +and a profile which is +written in a Gaussian form. The center of the profile is given by +the parameter $x sub 0$. The profile is not exactly Gaussian because the +scale, $s$, is not a constant but depends on $DELTA x$. The scale +function has three terms; a constant term, $s sub 0$, which is the scale +near the center of the profile, and even and odd terms, $s sub 1$ +and $s sub 2$, +which change the scale in the wings of the profile. +.PP +The characteristic of the profile which must be satisfied is that at +large distances from the profile center the scale is positive. This +requirement means that the profile will be monotonically decreasing at +large distances and will have a finite luminosity. This point was +crucial in determining the form of the scale function. A straight +power series in $DELTA x$ does not work because power series diverge. +Instead, the scale function is defined in terms of a separation +variable $y$ which is bounded by -1 and 1 at infinite separation and is +zero at zero separation. The parameter $x sub 1$ defines a characteristic +distance where the character of $y$ changes from going as $DELTA x$ to +asymptotic to 1. The parameters are, thus, $I sub 0$, $x sub 0$, $s sub 0$, +$s sub 1$, $s sub 2$, $x sub 1$. +.PP +An important property of this model is that the terms have a physical +interpretation. The profile scale and profile center are obvious and +any model must include them. It is the remaining terms, $s sub 0$, +$s sub 1$, $s sub 2$, +and $x sub 1$, which are called the shape parameters, which are interesting. +In an ideal aperture plate system the shape of a profile would be +given by the projection of the circular aperture into the cross dimension: + +.EQ +P( DELTA x ) = sqrt {1 - a DELTA x sup 2} +.EN + +where the constant a is related to the size of the hole by + +.EQ +a = 1 / r sup 2 +.EN + +For small $DELTA x$ the profile can be expressed in the Gaussian form with +a scale + +.EQ +s = a( 1/2 + a DELTA x sup 2 + ...) +.EN + +Thus, even in a perfect aperture plate system a Gaussian form shows the +scale increasing from a central value determined by the size of the hole. +In other words, the profile decreases more sharply than a Gaussian. +.PP +However, no aperture plate system is ideal because the thickness of +the aperture plate is finite and there is scattering and changes in +the focus of the system. One must +convolve the profile above with a scattering/focus function. One can show +that for reasonable functions, exponential and Gaussian, +with some scale b the final profile is a function of the ratio b/a. +If the ratio is less than 1 then the profile will be more like that of +the hole and the profile will be sharper than a Gaussian in the wings. +If the ratio is much greater than 1 then the profile will become the +scattering profile at large separations. Simulations using Gaussian +and exponential scattering profiles show behaviors very much like the +profile (1) with $s sub 1$ greater than zero when b/a < 1 +meaning the profile becomes sharper (than a Gaussian) in the wings +and $s sub 1$ < 0 when b/a > 1. +Thus, $s sub 1$ defines the scale of the scattering profile relative +to the hole size. +The size of the hole is incorporated into the parameter $x sub 1$. +The parameter $s sub 2$ allows an asymmetry in the profile. +.PP +An interesting property of the scale function is that it is all right +for it to be negative at small distances from the profile center. This +occurs when $s sub 0$ is negative. The effect of this, provided $s$ +becomes positive at large distances, is to give a two horned profile. +This, in fact, is observed when the focus of the system becomes very +poor. +.PP +The best fits (least chi-square or rms residual) are +obtained when each spectrum at each wavelength has independent +parameters. However, this sometimes gives rise to unphysical results. +If left entirely unconstrained the parameter fitting algorithm can +make one line broad and dominant and a neighboring line weak and +sharp. +This is not, of course, a property of the camera or detector. +Thus, constraints based on the physics of the +camera/detector are used. This means that the shape +parameters $s sub 0$, $s sub 1$, $s sub 2$, and $x sub 1$ +are coupled locally by making them vary as a polynomial of position +across the dispersion. One might also +constrain the variation of the shape along the spectrum as is done in +the Cyber. This is not needed because there are no drastic differences +between the fitted parameters at neighboring points along the spectra. +.PP +My experience with the Cyrogenic Camera system has shown the +following. The focus ripples twice across the CCD with the +propagation angle being approximately 30 degrees from the long dimension. +The change in focus is partly just a scale change. This is seen in +the correlation of $s sub 0$ with the image scale found by \fBap_plate\fR. +The shape parameter $s sub 1$ changes sign from positive to +negative indicating that when the focus is good the profile +decreases faster than a Gaussian and when the focus is bad it decreases +slower. Occassionally the focus is very bad and $s sub +0$ is negative and $s sub 1$ is small and positive causing a broad two +horned profile. The +assymetry parameter, $s sub 2$, is useful only when the signal is strong near +the peak of a quartz exposure. It is not really necessary to include +it in the model fits. The assymetry parameter was dominant, however, +in some Echelle data which were clearly asymmetric. The value of +$x sub 1$ is +not highly sensitive and can be fixed for a given hole size. Large +changes in the hole size would require resetting $x sub 1$. +The use of the four parameters, $I sub 0$, $x sub 0$, $s sub 0$, +and $s sub 1$, allow good fits +to all the data I've examined including those in which the peak/valley +intensity ratio across the spectra is about 1.1. It is the importance +of the parameter $s sub 1$ which improves the fitting dramatically over the +Cyber three parameter fitting (in addition to a different fitting +algorithm). +.PP +The heart of profile fitting is the solution of the multi-parameter +least-squares problem. In a blended multi-spectra image the profile +parameters of one spectra are affected by its neighbors which are, +in turn, affected by their neighbors and so forth. The key to this +type of problem is to realize that only nearest neighbors affect the +profile parameters and this leads to a "banded" least-squares matrix. +A banded matrix is one in which cross terms far from the diagonal are +zero. Solution of banded matrices is much more efficient than solving +the entire matrix. This allows solution for more than 100 parameters +simultaneously in a short time. +.PP +Use of the banded multi-parameter solution has the restriction, however, +that there can be no parameters in the model which are not local to +the profiles. This affects the way +global constraints are applied to the parameters. In particular, +the way the shape parameters are constrained to vary smoothly across the +detector. +The shape parameters are first found as independent parameters by the +banded matrix solution and then smoothed by a polynomial in x. +.PP +An area which was extensively investigated was the appropriate +weighting to use for the model fitting. The most likely choices are +weighting by $1 / sqrt data$ and unit weight corresponding to +$chi sup 2$ +and least squares fitting. It was found that the two methods +agreed fairly closely but that the least squares fitting was more +appropriate because the blending correction depends largely on the +value of the peak intensity and less on the exact fit of the wings. +With $chi sup 2$ the peak is fit with less accuracy in order to improve +the fit in the wings of the profile. In some cases this gave clear +errors in estimating the peak intensity and, hence, the proper contributions +between the blended spectra were not made. +.PP +Now follows the details of the fitting algorithm. +The algorithm is a series of script steps in \fBmultiap_extract\fR +which call the model fitting program \fBmodel_fit\fR with different +parameters. In the script all bands are fit, $x sub 1$ is fixed, +and the asymmetry shape parameter $s sub 2$ is ignored. +The four parameter fit is applied to bands of 32 lines. The band +solutions are linearly interpolated to the full image and then only +the intensity scale parameter is calculated for each line during the +extraction of the spectra with \fBmodel_extract\fR. +.PP +The general fitting scheme proceeds as follows: +.LP +1. Fit the three parameters $I sub 0$, $x sub 0$, $s sub 0$ with +$x sub 1$ fixed and $s sub 1$ and $s sub 2$ +zero. This is precisely a Gaussian fit. The three parameters are +determined simultaneously for all the lines at once using the banded +matrix method. Thus for 50 lines the solution has 150 variables. +After each fit the scale +parameter $s sub 0$ is smoothed by a polynomial in x. The polynomial is +taken with seven terms. +.LP +2. Once the improvement in each iteration becomes less than a +specified amount (2% in rms residual) the next parameter $s sub 1$ is added. +The solution has two steps: fit for $s sub 0$ and $s sub 1$ with $I sub 0$ +and $x sub 0$ fixed and +then fit $I sub 0$ and $x sub 0$ with $s sub 0$ and $s sub 1$ fixed. As before the scale terms +are smoothed by a seventh order polynomial. Attempts to solve for all +four parameters a once gave unstable results for reasons I don't +understand. +.LP +3. If desired, the last shape parameter $s sub 2$ can be added by solving +for $s sub 0$, $s sub 1$, and $s sub 2$ while holding $I sub 0$ and +$x sub 0$ fixed and then solving for +$I sub 0$ and $x sub 0$. This provides some improvement when the signal is very +strong but is generally not needed in the multi-aperture plate data. +It can be an important term in the Echelle data. +.LP +4. It is possible to then adjust $x sub 1$ followed by steps 2 or 3. +However, this gives very little improvement and $x sub 1$ should be fixed for +each type of data. +.LP +5. During the final extraction when individual lines are evaluated a one +parameter fit is used to find $I sub 0$ for each spectra. This is +rather slow, however, on the order of 3 hours per frame. By using +the pixel value near $x sub 0$ as the value for $I sub 0$ the extraction is reduced +to 13 minutes per frame (see section 12). +.PP +In addition to the preceeding steps the fitting algorithm applies some +heuristic constraints. These constraints limit how far the peak positions can +shift in each iteration, require the peak intensity to remain positive, and +limit the scale function to be positive at large values of y. +.NH 2 +STRIP_EXTRACT -- Unblended Profiles +.PP +For unblended multi-spectra data the profiles can be anything. The profiles +are obtained by averaging a number of lines (say 32) and normalizing +at some point like the peak value. These profiles are then used for +detecting bad pixels, such as cosmic rays, and correcting for them as +discussed in the section on cleaning. Modeling using the \fBmodel_fit\fR +program is only used on Echelle data to find peak positions +accurately in order to follow any curvature of the spectra. +.NH +Extraction and Cleaning +.PP +The extraction of spectra are done separately from the modeling. It is +possible to extract spectra without any modeling at all using +\fBstrip_extract\fR. The extraction step also allows the user to specify +if cleaning of the spectra for cosmic rays is desired. Also modifying +the image is an option. +.NH 2 +MODEL_EXTRACT +.PP +Extraction and cleaning using a model fit is described here. +First the $I sub 0$ values for the model profiles are determined for +all the spectra in a line either by multi-parameter fitting or by +taking the peak value. The pixel values are then compared to the +model in a chi-squared way: + +.EQ +r = (data - model) / sqrt model +.EN + +If the value of r is larger than a set amount, say 5, then the pixel +value is set to that of the model. Since the "bad" pixel may affect +the intensity scale $I sub 0$ the cleaning is iterated until no further +pixels are changed. +.PP +The fitting of the data from an individual line of data to the model profiles +is the key element in this scheme. The best method is to use all the +data in a multi-parameter least square fit. This minimizes the effect +of bad pixels on the estimated profile which is the essence of this +cleaning method. However, while the time required to do this for one +line is not great, it adds up to nearly three hours for the 800 lines +in a CRYOCAM frame. A quick alternative is to scale the model profile +by the data value at the peak position. This is +quite fast. However, if the peak has a cosmic ray event or is +otherwise bad then the estimated profile will not correspond closely +to the data profile and the cleaning procedure will make gross errors. +The limited experience I've had with the Echelle and MAP data +has worked well with using the peak estimate. However, the possible +problems make me nervous and some compromise based on using more than +the peak to estimate the intensity scale of the profile to the data +needs to be found. This is important because much of the feedback on +the \fBmultispec\fR package from Paul Hintzen and Caty Pilachowski +have dealt with +the particular usefulness of a good cosmic ray cleaning algorithm in +extracting multi-spectra data. +.NH 2 +STRIP_EXTRACT +.PP +Removing cosmic rays is the major part of Echelle extraction. +Because these are unblended spectra of arbitrary shape a strip +extraction is all that is needed. +The cleaning is done by the same algorithm used for the multi-aperture +plates except that the profiles are found, as described earlier, by +averaging a number of lines. +The intensity scaling is determined from either a least-square fit +or the peak value. +The same question about the appropriate way to +determine the fit of the profiles to the data discussed previously +applies here except since the spectra are not blended the spectra +can be treated separately in any least square fitting. +.NH +AP_PLATE -- Aperture Plate Correlation +.PP +The final step in the extraction of a multi-aperture plate image is to +correlate the spectra with the on-line database description of the +drilled hole positions. This allows for estimates of relative wavelength +offsets and the identification of the spectra with the ID, RA, and DEC +parameters. +.PP +The spectra are fitted to the known aperture plate drilled positions, given in +millimeters, to find an \fIangle\fR for the aperture plate relative to the +detector x dimension and the image \fIscale\fR in pixels / millimeter, + +.EQ +x sub fit = a + scale (x sub drill cos (angle) + y sub drill sin (angle))~. +.EN + +If the number of spectra is less than that given by the aperture plate drilled +positions then a correlation is done leaving out sequences of +consecutive holes until the fit residual is minimized. If the number of +spectra is greater than that supposedly drilled then sequences of +consecutive peaks are left out of the fit to minimize the residual. +The missing holes or extra peaks are printed out and, if allowed, the aperture +plate description file is modified, otherwise the program terminates. +In all cases if the final fit residual is greater than 1 +pixel the program will terminate. +The program prints out the \fIangle\fR of the aperture plate and the \fIscale\fR +which is also stored in the database. +.PP +An indication that a large part of the focus variation is purely a +scale change is that the derived image \fIscale\fR correlates very well with +the width of the spectra as derived from the profile fitting. I +estimate that at least 50% of the focus variation is a scale +variation. This is good news in the sense that a scale variation will +be taken out in the dispersion solution and lines in different parts +of the detector will become more similiar after the solution. +However, the scale variation does not account for all the profile +shape changes and there is indeed a change in the point spread function +across the detector. +.NH +Problems +.PP +There a few problems which I have not been able to resolve or have not +had the time to consider. The problems which are largely intractable +without a great deal of effort are the unexplained background +variations and deconvolving the spectra for the variation in the +point-spread-function. The background variations are abrupt increases +in the background in parts of the CRYOCAM detector. The step edge sometimes +occurs under the spectra and so any smooth polynomial fit to the +background will not be very good. The modeling of the multi-aperture +plate profiles provides information about the point-spread function +but a deconvolution of the variation in the PSF is a difficult problem +and not warrented at this time. +.PP +I had expected that the large scale response of the CRYOCAM could be +corrected by determining an overall average quartz spectrum from all the +extracted quartz spectra and then dividing the object spectra in each +hole by the ratio of the average quartz spectra from that hole to the +overall average quartz spectrum. This was attempted and it was found +to work only partially. Specifically, while there might be a 20% +difference between a spectrum on the edge and one near the center of +the detector the quartz correction left a 10% difference in the object +spectra. This is apparently due to a poor illumination by the quartz +light source which does not correspond to the telescope illumination. +This quartz correction technique may be made available to users if +desired. +.NH +Comparison with the Cyber Extraction Package +.PP +The discussion of this section must be qualified by the fact that I +have not used the Cyber Extraction Package and I base my understanding on the +algorithms from the Multi-Aperture Plate Data Reduction Manual and +conversations with knowledgable people. There are many differences +both major and minor and this section only seeks to mention the +some of the important differences. In the Cyber package: + +The background is subtracted from the images as a preliminary process. + +The background is either constant or linear across the spectra. + +The flat fields are produced by modeling the quartz and data from +several quartz exposures cannot be easily combined. + +The initial peak finding and aperture plate correlation algorithm is less +automated in determining missing or additional holes. + +The model fitting uses only a three parameter Gaussian model +and the algorithms do not yield results when the focus becomes poor. + +The fitting algorithm is neighbor subtraction rather than full +simultaneous solution for all the profiles. + +The model fitting is applied only to a quartz and the model is transfered to +object exposures. This does not allow the shape of the profiles to +change with time as the telescope moves. + +The modelling does not couple solutions for neighboring spectra +across the dispersion as is suggested in this design and it does smooth +along the spectra which is not done in this proposal. + +The extraction is only to some specified sigma in the model profile and +there is no attempt to correct for blending. + +There is no cleaning of the spectra. +.NH +Discussion +.PP +The only data which has passed beyond the extraction phase using the +algorithms described here was that of Paul Hintzen. +Comparison of data reduced by the TV package for +spectra extracted by both the Cyber package and the techniques of the +suggested \fBmultispec\fR package were quite comparable. To the level he +examined +the spectra there was no clear increase in accuracy though the \fBmultispec\fR +extractions generally had higher counts due to the full extraction of +the blended spectra. The big advantages found were +the ability to extract all the data even when the focus +became very poor and the success of the cosmic ray cleaning +algorithm. Thus, Hintzen feels that the need for speed in the extraction +(primarily dependent on the cleaning algorithm) +is modified significantly by the difficulty of dealing with cosmic +rays in the TV spectral analysis programs. More exploration +of techniques for determining the profile intensity scale from the +model without the full multi-parameter solution is warrented for this +reason. +.PP +I have extracted some Echelle data including field flattening. The +data had a considerable number of cosmic rays which were removed +quite well. The extracted spectra were put into a CAMERA format +for further analysis. +.PP +The programs were recently applied to a long slit analysis problem +being studied by Vesa Junkkarinen. The image was already flat fielded. +The data had two closely spaced and very faint diffuse objects and scattered +light from a nearby QSO. +The three spectra were so weak and closely spaced +that the automatic finding was not used. However, the rest of the modeling +and extraction were applied directly. +The \fBfind_bckgrnd\fR program, whose original purpose was to correct for +scattered light, worked well to extrapolate the sky across the +image. The model fitting accurately followed +the peaks of the spectra but the shape fitting was only moderately accurate +since the model shape parameters are not suited to modeling galaxies. +It successfully extracted spectra with a minimum of effort on my part. +Analysis of the extracted spectra and comparison with other techniques +must still be done. The conclusions to be drawn from this experiment are +that with sufficiently general multi-spectra tools multiple objects in +long slit spectroscopy can be handled. +.PP +One area in which I do not have practical experience is +the extraction of HGVS data. I believe +the proposed design will work on this type of data. +.PP +A point which needs to be considered in the final design are the +formats of the data files. The currently used one dimensional spectra +formats are an IIDS format and a CAMERA image format. +The formating of data files for the current spectral analysis packages by +\fBto_iids\fR starts from the \fBmultispec\fR database and throws away a lot +of information about the spectra. +Some refinement of this database should focus on the format +to be used by a new \fBIRAF\fR spectral analysis package. +.PP +It should be pointed out that many of the operations can have +alternate algorithms substituted. In particular, the smoothing +algorithm for the multi-aperture plate flat fields can be replaced by +some other scheme. The links between the multi-parameter fitting +program and the model have been made very general for investigating +a broad range of models. Thus, it is also possible to substitute +additional model profiles with relative ease. +.PP +Estimates of excution time are taken from the experimental C programs +implementing the algorithms of this design and they are only +approximate estimates. The steps corresponding +to \fBdebias\fR, \fBmultispec_flat\fR, and \fBflat_divide\fR for +the multi-aperture data from the CRYOCAM take +about 1 hour for a typical set of frames, say 5 to 15. This includes +debiasing, triming, computing a flat field from several quartz frames +and dividing the quartz into the object frames. +.PP +The CRYOCAM \fBmultiap_extract\fR phase takes about 40 minutes for the modeling of a frame using 32 lines per band and either 3 hours for an extraction +using the profile fitting +method or 14 minutes for extraction using the peak profile scaling +method. +.PP +Finally, the \fBto_iids\fR takes about 3 minutes per frame. It takes +this long because it has to convert the \fBmultispec\fR database organized across +the dispersion into formats in which the data is stored as consecutive +spectra; i.e. a type of rotation operation. |