diff options
Diffstat (limited to 'noao/twodspec/multispec/doc')
| -rw-r--r-- | noao/twodspec/multispec/doc/MSalgo.ms | 1032 | ||||
| -rw-r--r-- | noao/twodspec/multispec/doc/MSalgo_c.doc | 522 | ||||
| -rw-r--r-- | noao/twodspec/multispec/doc/MSalgo_c.hlp | 449 | ||||
| -rw-r--r-- | noao/twodspec/multispec/doc/MSspecs.doc | 698 | ||||
| -rw-r--r-- | noao/twodspec/multispec/doc/MSspecs.hlp | 659 | ||||
| -rw-r--r-- | noao/twodspec/multispec/doc/MSspecs_c.hlp | 243 | ||||
| -rw-r--r-- | noao/twodspec/multispec/doc/findpeaks.hlp | 88 | ||||
| -rw-r--r-- | noao/twodspec/multispec/doc/fitfunc.hlp | 73 | ||||
| -rw-r--r-- | noao/twodspec/multispec/doc/fitgauss5.hlp | 148 | ||||
| -rw-r--r-- | noao/twodspec/multispec/doc/modellist.hlp | 52 | ||||
| -rw-r--r-- | noao/twodspec/multispec/doc/msextract.hlp | 172 | ||||
| -rw-r--r-- | noao/twodspec/multispec/doc/mslist.hlp | 77 | ||||
| -rw-r--r-- | noao/twodspec/multispec/doc/msplot.hlp | 44 | ||||
| -rw-r--r-- | noao/twodspec/multispec/doc/msset.hlp | 104 | ||||
| -rw-r--r-- | noao/twodspec/multispec/doc/multispec.ms | 532 | ||||
| -rw-r--r-- | noao/twodspec/multispec/doc/newextract.hlp | 61 | ||||
| -rw-r--r-- | noao/twodspec/multispec/doc/newimage.hlp | 130 |
17 files changed, 5084 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. diff --git a/noao/twodspec/multispec/doc/MSalgo_c.doc b/noao/twodspec/multispec/doc/MSalgo_c.doc new file mode 100644 index 00000000..b3322dff --- /dev/null +++ b/noao/twodspec/multispec/doc/MSalgo_c.doc @@ -0,0 +1,522 @@ + MULTISPEC (Dec83) Multispec Algorithms MULTISPEC (Dec83) + + + + Algorithms for the Multi-Spectra Extraction Package + Analysis and Discussion + December 2, 1983 + + + +1. Disclaimer + + This should not be taken as a statement of how the algorithms of +the final package should function; this is merely an analysis and +discussion of the algorithms, and should be followed by further +discussion before we decide what course to follow in the final +package. We may very well decide that the level of effort required to +implement rigorously correct nonlinear fitting algorithms is not +justified by the expected scientific usage of the package. Before we +can decide that, though, we need an accurate estimate of the level of +effort required. + +In attacking nonlinear surface fitting problems it is important to +recognize that almost any techniques can be made to yield a result +without the program crashing. Production of a result (extraction of a +spectrum) does not mean that the algorithm converged, that the +solution is unique, that the model is accurate, or that the +uncertainties in the computed coefficients have been minimized. + + + +2. Multispec Flat (pg. 4) + + This sounds like a classical high pass filter and might be best +implemented via convolution. Using a convolution operator with a +numerical kernel has the advantage that the filter can be easily +modifed by resampling the kernel or by changing the size of the +kernel. It is also quite efficient. The boundary extension feature +of IMIO makes it easy to deal with the problem of the kernel +overlapping the edge of the image in a convolution. Since the +convolution is one-dimensional (the image is only filtered in Y), it +will always be desirable to transpose the image. + +The method used to detect and reject bad pixels (eqn 1) is not correct. +The rejection criteria should be invariant with respect to a scaling +of the pixel values. If the data has gone through very much +processing (i.e., dtoi on photographic data), the relation between +photon counts and pixel value may be linear, but the scale is +unknown. Rejection by comparison of a data value to a "smoothed" +value is more commonly done as follows: + + reject if: abs (observed - smoothed) > (K * sigma) + +where sigma is the noise sigma of the data, generally a function of +the signal. + +It is often desirable in rejection algorithms to be able to specify, + + + -1- + MULTISPEC (Dec83) Multispec Algorithms MULTISPEC (Dec83) + + + +as an option, that all pixels within a specified radius of a bad pixel +be rejected, rather than just the pixel which was detected. This is +only unnecessary if the bad pixels are single pixel events (no +wings). Region rejection makes an iterative rejection scheme converge +faster, as well as rejecting the faint wings of the contaminated +region. + + + +2.1 Dividing by the Flat (pg. 5) + + There is no mention of any need for registering the flat with the +data field. Is it safe to assume that the quartz and the object +frames are precisely registered? What if the user does in fact +average several quartz frames taken over a period of time? (Image +registration is a general problem that is probably best left until +solved in IMAGES). + + + +3. Multiap Extraction (pg. 5-6, 8-13) + + The thing that bothers me most about the modeling and extraction +process is that the high signal to noize quartz information is not +used to full advantage, and the background is not fitted very +accurately. The present algorithms will work well for high signal to +noise data, but will result in large (percentage) errors for faint +spectra. + +Basically, it seems to me that the high signal to noise quartz spectra +should, in many cases, be used to determine the position and shape of +the spectral lines. This is especially attractive since the quartz +and spectra appear to be closely registered. Furthermore, if the +position-shape solution and extraction procedures are separate +procedures, there is nothing to prevent one from applying both to the +object spectum if necessary for some reason (i.e., poor registration, +better signal to noise in the object spectrum in the region of +interest, signal dependent distortions, lack of a quartz image, etc., +would all justify use of the object frame). It should be possible to +model either the quartz or the object frame, and to reuse a model for +more than one extraction. + +Let us divide the process up into two steps, "modeling", and +"extraction" (as it is now). The "calibration frame" may be the +quartz, an averaged quartz, or the object frame. Ideally it will have +a high signal to noise ratio and any errors in the background should +be negligible compared to the signal. + +We do not solve for the background while modeling the calibration +frame; we assume that the background has been fitted by any of a +variety of techniques and a background frame written before the +calibration frame is modeled. A "swath" is the average of several +image lines, where an image line runs across the dispersion, and a + + + -2- + MULTISPEC (Dec83) Multispec Algorithms MULTISPEC (Dec83) + + + +column along the dispersion. + + + +3.1 Modeling + + I would set the thing up to start fitting at any arbitrary swath, +rather than the first swath, because it not much harder, and there is +no guarantee that the calibration frame will have adequate signal to +noise in the first swath (indeed often the lowest signal to noise will +be found there). We define the "center" swath as the first swath to +be fitted, corresponding to the highest signal to noise region of the +calibration frame. By default the center swath should be the swath +used by find_spectra, especially if there is significant curvature in +the spectra. + +algorithm model_calibration_frame + +begin + extract center swath + initialize coeff using centers from find_spectra + model center swath (nonlinear) + + for (successive swaths upward to top of frame) { + extract swath + initialize coeff to values from last fit + model swath (nonlinear) + save coeff in datafile + } + + set last-fit coeff to values for center swath + for (successive swaths downward to bottom of frame) { + extract swath + initialize coeff to values from last fit + model swath (nonlinear) + save coeff in datafile + } + + smooth model coeff (excluding intensity) along the dispersion + [high freq variations in spectra center and shape from line] + [to line are nonphysical] + variance of a coeff at line-Y from the smoothed model value is + a measure of the uncertainty in that coeff. +end + + +I would have the background fitting routine write as output a +background frame, the name of which would be saved in the datafile, +rather than saving the coeff of the bkg fit in the datafile. The +background frame may then be produced by any of a number of +techniques; storing the coefficients of the bkg fit in the datafile +limits the technique used to a particular model. For similar reasons, +the standard bkg fitting routine should be broken up into a module + + + -3- + MULTISPEC (Dec83) Multispec Algorithms MULTISPEC (Dec83) + + + +which determines the region to be fitted, and a module which fits the +bkg pixels and writes the bkg image. + +For example, if the default background fitting routine is a line by +line routine, the output frame could be smoothed to remove the +(nonphysical) fluctuations in the background from line to line. A +true two dimensional background fitting routine may be added later +without requiring modifications to the datafile. Second order +corrections could be made to the background by repeating the solution +using the background fitted by the extraction procedure. + + +procedure extract_swath + +begin + extract raw swath from calibration frame + extract raw swath from background frame + return (calib swath minus bkg swath) +end + + +The algorithm used to simultaneously model all spectra in a swath from +across the dispersion is probably the most difficult and time consuming +part of the problem. The problem is nonlinear in all but one of the +four or more parameters for each spectra. You have spent a lot of +time on this and we are probably not going to be able to improve on +your algorithms significantly, though the generation of the matrix in +each step can probably be optimized significantly. + +The analytic line-profile model is the most general and should work +for most instruments with small circular apertures, even in the +presence of severe distortions. It should be possible, however, to +fit a simpler model given by a lookup table, solving only for the +position and height of each spectra. This model may be adequate for +many instruments, should be must faster to fit, and may produce more +accurate results since there are fewer parameters in the fit. The +disadvantage of an empirical model is that it must be accurately +interpolated (including the derivatives), requiring use of spline +interpolation or a similar technique (I have tried linear and it is +not good enough). Vesa has implemented procedures for fitting splines +and evaluating their derivatives. + +Fitting the empirical model simultaneously to any number of spectra +should be straightforward provided the signal to noise is reasonable, +since there are few parameters in the fit and the matrix is banded +(the Marquardt algorithm would work fine). However, if you ever have +to deal with data where a very faint or nonexistent spectra is next to +a bright one, it may be difficult to constrain the fit. I doubt if +the present approach of smoothing the coeff across the dispersion and +iterating would work in such a case. The best approach might be to +fix the center of the faint spectra relative to the bright one once +the signal drops below a certain level, or to drop it from the fit +entirely. This requires that the matrix be able to change size during + + + -4- + MULTISPEC (Dec83) Multispec Algorithms MULTISPEC (Dec83) + + + +the fit. + +algorithm fit_empirical_model + +begin + [upon entry, we already have an initial estimate of the coeff] + + # Marquardt (gradient expansion) algorithm. Make 2nd order + # Taylor's expansion to chisquare near minimum and solve for + # correction vector which puts us at minimum (subject to + # Taylor's approx). Taylor's approximation rapidly becomes + # better as we near the minimum of the multidimensional + # chisquare, hence convergence is extremely rapid given a good + # starting estimate. + + repeat { + evaluate curvature matrix using current coeff. + solve banded curvature matrix + + compute error matrix + for (each spectra) + if (uncertainty in center coeff > tol) { + fix center by estimation given relative spacing + in higher signal region + remove spectra center from solution + } + + if (no center coeff were rejected) + tweak correction vector to accelerate convergence + new coeff vector = old coeff vector + correction vector + compute norm of correction vector + } until (no more center coeff rejected and norm < tolerance) + + compute final uncertainties +end + + +The following is close to what is currently done to fit the analytic +model, as far as I can remember (I have modified it slightly to +stimulate discussion). The solution is broken up into two parts to +reduce the number of free parameters and increase stability. If the +uncertainty in a free parameter becomes large it is best to fix the +parameter (it is particularly easy for this data to estimate all but +the intensity parameter). A fixed parameter is used in the model and +affects the solution but is not solved for (i.e., like the background). + +The analytic fit will be rather slow, even if the outer loop is +constrained to one iteration. If it takes (very rough estimates) .5 +sec to set up the banded matrix and .3 sec to solve it, 3 iterations +to convergence, we have 5 sec per swath. If we have an 800 lines +broken into swaths of 32 lines, the total is 125 sec per image (to +within a factor of 5). + + + + -5- + MULTISPEC (Dec83) Multispec Algorithms MULTISPEC (Dec83) + + + +algorithm fit_analytic_model + +begin + [upon entry, we already have an initial estimate of the coeff] + + repeat { + save coeff + solve for center,height,width of each line with second + order terms fixed (but not necessarily zero) + apply constraints on line centers and widths + repeat solution adding second order coeff (shape terms) + + compute error matrix + for (each coeff) + if (uncertainty in coeff > tol) { + fix coeff value to reasonable estimate + remove coeff from solution + } + + compute total correction vector given saved coeff + if (no coeff were rejected) + tweak correction vector to accelerate convergence + compute norm of correction vector + } until (no additional coeff rejected and norm < tolerance) + + compute final uncertainties +end + + + +3.2 Extraction + + The purpose of extraction is to compute the integral of the spectra +across the dispersion, producing I(y) for each spectra. An estimate of +the uncertainty U(y) should also be produced. The basic extraction +techniques are summarized below. The number of spectra, spectra +centers, spectra width and shape parameters are taken from the model +fitted to the calibration frame as outlined above. We make a +simultaneous solution for the profile heights and the background (a +linear problem), repeating the solution independently for each line in +the image. For a faint spectrum, it is essential to determine the +background accurately, and we can do that safely here since the matrix +will be very well behaved. + + (1) Aperture sum. All of the pixels within a specified radius of + the spectra are summed to produce the raw integral. The + background image is also summed and subtracted to yield the + final integral. The radius may be a constant or a function of + the width of the profile. Fractional pixel techniques should + be used to minimize sampling effects at the boundaries of the + aperture. Pixel rejection is not possible since there is no + fitted surface. The model is used only to get the spectra + center and width. This technique is fastest and may be best + + + -6- + MULTISPEC (Dec83) Multispec Algorithms MULTISPEC (Dec83) + + + + if the profile is difficult to model, provided the spectra are + not crowded. + + (2) Weighted aperture sum. Like (1), except that a weighting + function is applied to correct for the effects of crowding. + The model is fitted to each object line, solving for I + (height) and B (background) with all other parameters fixed. + This is a linear solution of a banded matrix and should be + quite fast provided the model can be sampled efficiently to + produce the matrix. It is possible to iterate to reject bad + pixels. The weight for a spectra at a data pixel is the + fractional contribution of that spectra to the integral of the + contributions of all spectra. + + (3) Fit and integrate the model. The model is fitted as in (2) to + the data pixels but the final integral is produced by + integrating the model. This technique should be more + resistant to noise in the data than is (2), because we are + using the high signal to noise information in the model to + constrain the integral. More accurate photometric results + should therefore be possible. + + +Method (2) has the advantage that the integral is invariant with +respect to scale errors in the fitted models, provided the same error +is made in each model. Of course, the same error is unlikely to be +made in all models contributing to a point; it is more likely that an +error will put more energy into one spectra at the expense of its +neighbors. In the limit as the spectra become less crowded, however, +the effects of errors in neighboring spectra become small and the +weighted average technique looks good; it becomes quite insensitive to +errors in the model and in the fit. For crowded spectra there seems +no alternative to a good multiparameter fit. For faint spectra method +(3) is probably best, and fitting the background accurately becomes +crucial. + +In both (2) and (3), subtraction of the scaled models yields a residual +image which can be used to evaluate at a glance the quality of the fit. +Since most all of the effort in (2) and (3) is in the least squares +solution and the pixel rejection, it might be desirable to produce two +integrals (output spectra), one for each algorithm, as well as the +uncertainty vector (computed from the covariance matrix, not the +residual). + + + +3.3 Smoothing Coefficient Arrays + + In several places we have seen the need for smoothing coefficient +arrays. The use of polynomials for smoothing is questionable unless +the order of the polynomial is low (3 or less). High order +polynomials are notoriously bad near the endpoints of the fitted +array, unless the data curve happens to be a noisy low order + + + -7- + MULTISPEC (Dec83) Multispec Algorithms MULTISPEC (Dec83) + + + +polynomial (rare, to say the least). Convolution or piecewise +polynomial functions (i.e., the natural cubic smoothing spline) should +be considered if there is any reason to believe that the coefficient +array being smoothed may have high frequency components which are +physical and must be followed (i.e., a bend or kink). + + + +3.4 Weighting (pg. 11) + + The first weighting scheme (1 / sqrt (data)) seems inverted to me. +The noise goes up as with the signal, to be sure, but the signal to +noise usually goes up faster. Seems to me the weight estimate should +be sqrt(data). It also make more sense to weight the least blended +(peak) areas most. + + + +3.5 Rejection criteria (pg. 13) + + The same comments apply to this rejection criterium as in section +2. I assume that "(data - model)" is supposed to be "abs (data - +model"). + + + +3.6 Uncertainties and Convergence Criteria + + I got the impression that you were using the residual of the data +minus the fitted surface both as the convergence criterium and as a +measure of the errors in the fit. It is neither; assuming a perfect +model, the residual gives only a measure of the noise in the data. + +Using the residual to establish a convergence criterium seems +reasonable except that it is hard to reliably say what the criterium +should be. Assuming that the algorithm converges, the value of the +residual when convergence is acheived is in general hard to predict, +so it seems to me to be difficult to establish a convergence +criterium. The conventional way to establish when a nonlinear fit +converges is by measuring the norm of the correction vector. When the +norm becomes less than some small number the algorithm is said to have +converged. The multidimensional chisquare surface is parabolic near +the minimum and a good nonlinear algorithm will converge very rapidly +once it gets near the minimum. + +The residual is a measure of the overall goodness of fit, but tells us +nothing about the uncertainties in the individual coefficients of the +model. The uncertainties in the coefficients are given by the +covariance or error matrix (see Bevington pg. 242). It is ok to push +forward and produce an extraction if the algorithm fails to converge, +but ONLY provided the code gives a reliable estimate of the +uncertainty. + + + + -8- + MULTISPEC (Dec83) Multispec Algorithms MULTISPEC (Dec83) + + + +3.6 Evaluating the Curvature Matrix Efficiently + + The most expensive part of the reduction process is probably +evaluating the model to form the curvature matrix at each iteration in +the nonlinear solution. The most efficient way to do this is to use +lookup tables. If the profile shape does not change, the profile can +be sampled, fitted with a spline, and the spline evaluated to get the +zero through second derivatives for the curvature matrix. This can be +done even if the width of the profile changes by adding a linear +term. If the shape of the profile has to change, it is still possible +to sample either the gaussian or the exponential function with major +savings in computation time. + + + +3.7 Efficient Extraction (pg. 12) + + The reported time of 3 cpu hours to extract the spectra from an +800 line image is excessive for a linear solution. I would estimate +the time required for the 800 linear banded matrix solutions at 4-8 +minutes, with a comparable time required for matrix setup if it is +done efficiently. I suspect that the present code is not setting up +the linear banded matrix efficiently (not sampling the model +efficiently). Pixel rejection should not seriously affect the timings +assuming that bad pixels are not detected in most image lines. + + + +4. Correcting for Variations in the PSF + + For all low signal to noise data it is desirable to correct for +variations in the point spread function, caused by variable focus, +scattering, or whatever. This does not seem such a difficult problem +since the width of the line profile is directly correlated with the +width of the PSF and the information is provided by the current model +at each point in each extracted spectrum. The extracted spectra can +be corrected for the variation in the PSF by convolution with a spread +function the width of which varies along the spectrum. diff --git a/noao/twodspec/multispec/doc/MSalgo_c.hlp b/noao/twodspec/multispec/doc/MSalgo_c.hlp new file mode 100644 index 00000000..4b9c3356 --- /dev/null +++ b/noao/twodspec/multispec/doc/MSalgo_c.hlp @@ -0,0 +1,449 @@ +.help multispec Dec83 "Multispec Algorithms" +.ce +Algorithms for the Multi-Spectra Extraction Package +.ce +Analysis and Discussion +.ce +December 2, 1983 + +.sh +1. Disclaimer + + This should not be taken as a statement of how the algorithms of the +final package should function; this is merely an analysis and discussion +of the algorithms, and should be followed by further discussion before we +decide what course to follow in the final package. We may very well decide +that the level of effort required to implement rigorously correct nonlinear +fitting algorithms is not justified by the expected scientific usage of +the package. Before we can decide that, though, we need an accurate estimate +of the level of effort required. + +In attacking nonlinear surface fitting problems it is important to recognize +that almost any techniques can be made to yield a result without the program +crashing. Production of a result (extraction of a spectrum) does not mean +that the algorithm converged, that the solution is unique, that the model +is accurate, or that the uncertainties in the computed coefficients have +been minimized. + +.sh +2. Multispec Flat (pg. 4) + + This sounds like a classical high pass filter and might be best implemented +via convolution. Using a convolution operator with a numerical kernel has +the advantage that the filter can be easily modifed by resampling the kernel +or by changing the size of the kernel. It is also quite efficient. The +boundary extension feature of IMIO makes it easy to deal with the problem of +the kernel overlapping the edge of the image in a convolution. Since the +convolution is one-dimensional (the image is only filtered in Y), it will +always be desirable to transpose the image. + +The method used to detect and reject bad pixels (eqn 1) is not correct. +The rejection criteria should be invariant with respect to a scaling of the +pixel values. If the data has gone through very much processing (i.e., +dtoi on photographic data), the relation between photon counts and pixel value +may be linear, but the scale is unknown. Rejection by comparison of a data +value to a "smoothed" value is more commonly done as follows: + + reject if: abs (observed - smoothed) > (K * sigma) + +where sigma is the noise sigma of the data, generally a function of the signal. + +It is often desirable in rejection algorithms to be able to specify, +as an option, that all pixels within a specified radius of a bad pixel +be rejected, rather than just the pixel which was detected. This is only +unnecessary if the bad pixels are single pixel events (no wings). Region +rejection makes an iterative rejection scheme converge faster, as well as +rejecting the faint wings of the contaminated region. + +.sh +2.1 Dividing by the Flat (pg. 5) + + There is no mention of any need for registering the flat with the data +field. Is it safe to assume that the quartz and the object frames are +precisely registered? What if the user does in fact average several quartz +frames taken over a period of time? (Image registration is a general +problem that is probably best left until solved in IMAGES). + +.sh +3. Multiap Extraction (pg. 5-6, 8-13) + + The thing that bothers me most about the modeling and extraction +process is that the high signal to noize quartz information is not used to +full advantage, and the background is not fitted very accurately. The +present algorithms will work well for high signal to noise data, but +will result in large (percentage) errors for faint spectra. + +Basically, it seems to me that the high signal to noise quartz spectra +should, in many cases, be used to determine the position and shape of the +spectral lines. This is especially attractive since the quartz and spectra +appear to be closely registered. Furthermore, if the position-shape solution +and extraction procedures are separate procedures, there is nothing to prevent +one from applying both to the object spectum if necessary for some reason +(i.e., poor registration, better signal to noise in the object spectrum in +the region of interest, signal dependent distortions, lack of a quartz image, +etc., would all justify use of the object frame). It should be possible to +model either the quartz or the object frame, and to reuse a model for more +than one extraction. + +Let us divide the process up into two steps, "modeling", and "extraction" +(as it is now). The "calibration frame" may be the quartz, an averaged +quartz, or the object frame. Ideally it will have a high signal to noise +ratio and any errors in the background should be negligible compared to +the signal. + +We do not solve for the background while modeling the calibration frame; +we assume that the background has been fitted by any of a variety of +techniques and a background frame written before the calibration frame +is modeled. A "swath" is the average of several image lines, where an +image line runs across the dispersion, and a column along the dispersion. + +.sh +3.1 Modeling + + I would set the thing up to start fitting at any arbitrary swath, rather +than the first swath, because it not much harder, and there is no guarantee +that the calibration frame will have adequate signal to noise in the first +swath (indeed often the lowest signal to noise will be found there). +We define the "center" swath as the first swath to be fitted, corresponding +to the highest signal to noise region of the calibration frame. By default +the center swath should be the swath used by find_spectra, especially if +there is significant curvature in the spectra. + +.ks +.nf +algorithm model_calibration_frame + +begin + extract center swath + initialize coeff using centers from find_spectra + model center swath (nonlinear) + + for (successive swaths upward to top of frame) { + extract swath + initialize coeff to values from last fit + model swath (nonlinear) + save coeff in datafile + } + + set last-fit coeff to values for center swath + for (successive swaths downward to bottom of frame) { + extract swath + initialize coeff to values from last fit + model swath (nonlinear) + save coeff in datafile + } + + smooth model coeff (excluding intensity) along the dispersion + [high freq variations in spectra center and shape from line] + [to line are nonphysical] + variance of a coeff at line-Y from the smoothed model value is + a measure of the uncertainty in that coeff. +end +.fi +.ke + + +I would have the background fitting routine write as output a background +frame, the name of which would be saved in the datafile, rather than saving +the coeff of the bkg fit in the datafile. The background frame may then +be produced by any of a number of techniques; storing the coefficients of +the bkg fit in the datafile limits the technique used to a particular model. +For similar reasons, the standard bkg fitting routine should be broken up +into a module which determines the region to be fitted, and a module which +fits the bkg pixels and writes the bkg image. + +For example, if the default background fitting routine is a line by line +routine, the output frame could be smoothed to remove the (nonphysical) +fluctuations in the background from line to line. A true two dimensional +background fitting routine may be added later without requiring modifications +to the datafile. Second order corrections could be made to the background +by repeating the solution using the background fitted by the extraction +procedure. + + +.ks +.nf +procedure extract_swath + +begin + extract raw swath from calibration frame + extract raw swath from background frame + return (calib swath minus bkg swath) +end +.fi +.ke + + +The algorithm used to simultaneously model all spectra in a swath from +across the dispersion is probably the most difficult and time consuming +part of the problem. The problem is nonlinear in all but one of the four +or more parameters for each spectra. You have spent a lot of time on this +and we are probably not going to be able to improve on your algorithms +significantly, though the generation of the matrix in each step can +probably be optimized significantly. + +The analytic line-profile model is the most general and should work for most +instruments with small circular apertures, even in the presence of severe +distortions. It should be possible, however, to fit a simpler model given +by a lookup table, solving only for the position and height of each spectra. +This model may be adequate for many instruments, should be must faster to +fit, and may produce more accurate results since there are fewer parameters +in the fit. The disadvantage of an empirical model is that it must be +accurately interpolated (including the derivatives), requiring use of spline +interpolation or a similar technique (I have tried linear and it is not good +enough). Vesa has implemented procedures for fitting splines and evaluating +their derivatives. + +Fitting the empirical model simultaneously to any number of spectra should +be straightforward provided the signal to noise is reasonable, since there +are few parameters in the fit and the matrix is banded (the Marquardt +algorithm would work fine). However, if you ever have to deal with data +where a very faint or nonexistent spectra is next to a bright one, it may +be difficult to constrain the fit. I doubt if the present approach of +smoothing the coeff across the dispersion and iterating would work in such +a case. The best approach might be to fix the center of the faint spectra +relative to the bright one once the signal drops below a certain level, +or to drop it from the fit entirely. This requires that the matrix be able +to change size during the fit. + +.ks +.nf +algorithm fit_empirical_model + +begin + [upon entry, we already have an initial estimate of the coeff] + + # Marquardt (gradient expansion) algorithm. Make 2nd order + # Taylor's expansion to chisquare near minimum and solve for + # correction vector which puts us at minimum (subject to + # Taylor's approx). Taylor's approximation rapidly becomes + # better as we near the minimum of the multidimensional + # chisquare, hence convergence is extremely rapid given a good + # starting estimate. + + repeat { + evaluate curvature matrix using current coeff. + solve banded curvature matrix + + compute error matrix + for (each spectra) + if (uncertainty in center coeff > tol) { + fix center by estimation given relative spacing + in higher signal region + remove spectra center from solution + } + + if (no center coeff were rejected) + tweak correction vector to accelerate convergence + new coeff vector = old coeff vector + correction vector + compute norm of correction vector + } until (no more center coeff rejected and norm < tolerance) + + compute final uncertainties +end +.fi +.ke + + +The following is close to what is currently done to fit the analytic +model, as far as I can remember (I have modified it slightly to stimulate +discussion). The solution is broken up into two parts to reduce the number +of free parameters and increase stability. If the uncertainty in a free +parameter becomes large it is best to fix the parameter (it is particularly +easy for this data to estimate all but the intensity parameter). A fixed +parameter is used in the model and affects the solution but is not solved +for (i.e., like the background). + +The analytic fit will be rather slow, even if the outer loop is constrained +to one iteration. If it takes (very rough estimates) .5 sec to set up the +banded matrix and .3 sec to solve it, 3 iterations to convergence, we have +5 sec per swath. If we have an 800 lines broken into swaths of 32 lines, +the total is 125 sec per image (to within a factor of 5). + + +.ks +.nf +algorithm fit_analytic_model + +begin + [upon entry, we already have an initial estimate of the coeff] + + repeat { + save coeff + solve for center,height,width of each line with second + order terms fixed (but not necessarily zero) + apply constraints on line centers and widths + repeat solution adding second order coeff (shape terms) + + compute error matrix + for (each coeff) + if (uncertainty in coeff > tol) { + fix coeff value to reasonable estimate + remove coeff from solution + } + + compute total correction vector given saved coeff + if (no coeff were rejected) + tweak correction vector to accelerate convergence + compute norm of correction vector + } until (no additional coeff rejected and norm < tolerance) + + compute final uncertainties +end +.fi +.ke + +.sh +3.2 Extraction + + The purpose of extraction is to compute the integral of the spectra +across the dispersion, producing I(y) for each spectra. An estimate of +the uncertainty U(y) should also be produced. The basic extraction techniques +are summarized below. The number of spectra, spectra centers, spectra width +and shape parameters are taken from the model fitted to the calibration +frame as outlined above. We make a simultaneous solution for the profile +heights and the background (a linear problem), repeating the solution +independently for each line in the image. For a faint spectrum, it is +essential to determine the background accurately, and we can do that safely +here since the matrix will be very well behaved. +.ls 4 +.ls (1) +Aperture sum. All of the pixels within a specified radius of the spectra +are summed to produce the raw integral. The background image is also summed +and subtracted to yield the final integral. The radius may be a constant or a +function of the width of the profile. Fractional pixel techniques should +be used to minimize sampling effects at the boundaries of the aperture. +Pixel rejection is not possible since there is no fitted surface. The model +is used only to get the spectra center and width. This technique is fastest +and may be best if the profile is difficult to model, provided the spectra +are not crowded. +.le +.ls (2) +Weighted aperture sum. Like (1), except that a weighting function is +applied to correct for the effects of crowding. The model is fitted +to each object line, solving for I (height) and B (background) with all +other parameters fixed. This is a linear solution of a banded matrix and +should be quite fast provided the model can be sampled efficiently to +produce the matrix. It is possible to iterate to reject bad pixels. +The weight for a spectra at a data pixel is the fractional contribution +of that spectra to the integral of the contributions of all spectra. +.le +.ls (3) +Fit and integrate the model. The model is fitted as in (2) to the data +pixels but the final integral is produced by integrating the model. +This technique should be more resistant to noise in the data than is (2), +because we are using the high signal to noise information in the model to +constrain the integral. More accurate photometric results should therefore +be possible. +.le +.le + + +Method (2) has the advantage that the integral is invariant with respect +to scale errors in the fitted models, provided the same error is made in +each model. Of course, the same error is unlikely to be made in all +models contributing to a point; it is more likely that an error will put +more energy into one spectra at the expense of its neighbors. In the limit +as the spectra become less crowded, however, the effects of errors in +neighboring spectra become small and the weighted average technique looks +good; it becomes quite insensitive to errors in the model and in the fit. +For crowded spectra there seems no alternative to a good multiparameter +fit. For faint spectra method (3) is probably best, and fitting the +background accurately becomes crucial. + +In both (2) and (3), subtraction of the scaled models yields a residual +image which can be used to evaluate at a glance the quality of the fit. +Since most all of the effort in (2) and (3) is in the least squares solution +and the pixel rejection, it might be desirable to produce two integrals +(output spectra), one for each algorithm, as well as the uncertainty vector +(computed from the covariance matrix, not the residual). + +.sh +3.3 Smoothing Coefficient Arrays + + In several places we have seen the need for smoothing coefficient arrays. +The use of polynomials for smoothing is questionable unless the order of +the polynomial is low (3 or less). High order polynomials are notoriously +bad near the endpoints of the fitted array, unless the data curve happens +to be a noisy low order polynomial (rare, to say the least). Convolution or +piecewise polynomial functions (i.e., the natural cubic smoothing spline) +should be considered if there is any reason to believe that the coefficient +array being smoothed may have high frequency components which are physical and +must be followed (i.e., a bend or kink). + +.sh +3.4 Weighting (pg. 11) + + The first weighting scheme (1 / sqrt (data)) seems inverted to me. +The noise goes up as with the signal, to be sure, but the signal to noise +usually goes up faster. Seems to me the weight estimate should be sqrt(data). +It also make more sense to weight the least blended (peak) areas most. + +.sh +3.5 Rejection criteria (pg. 13) + + The same comments apply to this rejection criterium as in section 2. +I assume that "(data - model)" is supposed to be "abs (data - model"). + +.sh +3.6 Uncertainties and Convergence Criteria + + I got the impression that you were using the residual of the data minus +the fitted surface both as the convergence criterium and as a measure of the +errors in the fit. It is neither; assuming a perfect model, the residual gives +only a measure of the noise in the data. + +Using the residual to establish a convergence criterium seems reasonable +except that it is hard to reliably say what the criterium should be. +Assuming that the algorithm converges, the value of the residual when +convergence is achieved is in general hard to predict, so it seems to me to +be difficult to establish a convergence criterium. The conventional way +to establish when a nonlinear fit converges is by measuring the norm of +the correction vector. When the norm becomes less than some small number +the algorithm is said to have converged. The multidimensional chisquare +surface is parabolic near the minimum and a good nonlinear algorithm will +converge very rapidly once it gets near the minimum. + +The residual is a measure of the overall goodness of fit, but tells us +nothing about the uncertainties in the individual coefficients of the model. +The uncertainties in the coefficients are given by the covariance or error +matrix (see Bevington pg. 242). It is ok to push forward and produce an +extraction if the algorithm fails to converge, but ONLY provided the code +gives a reliable estimate of the uncertainty. + +.sh +3.6 Evaluating the Curvature Matrix Efficiently + + The most expensive part of the reduction process is probably evaluating +the model to form the curvature matrix at each iteration in the nonlinear +solution. The most efficient way to do this is to use lookup tables. +If the profile shape does not change, the profile can be sampled, fitted +with a spline, and the spline evaluated to get the zero through second +derivatives for the curvature matrix. This can be done even if the width +of the profile changes by adding a linear term. If the shape of the profile +has to change, it is still possible to sample either the gaussian or the +exponential function with major savings in computation time. + +.sh +3.7 Efficient Extraction (pg. 12) + + The reported time of 3 cpu hours to extract the spectra from an 800 line +image is excessive for a linear solution. I would estimate the time required +for the 800 linear banded matrix solutions at 4-8 minutes, with a comparable +time required for matrix setup if it is done efficiently. I suspect that the +present code is not setting up the linear banded matrix efficiently (not +sampling the model efficiently). Pixel rejection should not seriously affect +the timings assuming that bad pixels are not detected in most image lines. + +.sh +4. Correcting for Variations in the PSF + + For all low signal to noise data it is desirable to correct for variations +in the point spread function, caused by variable focus, scattering, or +whatever. This does not seem such a difficult problem since the width of +the line profile is directly correlated with the width of the PSF and the +information is provided by the current model at each point in each extracted +spectrum. The extracted spectra can be corrected for the variation in the +PSF by convolution with a spread function the width of which varies along +the spectrum. +.endhelp diff --git a/noao/twodspec/multispec/doc/MSspecs.doc b/noao/twodspec/multispec/doc/MSspecs.doc new file mode 100644 index 00000000..09955e9c --- /dev/null +++ b/noao/twodspec/multispec/doc/MSspecs.doc @@ -0,0 +1,698 @@ + MULTISPEC (Oct83) Multi-Spectra Extraction Package MULTISPEC (Oct83) + + + + Detailed Specifications for the Multi-Spectra Extraction Package + F. Valdes + December 8, 1983 + + +1. Introduction + + The multi-spectra extraction package (MULTISPEC) provides the +basic tools for modeling, cleaning, and extracting spectra from images +containing multiple aperture spectra running roughly parallel. These +tools will generally be combined in reduction script tasks but may +also be used directly for non-standard analysis. + + This design presents the requirements and specifications for the +MULTISPEC package. Details concerning the algorithms are given in a +separate document, Algorithms for the Multi-Spectra Extraction Package. + + +2. Input Data Requirements + + The input data for the MULTISPEC package consists of image files +containing one or more aperture spectra. The spectra are required to +run roughly parallel to each other and parallel to the second +digitization axis. The latter requirement may require a general +rotation and interpolation image operator. The images are assumed to +be corrected to linear relative intensity. Thus, the steps of +correcting digital detector images for dark current, bias, and +pixel-to-pixel sensitivity variations must be performed before using +the MULTISPEC tasks. + + Because the the MULTISPEC package is being developed concurrently +with the IRAF standard image processing tools this document specifies +the requirements for the preliminary image processing needed to +prepare digital detector images for the MULTISPEC package. + + +2.1 Basic Digital Detector Reduction Tasks + + The prelimary reduction of multi-spectra images uses CL scripts +based on general image operators. Some of the scripts are for +specific instruments or specific reduction applications and some are +generally useful image processing tasks. The scripts allow the +specification of many images for which the operations will be +repetitively applied. + + The following CL scripts are required to reduce multi-spectra +images from digital detectors. + + + debias multispec_flat flat_divide + + + + + + -1- + MULTISPEC (Oct83) Multi-Spectra Extraction Package MULTISPEC (Oct83) + + + +debias + The files in a list of filenames are automatically debiased and + trimmed. This routine will be instrument specific but used by + other reduction tasks beyond MULTISPEC. + +multispec_flat + The files in a list of quartz multi-spectra filenames are added, + the result is smoothed along the dispersion dimension, and then + the original image is divided by the smoothed image to produce a + flat field image. The unsmoothed to smoothed ratio is computed + only if the value of the smoothed pixel is greater than a + specified amount. Otherwise, the ratio is set to unity. This + routine is not instrument specific but is used only for MULTISPEC + reductions. + +flat_divide + The files in a list of filenames are each divided by a specified + flat field image. This routine is not instrument or application + specific. + + The required general image processing programs needed to implement +these scripts are specified below. + + +(1) A routine to compute the average value from a specified area of the + image. Used to determine the average bias value from a bias strip. + +(2) A routine to trim a specified portion of an image. Used to trim + the bias strip. + +(3) Routines to multiply and subtract images by a constant. Used to + scale images such as dark exposures and to remove the average bias + value as obtained by (1) above. + +(4) Routines to subtract, add, and divide images. Used to subtract + dark current and bias exposures, to add several exposures to + increase the signal-to-noise, and to divide by a flat field image. + The divide routine must give the user the option to substitute a + constant or ignore any divisions in which the denominator is less + than a specified value. + +(5) A routine to rotate or transpose an image. Used to align the + spectra along lines or columns. + +(6) A routine to apply a filter to lines of the image. For + multi-spectra images a smooth quartz is produced by using a + running quadratic filter along each line of the dispersion + dimension. The filter must be able to recognize bad pixels + (specified by a user defined threshold) and remove them from the + filtering operation. + + + + + + -2- + MULTISPEC (Oct83) Multi-Spectra Extraction Package MULTISPEC (Oct83) + + + +3. Requirements for the MULTSPEC Package + + The MULTISPEC package shall satisfy the following requirements. + +(1) The component programs shall be CL callable. + +(2) The programs shall interact only through image files and MULTISPEC + data files. + +(3) It shall be possible to extract spectra without modeling. + +(4) The entire image shall be extracted and not limited by failures in + the algorithms. + +(5) It shall be possible to specify specific lines or swaths in the + image on which to operate. + +(6) CL scripts shall be provided for the common data sources. These + scripts will work automatically. + +The follow functions shall be provided: + +o Make an initial rough but automated identification of the spectra + locations. + +o Provide for a user identification list for the spectra locations. + This list shall be of the standard image cursor type to allow + generation of the list with the standard image cursor programs. + +o Determine and correct for a slowly varying background. + +o Reliably and accurately trace spectra in the presence of geometric + distortions (pincushion, s, shear, etc.). + +o Extract spectra by one of: + + a. Strips of constant width about the located spectra. The width + may be specified to fractions of a pixel and the extraction + will use fractional pixel interpolation. l + + b. Strips of width proportional to a Gaussian width parameter. + + c. Modeling to obtain estimates of the total luminosity. The + estimate will be the integral of the model. + + d. Summation of the data pixel values with fractional + contributions of the pixel value to the spectra based on + modeling. + +o An option shall be available to specify whether to ignore blank + pixels or use interpolated values. + + o Programs shall be provided to produce data files which can be + + + -3- + MULTISPEC (Oct83) Multi-Spectra Extraction Package MULTISPEC (Oct83) + + + + accessed by one dimensional spectroscopic reduction routines. + At a minimum these formats shall include: + + a. Reduction to an image file consisting of one line per + extracted spectrum + + b. The standard IIDS format available with the CYBER + Multi-Aperture Plate programs + + +3.2 Modeling Requirements + + The modeling of multi-spectra images, particularly in the case of +blended spectra, shall: + +(1) Model blended spectra with sufficient reliability and robustness + that a reasonable solution is always obtained, though of possibly + limited usefulness. + +(2) The modeling shall provide estimates for the uncertainties in the + fitted parameters as a function of position along the spectrum. + +(3) Remove cosmic rays and other defective pixels by reference to the + model. + +(4) Allow the transfer of a model solution for one image to another + image. + +(5) Display numerically and graphically the data, the fitted model, and + the residuals. + + +4. Program Specifications + + +4.1 Basic Programs + + The basic programs of the package are general purpose tools which +initialize a MULTISPEC data file and perform a single fundamental +operation on the data in the MULTISPEC data file. There is one data +file associated with each image. The data file is hidden from the +user and so the user need not be aware of the data file. The data +files are referenced only the image filename specified in the program +parameters. The data files contain such information as a processing +history, the spectra positions and extracted luminosities, the model +parameters (one set for each spectra for each modelled image line (or +swath), etc. The programs generally are allowed to specify specific +lines, columns, and/or spectra on which to operate. The line, column +and spectra specifications are given as strings which contain numbers +separated by whitespace, commas, and the range indicator "-". The +script tasks of section 4.2 will combine these basic programs to +perform a general multi-spectra extraction. + + + + -4- + MULTISPEC (Oct83) Multi-Spectra Extraction Package MULTISPEC (Oct83) + + + + ap_plate copy_params find_spectra convolve + fit_bckgrnd find_bckgrnd line_list model_extrac + model_fit model_image model_list sigma_extract + strip_extract to_iids to_image to_onedspec + +ap_plate + The information from an on-line data file containing descriptions + of all the aperture plates prepared at Kitt Peak is read to find a + specified aperture plate. The drilled aperture positions are + correlated with the spectra in the image to deduce relative + wavelength offsets. The identifications for the spectra as well + as other auxiliary information is recorded in the data file. If + no image file is specified then only the aperture plate + information is printed. This program is used in the + MULTIAP_EXTRACT program. This program is not essential to the + operation of the MULTISPEC package. + + Multi-Spectra image image = + Aperture plate plate = + (mode = ql) + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + -5- + MULTISPEC (Oct83) Multi-Spectra Extraction Package MULTISPEC (Oct83) + + + +The Background + The are two possibilities for dealing with the background. In the + first case, FIT_BCKGRND, the background will be fitted by + polynomials and the coefficients stored in the MULTISPEC data + file. These coefficients are then used by the other programs to + estimate the background at the spectra. The second option, + FIND_BCKGRND, generates a background image in which the spectra + and other selected areas are set to blank pixels. Then a general + image interpolator is used fill in the blank pixels with background + estimates. The other MULTISPEC programs will then access this + background frame. The background frame image name will be stored + in the MULTISPEC data file and the image header. + + + fit_bckgrnd + Fit a background in a MULTISPEC image by a polynomial using + pixels not near the spectra and in the user specified swaths + and columns. The buffer distance is in pixels and refers to a + minimum distance from the center of any spectrum beyond which + the background pixels are found. Blank pixels are ignored in + the background fit. Deviant pixels will be rejected. + + Multi-Spectra image image = + Buffer from spectra buffer = 12 + Polynomial order order = 3 + Lines per swath (lines_per_swath = 32) + Swaths to fit (swaths = 1-1000) + Columns to fit (columns = 1-1000) + Rejection threshold (threshold = 5) + Print general diagnostics (verbose = no) + (mode = ql) + + find_bckgrnd + The spectra within a buffer distance and specified areas are + set to blank pixels and the remaining pixels copied to a + background image file. + + Multi-Spectra image image = + Background image background = + Buffer from spectra buffer = 12 + Lines to ignore (lines = ) + Columns to ignore (columns = ) + (mode = ql) + +convolve + A program will be provided to reduce either the extracted spectrum + or the modeled image to a common point-spread function. + + + + + + + + + -6- + MULTISPEC (Oct83) Multi-Spectra Extraction Package MULTISPEC (Oct83) + + + +copy_params + Create a MULTISPEC data file for a new image using appropriate + MULTISPEC parameters from an old image. The old image must have + been processed to find the spectra using FIND_SPECTRA and possibly + model fit. + + Old Multi-Spectra image old_image = + New Multi-Spectra image new_image = + (mode = ql) + +find_spectra + Initially locate the spectra in a MULTISPEC image. The positions + of the spectra within the range of columns are determined for the + starting line and then the spectra are tracked within the range of + lines. The minimum separation and minimum width would generally + be set for a particular instrument. If the automatic search is + not used then a list of cursor positions is read from the standard + input. + + Multi-Spectra image image = + Automatic search auto = yes + Starting line start_line = + Minimum separation (min_sep = 1) + Minimum width (min_width = 1) + Averaging width (average = 32) + Lines to search (lines = 1-1000) + Columns to search (columns = 1-1000) + Print general diagnostics (verbose = no) + (mode = ql) + +line_list + For the specified lines in the image print the image column + number, data value (possibly as a swath average), the model value + at that point (i.e. the sum of the model contributions from all + the spectra), the background value, and the residual. Plotting + scripts may be written using this routine to show the quality of a + model fit. + + Multi-Spectra image image = + Lines to list (lines = 1-1000) + (mode = ql) + + + + + + + + + + + + + + + -7- + MULTISPEC (Oct83) Multi-Spectra Extraction Package MULTISPEC (Oct83) + + + +model_extract + A previously fitted model is used to extract the spectra total + luminosity by apportioning the data values to spectra in the ratio + indicated by the model. If the clean option is specified then the + model is used to detect pixels which deviate from the model by a + specified amount. The model value replaces the deviant pixel in + the extraction and, if specified, also in the image file. + + Multi-Spectra image image = + Lines to extract (lines = 1-1000) + Clean spectra (clean = yes) + Cleaning threshold (threshold = 5) + Modify image (modify = yes) + Print general diagnostics (verbose = no) + (mode = ql) + +model_fit + A specified model is iteratively fitted to the data in each of the + specified lines (or swaths) until the RMS residual fails to + decrease. The models are selected by a string. The possible + values are + + (null string) - initialize the model + i - fit only the intensity scale + ip - fit the intensity scale and the position + ips1 - fit the intensity scale, position, and one parameter shape + ips2 - fit the intensity scale, position, and two parameter shape + ips3 - fit the intensity scale, position, and three parameter shape + ips4 - fit the intensity scale, position, and four parameter shape + These models will be combined in a script to search for the best + fit. + + The initial shape parameters will generally be set by scripts for a + particular data reduction. + + Multi-Spectra image image = + Model type model = + Lines per swath (lines_per_swath = 32) + Swaths to model (swaths = 1-1000) + Initial shape1 (shape1 = .1 ) + Initial shape2 (shape2 = 0 ) + Initial shape3 (shape3 = 0 ) + Initial shape4 (shape4 = 5 ) + Print general diagnostics (verbose = no) + (mode = ql) + + + + + + + + + + + -8- + MULTISPEC (Oct83) Multi-Spectra Extraction Package MULTISPEC (Oct83) + + + +model_image + An image file of the fitted model is created. This image may then + be displayed or a residual image may be calculated and displayed. + + Multi-Spectra image image = + Model image model = + (mode = ql) + .nf + .le + .ls model_list + For the specified lines and spectra the model is listed. + The listing gives, for each spectra, + the spectrum number, the line number, the fitted position, + the estimated wavelength, the + extracted luminosity, the intensity scale, model width parameters, and + the background polynomial coefficients. This routine can be used in scripts + to plot the extracted spectra, the trend of width with wavelength, and so + forth. + + .nf + Multi-Spectra image image = + Lines to list (lines = 1-1000) + Spectra to list (spectra = 1-1000) + (mode = ql) + +sigma_extract + A previously fitted model is used to extract the spectra luminosity + within a specified sigma of the peak. Because the model is not + necessarily a Gaussian the sigma is used to compute the intensity + ratio of the cutoff to the peak assumining a Gaussian profile and + then the data is extracted to the point the model intensity falls + below that cutoff. If the clean option is specified then the + model is used to detect pixels which deviate from the model by a + specified amount. The model value replaces the deviant pixel in + the extraction and, if specified, also in the image file. + + Multi-Spectra image image = + Sigma extraction width width = 1. + Lines to extract (lines = 1-1000) + Clean spectra (clean = yes) + Cleaning threshold (threshold = 5) + Modify image (modify = yes) + Print general diagnostics (verbose = no) + (mode = ql) + + + + + + + + + + + + -9- + MULTISPEC (Oct83) Multi-Spectra Extraction Package MULTISPEC (Oct83) + + + +strip_extract + A strip of constant width about the spectra positions is extracted. + If cleanning is desired a smoothed estimate of the profile is + obtained by averaging a number of lines about the line to be + cleaned. After fitting for the intensity scale pixels are found + which deviate from the profile by a specified amount. The profile + value replaces the deviant pixel in the extraction and, if + specified, also in the image file. No prior modeling is required + to use this extraction routine. + + Multi-Spectra image image = + Strip extraction width width = 1. + Lines to extract (lines = 1-1000) + Clean spectra (clean = yes) + Cleaning threshold (threshold = 5) + Lines per profile average (averge_lines = 32) + Modify image (modify = yes) + Print general diagnostics (verbose = no) + (mode = ql) + +to_iids + For a specified prefix, files of the form prefix.nn, where nn is a + specified spectra number, are created containing the extracted + spectra for all the specified image files. The format of the + files is the IIDS format developed for the CYBER Multi-Aperture + Plate Extractions. + + Multi-Spectra image images = + IIDS filename prefix iids_file = + Spectra to format (spectra = 1-1000) + (mode = ql) + +to_image + An image file containing one line of the extracted luminosities + for each specified spectra in the specified MULTISPEC image. + + Multi-Spectra image in_image = + Extracted spectra image out_image = + Spectra (spectra = 1-1000) + (mode = ql) + +to_onedspec + The extractions are converted to an as yet to be specified format + for use in the ONEDSPEC reduction package. + + Multi-Spectra images images = + ONEDSPEC data file onedspec_file = + Spectra (spectra = 1-1000) + (mode = ql) + + + + + + + -10- + MULTISPEC (Oct83) Multi-Spectra Extraction Package MULTISPEC (Oct83) + + + +4.2 General MULTISPEC CL Scripts + + The general MULTISPEC CL scripts perform a series of steps needed +to extract the spectra from a specified list of image files. These +steps have been found to generally perform the desired extraction task +fully. + + + multiap_extract echelle_extract + +multiap_extract + The specified multi-aperture plate images are extracted. If no + starting solution image, one which has previously been extracted, + is specified then the script performs an automatic search for the + specified number of spectra. Otherwise the solution from the + starting image is used as the initial model. The background is + then determined. This is followed by a series of fitting steps on + swaths of data. (For further details on the fitting steps see the + Algorithms paper). A MODEL_EXTRACT and cleaning follows. + Finally, the extraction is correlated with the specified aperture + plate using AP_PLATE. If there was no starting image then this + extraction becomes the initial solution image. Subsequent images + are extracted starting from the initial solution image. + + Multi-Aperture images images = + Initial solution image initial = + Aperture plate number plate = + Number of spectra nspectra = + (mode = ql) + +echelle_extract + The specified echelle images are extracted. If no starting + solution image, one which has previously been extracted, is + specified then the script performs an automatic search for the + specified number of orders. Otherwise the solution from the + starting image is used as the initial starting point. The + background is then determined. Finally a STRIP_EXTRACT and + cleaning is performed. If there was no starting image then this + extraction becomes the initial solution image. Subsequent images + are extracted starting from the initial solution image. + + Echelle images images = + Initial solution image initial = + Number of orders norders = + Extraction width width = + (mode = ql) + + +5. Outline of a MULTISPEC Reduction + + The following outline is for the reduction of a cryogenic camera +multi-aperture plate. All the programmer supplied default values are +used. + + + -11- + MULTISPEC (Oct83) Multi-Spectra Extraction Package MULTISPEC (Oct83) + + + + (1) rcamera mtb, "ap165.", "s", "3-9" + (2) debias "ap165.*" + (3) multispec_flat "ap165.[36]", "ap165.flat" + (4) flat_divide "ap165.*", "ap165.flat" + (5) multiap_extract "ap165.*", "", 165, 50 + (6) to_onedspec "ap165.*", oned165 + + +(1) The data is read from the observing tape(s) using RCAMERA. The + image files created are ap165.3, ap165.4, ..., ap165.9. This is + easily accomplished by using the filename prefix "ap165." in the + RCAMERA program. The raw images may be examined at this point on + a display. + +(2) The images are debiased using DEBIAS with all the "ap165." files + specified. The debias program knows about the location of the + bias strip for the cryogenic camera. + +(3) A a flat field is created using MULTISPEC_FLAT in which the + desired quartz frames are specified and a flat field image + filename is defined. The created flat field image may be examined + on an image display if desired. + +(4) All the debiased images are divided by the flat field using + FLAT_DIVIDE. + +(5) The script MULTIAP_EXTRACT is run in which the aperture plate + number, the number of spectra, and the image files to be extracted + are specified. The number of spectra is found by examining an + image on an image display or by plotting a cut across the spectra + using a general image profile program. + +(6) Finally, the extracted spectra are formatted for the ONEDSPEC + package using TO_ONEDSPEC with the extracted images specified. diff --git a/noao/twodspec/multispec/doc/MSspecs.hlp b/noao/twodspec/multispec/doc/MSspecs.hlp new file mode 100644 index 00000000..92f285ed --- /dev/null +++ b/noao/twodspec/multispec/doc/MSspecs.hlp @@ -0,0 +1,659 @@ +.help multispec Oct83 "Multi-Spectra Extraction Package" +.sp 3 +.ce +Detailed Specifications for the Multi-Spectra Extraction Package +.ce +F. Valdes +.ce +December 8, 1983 +.sh +1. Introduction + + The multi-spectra extraction package (MULTISPEC) provides the basic tools +for modeling, cleaning, and extracting spectra from images +containing multiple aperture spectra running roughly parallel. +These tools will generally be combined in reduction script tasks +but may also be used directly for non-standard analysis. + + This design presents the requirements and specifications +for the MULTISPEC package. Details concerning the +algorithms are given in a separate document, Algorithms for the +Multi-Spectra Extraction Package. +.sh +2. Input Data Requirements + + The input data for the MULTISPEC package consists of image +files containing one or more aperture spectra. The spectra are +required to run roughly parallel to each other and parallel to the +second digitization axis. The latter requirement may require a +general rotation and interpolation image operator. The images are +assumed to be corrected to linear relative intensity. Thus, the +steps of correcting digital detector images for dark current, bias, and +pixel-to-pixel sensitivity variations must be performed before using +the MULTISPEC tasks. + + Because the MULTISPEC package is being developed +concurrently with the IRAF standard image processing +tools this document specifies the requirements for the preliminary +image processing needed to prepare digital detector images for the MULTISPEC +package. +.sh +2.1 Basic Digital Detector Reduction Tasks + + The prelimary reduction of multi-spectra images uses CL scripts +based on general image operators. +Some of the scripts are for specific instruments or specific +reduction applications and some are generally useful image processing +tasks. The scripts allow the specification of many images for which +the operations will be repetitively applied. + + The following CL scripts are required to reduce multi-spectra images +from digital detectors. + + +.nf + debias multispec_flat flat_divide +.fi +.ke +.ks +.ls 4 debias +The files in a list of filenames are automatically debiased and trimmed. +This routine will be instrument specific but used by other reduction +tasks beyond MULTISPEC. +.le +.ke +.ks +.ls multispec_flat +The files in a list of quartz multi-spectra filenames are added, +the result is smoothed +along the dispersion dimension, and then the original image is divided +by the smoothed image to produce a flat field image. The unsmoothed +to smoothed ratio is computed only if the value of the smoothed +pixel is greater than a specified amount. Otherwise, the ratio is set +to unity. This routine is not instrument specific but is used only +for MULTISPEC reductions. +.le +.ke +.ks +.ls flat_divide +The files in a list of filenames are each divided by a specified flat +field image. This routine is not instrument or application specific. +.le +.ke + + The required general image processing programs needed to implement these +scripts are specified below. + +.ls (1) +A routine to compute the average value from a specified area of the +image. Used to determine the average bias value from a bias strip. +.le +.ls (2) +A routine to trim a specified portion of an image. Used to trim the +bias strip. +.le +.ls (3) +Routines to multiply and subtract images by a constant. Used to scale +images such as dark exposures and to remove the average bias value as +obtained by (1) above. +.le +.ls (4) +Routines to subtract, add, and divide images. Used to subtract dark +current and bias exposures, to add several exposures to increase the +signal-to-noise, and to divide by a flat field image. +The divide routine must give the user the option to substitute a constant or +ignore any divisions in which the denominator is less than a specified value. +.le +.ls (5) +A routine to rotate or transpose an image. Used to align the spectra +along lines or columns. +.le +.ls (6) +A routine to apply a filter to lines of the image. For multi-spectra images +a smooth quartz is produced by using a running quadratic filter along each +line of the dispersion dimension. The filter must be able to recognize +bad pixels (specified by a user defined threshold) and remove them from the +filtering operation. +.le +.sh +3. Requirements for the MULTSPEC Package + + The MULTISPEC package shall satisfy the following requirements. +.ls (1) +The component programs shall be CL callable. +.le +.ls (2) +The programs shall interact only through image files and MULTISPEC data files. +.le +.ls (3) +It shall be possible to extract spectra without modeling. +.le +.ls (4) +The entire image shall be extracted and not limited by failures in the +algorithms. +.le +.ls (5) +It shall be possible to specify specific lines or swaths in the image +on which to operate. +.le +.ls (6) +CL scripts shall be provided for the common data sources. These scripts +will work automatically. +.le + +The follow functions shall be provided: +.ls o +Make an initial rough but automated identification of the spectra +locations. +.le +.ls o +Provide for a user identification list for the spectra locations. +This list shall be of the standard image cursor type to allow generation +of the list with the standard image cursor programs. +.le +.ls o +Determine and correct for a slowly varying background. +.le +.ls o +Reliably and accurately trace spectra in the presence of geometric +distortions (pincushion, s, shear, etc.). +.le + +Extract spectra by one of: +.ls a. +Strips of constant width about the located spectra. The width may be specified +to fractions of a pixel and the extraction will use fractional pixel +interpolation. +l +.le +.ls b. +Strips of width proportional to a Gaussian width parameter. +.le +.ls c. +Modeling to obtain estimates of the total luminosity. The estimate will +be the integral of the model. +.le +.ls d. +Summation of the data pixel values with fractional contributions of the +pixel value to the spectra based on modeling. +.le +.le +.ls o +An option shall be available to specify whether to ignore blank pixels +or use interpolated values. +.ls o +Programs shall be provided to produce data files which can be accessed +by one dimensional spectroscopic reduction routines. At a minimum +these formats shall include: +.ls a. +Reduction to an image file consisting of one line per extracted +spectrum +.le +.ls b. +The standard IIDS format available with the CYBER Multi-Aperture Plate +programs +.le +.le +.sh +3.2 Modeling Requirements + + The modeling of multi-spectra images, particularly in the case of +blended spectra, shall: +.ls (1) +Model blended spectra with sufficient reliability and robustness that +a reasonable solution is always obtained, though of possibly limited +usefulness. +.le +.ls (2) +The modeling shall provide estimates for the uncertainties in the fitted +parameters as a function of position along the spectrum. +.le +.ls (3) +Remove cosmic rays and other defective pixels by reference to the model. +.le +.ls (4) +Allow the transfer of a model solution for one image to another image. +.le +.ls (5) +Display numerically and graphically the data, the fitted model, and +the residuals. +.le +.sh +4. Program Specifications +.sh +4.1 Basic Programs + + The basic programs of the package are general purpose tools which +initialize a MULTISPEC data file and perform a single fundamental operation +on the data in the MULTISPEC data file. There is one data file associated +with each image. The data file is hidden from the user and so the user +need not be aware of the data file. +The data files are referenced only the image filename specified in the +program parameters. +The data files contain such information as a processing history, the +spectra positions and extracted luminosities, the model parameters (one +set for each spectra for each modelled image line (or swath), etc. +The programs generally are allowed to specify specific +lines, columns, and/or spectra on which to operate. +The line, column and spectra specifications are given as strings which +contain numbers separated by whitespace, commas, and the range indicator +"-". The script tasks +of section 4.2 will combine these basic programs to perform a general +multi-spectra extraction. + +.ks +.nf + ap_plate copy_params find_spectra convolve + fit_bckgrnd find_bckgrnd line_list model_extrac + model_fit model_image model_list sigma_extract + strip_extract to_iids to_image to_onedspec +.fi +.ke +.ks +.ls ap_plate +The information from an on-line data file containing descriptions of all +the aperture plates prepared at Kitt Peak is read to find a specified +aperture plate. The drilled aperture positions are correlated with the +spectra in the image to deduce relative wavelength offsets. The +identifications for the spectra as well as other auxiliary information +is recorded in the data file. +If no image file is specified then only the aperture +plate information is printed. This program is used in the MULTIAP_EXTRACT +program. This program is not essential to the operation of the MULTISPEC +package. + +.nf + Multi-Spectra image image = + Aperture plate plate = + (mode = ql) +.fi +.le +.ke +.ks +.ls The Background +The are two possibilities for dealing with the background. In the first +case, FIT_BCKGRND, the background will be fitted by polynomials and +the coefficients stored in the MULTISPEC data file. These coefficients +are then used by the other programs to estimate the background at the +spectra. The second option, FIND_BCKGRND, generates a background image in which +the spectra and other selected areas are set to blank pixels. Then a +general image interpolator is used fill in the blank pixels with background +estimates. The other MULTISPEC programs will then access this background +frame. The background frame image name will be stored in the MULTISPEC +data file and the image header. + +.ls fit_bckgrnd +Fit a background in a MULTISPEC image by a polynomial using pixels +not near the spectra and in the user specified swaths and columns. +The buffer distance is in pixels and refers to a minimum distance from +the center of any spectrum beyond which the background pixels are found. +Blank pixels are ignored in the background fit. Deviant pixels will be +rejected. + +.nf + Multi-Spectra image image = + Buffer from spectra buffer = 12 + Polynomial order order = 3 + Lines per swath (lines_per_swath = 32) + Swaths to fit (swaths = 1-1000) + Columns to fit (columns = 1-1000) + Rejection threshold (threshold = 5) + Print general diagnostics (verbose = no) + (mode = ql) +.fi +.le +.ls find_bckgrnd +The spectra within a buffer distance and specified areas are set to blank +pixels and the remaining pixels copied to a background image file. + +.nf + Multi-Spectra image image = + Background image background = + Buffer from spectra buffer = 12 + Lines to ignore (lines = ) + Columns to ignore (columns = ) + (mode = ql) +.fi +.le +.le +.ke +.ks +.ls convolve +A program will be provided to reduce either the extracted spectrum or +the modeled image to a common point-spread function. +.le +.ke +.ks +.ls copy_params +Create a MULTISPEC data file for a new image using +appropriate MULTISPEC parameters from an old image. +The old image must have been processed to find the spectra using FIND_SPECTRA +and possibly model fit. + +.nf + Old Multi-Spectra image old_image = + New Multi-Spectra image new_image = + (mode = ql) +.fi +.le +.ke +.ks +.ls find_spectra +Initially locate the spectra in a MULTISPEC image. +The positions of the spectra within the range of columns are determined +for the starting line and then the spectra are tracked within the +range of lines. The minimum separation +and minimum width would generally be set for a particular instrument. +If the automatic search is not used then a list of cursor positions is +read from the standard input. + +.nf + Multi-Spectra image image = + Automatic search auto = yes + Starting line start_line = + Minimum separation (min_sep = 1) + Minimum width (min_width = 1) + Averaging width (average = 32) + Lines to search (lines = 1-1000) + Columns to search (columns = 1-1000) + Print general diagnostics (verbose = no) + (mode = ql) +.fi +.le +.ke +.ks +.ls line_list +For the specified lines in the image print the image column +number, data value (possibly as a swath average), the model value at that +point (i.e. the sum of the model contributions from all the spectra), +the background value, and the residual. +Plotting scripts may be written using this routine to +show the quality of a model fit. + +.nf + Multi-Spectra image image = + Lines to list (lines = 1-1000) + (mode = ql) +.fi +.le +.ke +.ks +.ls model_extract +A previously fitted model is used to extract the spectra total luminosity +by apportioning the data values to spectra in the ratio indicated by the +model. If the clean option is specified then the model is used to detect +pixels which deviate from the model by a specified amount. +The model value replaces the deviant pixel in the extraction and, if specified, +also in the image file. + +.nf + Multi-Spectra image image = + Lines to extract (lines = 1-1000) + Clean spectra (clean = yes) + Cleaning threshold (threshold = 5) + Modify image (modify = yes) + Print general diagnostics (verbose = no) + (mode = ql) +.fi +.le +.ke +.ks +.ls model_fit +A specified model is iteratively fitted to the data in each of the specified +lines (or swaths) until the RMS residual fails to decrease. The models +are selected by a string. The possible values are + +.nf + (null string) - initialize the model + i - fit only the intensity scale + ip - fit the intensity scale and the position + ips1 - fit the intensity scale, position, and one parameter shape + ips2 - fit the intensity scale, position, and two parameter shape + ips3 - fit the intensity scale, position, and three parameter shape + ips4 - fit the intensity scale, position, and four parameter shape +.fi +These models will be combined in a script to search for the best fit. + +The initial shape parameters will generally be set by scripts for a +particular data reduction. + +.nf + Multi-Spectra image image = + Model type model = + Lines per swath (lines_per_swath = 32) + Swaths to model (swaths = 1-1000) + Initial shape1 (shape1 = .1 ) + Initial shape2 (shape2 = 0 ) + Initial shape3 (shape3 = 0 ) + Initial shape4 (shape4 = 5 ) + Print general diagnostics (verbose = no) + (mode = ql) +.fi +.le +.ke +.ks +.ls model_image +An image file of the fitted model is created. This image may then be displayed +or a residual image may be calculated and displayed. + +.nf + Multi-Spectra image image = + Model image model = + (mode = ql) +.fi +.le +.ke +.ks +.ls model_list +For the specified lines and spectra the model is listed. +The listing gives, for each spectra, +the spectrum number, the line number, the fitted position, +the estimated wavelength, the +extracted luminosity, the intensity scale, model width parameters, and +the background polynomial coefficients. This routine can be used in scripts +to plot the extracted spectra, the trend of width with wavelength, and so +forth. + +.nf + Multi-Spectra image image = + Lines to list (lines = 1-1000) + Spectra to list (spectra = 1-1000) + (mode = ql) +.fi +.le +.ke +.ks +.ls sigma_extract +A previously fitted model is used to extract the spectra luminosity +within a specified sigma of the peak. Because the model is not necessarily +a Gaussian the sigma is used to compute +the intensity ratio of the cutoff to the peak assuming a Gaussian profile +and then the data is extracted to the point the model intensity falls below that +cutoff. If the clean option is specified then the model is used to detect +pixels which deviate from the model by a specified amount. +The model value replaces the deviant pixel in the extraction and, if specified, +also in the image file. + +.nf + Multi-Spectra image image = + Sigma extraction width width = 1. + Lines to extract (lines = 1-1000) + Clean spectra (clean = yes) + Cleaning threshold (threshold = 5) + Modify image (modify = yes) + Print general diagnostics (verbose = no) + (mode = ql) +.fi +.le +.ke +.ks +.ls strip_extract +A strip of constant width about the spectra positions is extracted. +If cleanning is desired a smoothed estimate of the profile is obtained +by averaging a number of lines about the line to be cleaned. After fitting +for the intensity scale pixels are found which deviate from the profile by +a specified amount. +The profile value replaces the deviant pixel in the extraction and, +if specified, also in the image file. No prior modeling is required +to use this extraction routine. + +.nf + Multi-Spectra image image = + Strip extraction width width = 1. + Lines to extract (lines = 1-1000) + Clean spectra (clean = yes) + Cleaning threshold (threshold = 5) + Lines per profile average (averge_lines = 32) + Modify image (modify = yes) + Print general diagnostics (verbose = no) + (mode = ql) +.fi +.le +.ke +.ks +.ls to_iids +For a specified prefix, files of the form prefix.nn, where nn is a specified +spectra number, are created containing the extracted spectra for all +the specified image files. The format of the files is the IIDS format +developed for the CYBER Multi-Aperture Plate Extractions. + +.nf + Multi-Spectra image images = + IIDS filename prefix iids_file = + Spectra to format (spectra = 1-1000) + (mode = ql) +.fi +.le +.ke +.ks +.ls to_image +An image file containing one line of the extracted luminosities for each +specified spectra in the specified MULTISPEC image. + +.nf + Multi-Spectra image in_image = + Extracted spectra image out_image = + Spectra (spectra = 1-1000) + (mode = ql) +.fi +.le +.ke +.ks +.ls to_onedspec +The extractions are converted to an as yet to be specified format for +use in the ONEDSPEC reduction package. + +.nf + Multi-Spectra images images = + ONEDSPEC data file onedspec_file = + Spectra (spectra = 1-1000) + (mode = ql) +.fi +.le +.ke +.sh +4.2 General MULTISPEC CL Scripts + + The general MULTISPEC CL scripts perform a series of steps needed to +extract the spectra from a specified list of image files. These steps have +been found to generally perform the desired extraction task fully. + + +.nf + multiap_extract echelle_extract +.fi +.ks +.ls multiap_extract +The specified multi-aperture plate images are extracted. +If no starting solution image, one which has previously been extracted, +is specified then the script performs an automatic search for the +specified number of spectra. +Otherwise the solution from the starting image is used as the initial +model. The background is then determined. +This is followed by a series of fitting steps on swaths of data. +(For further details on the fitting steps see the Algorithms paper). +A MODEL_EXTRACT and cleaning follows. +Finally, the extraction is correlated with the specified aperture plate +using AP_PLATE. +If there was no starting image then this extraction becomes the +initial solution image. +Subsequent images are extracted starting from the initial solution image. + +.nf + Multi-Aperture images images = + Initial solution image initial = + Aperture plate number plate = + Number of spectra nspectra = + (mode = ql) +.fi +.le +.ke +.ks +.ls echelle_extract +The specified echelle images are extracted. +If no starting solution image, one which has previously been extracted, +is specified then the script performs an automatic search for the +specified number of orders. +Otherwise the solution from the starting image is used as the initial +starting point. The background is then determined. +Finally a STRIP_EXTRACT and cleaning is performed. +If there was no starting image then this extraction becomes the +initial solution image. +Subsequent images are extracted starting from the initial solution image. + +.nf + Echelle images images = + Initial solution image initial = + Number of orders norders = + Extraction width width = + (mode = ql) +.fi +.le +.sh +5. Outline of a MULTISPEC Reduction + + The following outline is for the reduction of a cryogenic camera +multi-aperture plate. All the programmer supplied default values are +used. + +.nf + (1) rcamera mtb, "ap165.", "s", "3-9" + (2) debias "ap165.*" + (3) multispec_flat "ap165.[36]", "ap165.flat" + (4) flat_divide "ap165.*", "ap165.flat" + (5) multiap_extract "ap165.*", "", 165, 50 + (6) to_onedspec "ap165.*", oned165 +.fi + +.ls (1) +The data is read from the observing tape(s) using RCAMERA. +The image files created are ap165.3, ap165.4, ..., ap165.9. This is +easily accomplished by using the filename prefix "ap165." in the RCAMERA +program. The raw images may be examined at this point on a display. +.le +.ls (2) +The images are debiased using DEBIAS with all the "ap165." files specified. +The debias program knows about the location of the bias strip for the +cryogenic camera. +.le +.ls (3) +A a flat field is created +using MULTISPEC_FLAT in which the desired quartz frames are specified +and a flat field image filename is defined. The created flat field +image may be examined on an image display if desired. +.le +.ls (4) +All the debiased images are divided by the flat field using FLAT_DIVIDE. +.le +.ls (5) +The script MULTIAP_EXTRACT is run in which the aperture plate number, +the number of spectra, and the image files to be extracted are specified. +The number of spectra is found by examining an image on an image display +or by plotting a cut across the spectra using a general image profile +program. +.le +.ls (6) +Finally, the extracted spectra are formatted for the ONEDSPEC package +using TO_ONEDSPEC with the extracted images specified. +.le +.endhelp diff --git a/noao/twodspec/multispec/doc/MSspecs_c.hlp b/noao/twodspec/multispec/doc/MSspecs_c.hlp new file mode 100644 index 00000000..848d589d --- /dev/null +++ b/noao/twodspec/multispec/doc/MSspecs_c.hlp @@ -0,0 +1,243 @@ + +.help multispec Nov82 "Multispec Specifications" +.ce +Comments on Multispec Package Specifications +.ce +November 8, 1983 + + + + The basic package structure and the decomposition of the package into +tasks looks good. The requirements for both general operators and canned +procedures are addressed well. I got the impression that you have a pretty +clear idea of what you want to do (which is the thing I am most looking for +when I read a specs document), but I confess to having to reread the document +several times to figure out what you have in mind. Your writing style is +very terse and leaves much up to the reader! + +Most of my comments have to do with details. These are presented in the +order in which they occurred while reading the document. These comments +apply only to the specs document. I have started going over the algorithms +paper, mostly when I could not understand a section of the specs document, +but I have not finished it yet. + +.sh +General Comments +.ls 4 +.ls (1) +When eventually we write the user documentation, the nomenclature +should be carefully explained up front. Users will tend to confuse +image lines and spectral lines, but there is little we can do about +that other than to make the distinction clear. The term "band" is +confusing because it normally refers to the third dimension of an +image and that is not how it is used here. A better term might be +"swath". In what follows I will continue to use the term band, but +it is definitely not too late to change. +.le +.ls (2) +It seems to me that the concept of a band or swath is a detail of how +the algorithm works and should not have such a prominent place in the +user interface to the package. Several of the routines require that +image coordinates be entered in units of band number and column. +This introduces an unnecessary coupling between two input parameters +and forces the user to convert from line number to band number. The +result will be that the user will be reluctant to change the number +of lines per band (I'll bet that you have kept this a constant in +using the prototype). My inclination would be to have the user enter +all coordinates in units of lines and columns, and have the program +select the nearest band depending on the band width parameter. +The band width could then be easily changed depending on the data, +without need to respecify the region of the image to be processed. +.le +.ls (3) +Routines all over the system will have an option for printing extra +information, i.e., a verbose mode of execution. I think we should +standardize on the name of this parameter. "Verbose" seems to me +more descriptive than "print", and is consistent with UNIX terminology. +.le +.le + +.sh +Pages 3,4 +.ls +.ls (1) +Functions for extracting spectra. I assume "strips of constant +width" means aperture sum out to a specified x-radius from the +center of a spectra. Can the radius be specified in fractional +pixels, and if so, does the routine do fractional pixel interpolation. +What happens if there are blank pixels in the aperture? + +If extraction is based on the model, I gather that you are still +summing data pixel values, using a weight for each spectra based +on the modeled contribution of each spectra to the data pixel. In +other words we are still taking an aperture sum, but with allowances +for crowding. This has the disadvantage that if we sum way out into +the wings, we will be adding noise to the aperture sum, degrading signal +to noise. + +Extraction based on integration of the model rather than +the data should be available as another extraction procedure; this may +yield better photometric results. I would eventually like to compare +the two approaches with artificial data. Also by integrating the model +there is no need to "clean" (I assume that deviant pixels are detected +and rejected when the model is fitted, or the model will not be +accurate). Blank pixels should be recognized and ignored when fitting +the model. +.le + +.ls (2) +I gather that all extracted spectra for an image are put into a single +imagefile. This is fine, even desirable, as long as it is ok if all +spectra share the same header, and as long as all we want to output +is intensity versus wavelength. If it is desired to also output the +signal to noise or whatever than another scheme may be needed. +.le +.ls (3) +The text file output form ('c'pg.4) should be engineered with the idea +that the user will take the data away in cardimage form. From the +description it sounds like there is one pixel (wavelength bin) per +line in the text file. This has its advantages, but is not what one +wants for a cardimage file, which always writes 80 chars per line. +Also, the detailed technical specs should give some details about +such a format; it is a major part of the user interface and people +will want to know what this format is going to look like. In a way +it is more important to specify formats like this than the calling +sequences of the tasks, because it is harder to change after the +package is released, and other program are written to read the +text format spectra. +.le +.ls (4) +To item 3.2 (2) (on uncertainty estimates) I would add "as a function +of position along the spectrum". +.le +.le + +.sh +4.1 Basic Programs +.ls +.ls (1) +Evidently there is a datafile associated with each image. What is +the function of the datafile? Is it transparent to the user? How +much is stored in the image header and how much in the datafile? +.le +.ls (2) +The distinction between "line_list" and "model_list" is confusing. +Does "line_list" print the sum of the models for all the spectra +a each column? Please specify the form of the output for this +procedure in more detail. The line_list and model_list procedures +are natural candidates for use with the "lists" utilities for +extracting columns, plotting one column against another, etc. I +could not tell whether or not this would work well from the info +given. +.le +.ls (3) +"ap_plate": "The identifications for the spectra ... is recorded." +Is recorded where? In the datafile? Is this information essential +to the operation of multispec, or is it merely passed on through +multispec? +.le +.ls (4) +"find_background": Might be more aptly named "fit_background". +I would expect "find" to mean find which regions of the image are +background and which are spectra. Find is spatial, fit is grayscale. + +We need to decide whether we want to specify polynomials in IRAF by +the order (0,1,2, etc.) or by the number of coefficients or terms. +It seems to me that people are most used to talking about second, +third, fifth etc. order polynomials and that we might better specify +polynomials with an "order" parameter rather than a "terms" param. + +Buffer radius or diameter? I would assume radius, but it is not +clear from the docs. What is being "searched"? Shouldn't that read +"bands to be fitted". The "colummns" parameter should permit a list +of ranges of columns; I couldn't tell whether this was the case +from the specs. Cursor input may be desirable here. + +Blank pixels should be detected and ignored when fitting the +background. Are deviant pixels detected and rejected? This is +generally a desirable option in a bkg fit. You may be able to +decompose this routine (internally) into a find_background and +a fit_background, making use of the Images background fitting +routines, though these generate an image as output rather than the +coeff of the fitted functions. I wuld guess that you are storing +the bkg coeff for each band in the datafile from the description, +and that the fit is strictly one-dimensional. + +If only a limited number of bands are fitted, what do you do about +the other bands if the bkg fit is one-dimensional? Is the user +req'd to use the same bands range when they do the extraction? +.le + +.ls (5) +"find_spectra". It is not clear how this routine uses cursor input. +Perhaps you should have a gcur type parameter. Reading cursor +coordinates from the standard input may be the way to go, but you +should explain how this is going to work. +.le +.ls (6) +"line_list". One output line per image line? One or more spectra +per output line? Output should be suitable for further processing +with the LISTS package utilities (i.e., getcol, and the graphics +utility which will plot or overplot lists). The specs should +specify the form of the output. +.le +.ls (7) +I assume that the extraction procedures extract spectra which +are put somewhere. Where, in the datafile? If the image is +to be cleaned, it would be safer to write a new output image, +or at least to rename the original. It is strange to have these +two quite different functions in the same module. +.le +.ls (8) +"model_fit". The range of modeling options is impressive, good +stuff. However, there must be something better than magic integer +numbers for specifying the model to be fitted. Perhaps the +strings "i, ip, ipw, ipw2, ipw3, ipw4", where 'i' is for intensity, +'p' for position, and 'w' for width. + +How are the "initial parameters" specified? +.le +.ls (9) +"model_list". Again, I can only guess from the description what the +output will look like. It sounds like it might be best to have +this routine print data for only one spectra at a time, particularly +if the lists package is to be used for analysis. It might be good +to have the line number in the output somewhere, especially if the +wavelength information is not available. +.le +.le + +.sh +4.2 Scripts +.ls +.ls (1) +It sounds like there is no easy alternative to an automatic search +for the line centers. This is best as long as it works, but the +users will want easy way to use the cursor available as an option. +A script such as this can easily use the line plot routine Images +to make a plot and generate a list of line centers, without even +requiring find_spectra to be able to access the cursor (and perhaps +it should not if the script can do it). The graphics cursor should +be used here rather than the image cursor. +.le +.le + +.sh +5. Example +.ls +.ls (1) +The rcamera example is in error. Rcamera, as implemented, has only +three query mode params, while you show four in the example. +I believe the ranges string should be quoted and should be the second +argument. + +The last command should be "to_onedspec", not "onedspec". +.le +.ls (2) +5.(5): It seems strange to make the user manually count 50 spectra +by examining a plot. If the program automatically finds centers, +this should not be necessary; if the user interactively marks centers, +it is not necessary. +.le +.le +.endhelp diff --git a/noao/twodspec/multispec/doc/findpeaks.hlp b/noao/twodspec/multispec/doc/findpeaks.hlp new file mode 100644 index 00000000..f6118281 --- /dev/null +++ b/noao/twodspec/multispec/doc/findpeaks.hlp @@ -0,0 +1,88 @@ +.help findpeaks Jul84 noao.twodspec.multispec +.ih +NAME +findpeaks -- Find peaks in a multi-spectra image +.ih +USAGE +findpeaks image lines contrast +.ih +PARAMETERS +.ls image +Image to be searched. +.le +.ls lines +Sample image lines in which the peaks are to be found. +.le +.ls contrast +Maximum contrast between the highest peak and the lowest peak. +.le +.ls separation = 5 +Minimum separation in pixels between acceptable peaks. +.le +.ls edge = 0 +Minimum distance in pixels to the edge of the image for acceptable peaks. +.le +.ls threshold = 0. +The minimum acceptable peak pixel value. +.le +.ls min_npeaks = 1 +Minimum number of peaks to be found. It is an error for fewer than +this number of peaks to be found. +.le +.ls max_npeaks = 1000 +Maximum number of peaks to be found. If more than this number of peaks +is found then only the those with the highest peak values are accepted. +.le +.ls columns = '*' +Columns to be searched. +.le +.ls naverage = 20 +Number of image lines around the sample line to be averaged before +finding the peaks. +.le +.ls debug = no +Print detailed information on the progress of the peak finding algorithm. +.le +.ih +DESCRIPTION +For each specified sample image line the number of peaks and their column +positions in the image are determined. +The number of peaks and their positions are assumed to correspond to points +along the spectra. This information is entered in the MULTISPEC database. + +The \fInaverage\fR image lines about the specified sample line are first +averaged. The local maxima in the average line are then located +in the specified columns more than the minimum distance from the edge of the +image. A minimum peak pixel value cutoff is determined as the maximum of +the specified \fIthreshold\fR and \fIcontrast\fR times the largest peak pixel +value. All local maxima with pixel values below the cutoff are rejected. +Next all peaks with separations less than \fIseparation\fR from a stronger +peak are rejected. Finally, if there are more than \fImax_npeaks\fR remaining +only the \fImax_npeaks\fR strongest peaks are accepted. If fewer +than \fImin_npeaks\fR are found then the task quits with an error. + +If the number of spectra has been previously determined, such as by an earlier +use of \fBfindpeaks\fR, then it is an error if a different number of +peaks is found. +.ih +EXAMPLES +The parameters of this task provide a great deal of flexibility in +automatically determining the number and positions of the peaks. +The most automatic method just uses the contrast to limit the acceptable +peaks: + + cl> findpeaks image.db 1 .1 + +However, if the number of spectra in the image is known: + + cl> findpeaks image.db 1 0 min=10 max=10 + +or if a threshold is known: + + cl> findpeaks image.db 1 0 threshold = 1000 + +For a noisy image the separation parameter can be set to eliminate spurious +noise peaks near the peaks to be found: + + cl> findpeaks image.db 1 .1 sep=20 +.endhelp diff --git a/noao/twodspec/multispec/doc/fitfunc.hlp b/noao/twodspec/multispec/doc/fitfunc.hlp new file mode 100644 index 00000000..09510bb7 --- /dev/null +++ b/noao/twodspec/multispec/doc/fitfunc.hlp @@ -0,0 +1,73 @@ +.help fitfunction Jul84 noao.twodspec.multispec +.ih +NAME +fitfunction -- Fit a function to the spectra parameter values +.ih +USAGE +fitfunction image +.ih +PARAMETERS +.ls image +Image in which the parameter values are to be fitted. +.le +.ls parameter = "x0" +Parameter to be fit. The legal minimum match abbreviated parameters +are x0, s0, s1, s2. +.le +.ls lines = "*" +Sample image lines to be used in the function fit. +.le +.ls spectra = "*" +Spectra for which the parameters are to be fit. +.le +.ls function = "interpolation spline" +Fitting function to be used. The function is specified as a string +which may be minimum match abbreviated. The functions currently available +are: +.ls interpolation spline +Interpolation spline of specified order. +.le +.ls smoothing spline +Smoothing spline of specified order and number of polynomial pieces. +.le +.le +.ls spline_order = 4 +Order of the fitting spline. The order must be even. +The minimum value is 2 and maximum value is determined from the number of +sample lines in the fit. +.le +.ls spline_pieces = 1 +The number of polynomial pieces in a smoothing spline. +The minimum value is 1 and the maximum value is determined from the number of +sample lines in the fit. +.le +.ih +DESCRIPTION +A function is fit to the parameter values previously determined at the sample +lines for each spectrum. The function coefficients are stored in the +database and the fitted values replace the original values at all the sample +lines (not just the sample lines used in the fit). The type of function, +the parameter to be fitted, the sample lines used in the fit, and the +spectra to be fitted are all selected by the user. The function is +extrapolated to cover all image lines. + +The values of the function fit at arbitrary image lines may be listed +with \fBmslist\fR. +.ih +EXAMPLES +The extraction of the spectra requires that a fitting function be +determined for the spectra positions. This is done by: + + cl> fitfunction image + +To smooth the parameter "s0" in model \fIgauss5\fR with a cubic spline +and leave out a bad point at sample line 7: + +.nf + cl> fitfunction image parmeter=s0 function=smooth \ + >>> lines="1-6,8-" +.fi +.ih +SEE ALSO +mslist +.endhelp diff --git a/noao/twodspec/multispec/doc/fitgauss5.hlp b/noao/twodspec/multispec/doc/fitgauss5.hlp new file mode 100644 index 00000000..bcb37276 --- /dev/null +++ b/noao/twodspec/multispec/doc/fitgauss5.hlp @@ -0,0 +1,148 @@ +.help fitgauss5 Jul84 noao.twodspec.multispec +.ih +NAME +fitgauss5 -- Fit spectra profiles with five parameter Gaussian model +.ih +USAGE +fitgauss5 image start +.ih +PARAMETERS +.ls image +Image to be modeled. +.le +.ls start +Starting sample line containing the initial model parameters. +.le +.ls lower = -10 +Lower limit for the profile fit relative to each spectrum position. +.le +.ls upper = 10 +Upper limit for the profile fit relative to each spectrum position. +.le +.ls lines = "*" +Sample image lines to be fit. +.le +.ls spectra = "*" +Spectra to be fit. +.le +.ls naverage = 20 +Number of data lines to be averaged about each sample image line before +model fitting. +.le +.ls factor = 0.05 +The model fit to each line is iterated until the RMS error between the +model line and the data line improves by less than this factor. +.le +.ls track = yes +Track the model solution from the starting line to the other sample lines? +.le +.ls algorithm = 1 +Parameter fitting algorithm to use. Legal values are 1 and 2. +.le +.ls fit_i0 = yes +Fit the profile scale parameters i0? +.le +.ls fit_x0 = yes +Fit the spectra position parameters x0? +.le +.ls fit_s0 = yes +Fit the spectra shape parameters s0? +.le +.ls fit_s1 = no +Fit the spectra shape parameters s1? +.le +.ls fit_s2 = no +Fit the spectra shape parameters s2? +.le +.ls smooth_s0 = yes +Fit a smoothing spline to the shape parameters s0 after each iteration? +.le +.ls smooth_s1 = yes +Fit a smoothing spline to the shape parameters s1 after each iteration? +.le +.ls smooth_s2 = yes +Fit a smoothing spline to the shape parameters s2 after each iteration? +.le +.ls spline_order = 4 +Order of the smoothing spline to be fit to the shape parameters. +.le +.ls spline_pieces = 3 +Number of polynomial pieces for the smoothing spline. +.le +.ls verbose = no +Print general information about the progress of the model fitting. +.le +.ih +DESCRIPTION +The spectra profiles in the interval (\fIlower, upper\fR) about each +spectrum position are fit with a five parameter Gaussian model for +the specified sample lines of the image. For a description of +the model see \fBgauss5\fR. The model fitting is performed using +simultaneous linearized least squares on the selected model profile +parameters as determined by the \fIalgorithm\fR for the specified +\fIspectra\fR. The parameter fitting technique computes correction +vectors for the parameters until the RMS error of the model image line +to the data image line, which is an average of \fInaverage\fR lines +about the sample line, improves by less than \fIfactor\fR. +A solution which increases the RMS error of the model is not allowed. + +If the parameter \fItrack\fR is yes then the initial model parameters are +those given in the database for the sample line \fIstart_line\fR. From +this starting point the model parameters are iterated to a best fit at +each specified sample line and then the best fit is used as the starting +point at the next line. The tracking sequence is from the starting line +to the last line and then, starting again from the starting line, to +the first line. Note that the model parameters, including the starting +spectra positions, need be set only at the starting line. + +If \fItrack\fR is no then each specified sample line is fitted independently +from the initial model parameters previously set for that line. This option +is used to add additional parameters to the model after an +initial solution has been obtained or to refit a new image whose database +was created as a copy of the database of a previously fit image. + +The shape parameters s0, s1, and s2 can be smoothed by fitting a spline of +specified \fIorder\fR and number of spline pieces, \fInpp\fR to the +parameters as a function of spectra position. +The smoothing is performed after each iteration and before +computing the next RMS error. The smoothing is a form of local constraint +to keep neighboring spectra from having greatly different shapes. +The possibility of such erroneous solutions being obtained is present in +very blended data. + +In \fIverbose\fR mode the RMS errors of each iteration are printed on the +standard output. + +The selection of the parameters to be fit and the order in which they are +fit is determined by \fIalgorithm\fR. These algorithms are: + +.ls 4 1 +This algorithm fits the selected parameters (\fIfit_i0, fit_x0, +fit_s0, fit_s1, fit_s2\fR) for the selected \fIspectra\fR simultaneously. +.le +.ls 4 2 +This algorithm begins by fitting the parameters i0, x0, and s0 +simultaneously. Note that the values of s1 and s2 are used but are +kept fixed. Next the parameters s0 and s1 (the shape) are fit simultaneously +keeping i0, x0, and s2 fixed followed by fitting i0 and x0 while +keeping s0, s1, and s2 (the shape) fixed. If either of these fits +fails to improve the RMS then the algorithm terminates. +Also, if after the two steps (the fit of s0 and s1 followed by the fit +of i0 and x0), the RMS of the fit has not improved by more than the +user specified factor the algorithm also terminates. This algorithm has been +found to be the best way to fit highly blended spectra. +.le +.ih +EXAMPLES +The default action is to fit Gaussian profiles to the spectra and trace +the fit from the starting line. An example of this is: + + cl> fitgauss5 image 1 + +To fit heavily blended spectra with the four parameter model (i0, x0, s0, s1): + + cl> fitgauss5 image 1 algorithm=2 +.ih +SEE ALSO +findspectra +.endhelp diff --git a/noao/twodspec/multispec/doc/modellist.hlp b/noao/twodspec/multispec/doc/modellist.hlp new file mode 100644 index 00000000..70e95ce4 --- /dev/null +++ b/noao/twodspec/multispec/doc/modellist.hlp @@ -0,0 +1,52 @@ +.help modellist Jul84 noao.twodspec.multispec +.ih +NAME +modellist -- List data and model pixel values +.ih +USAGE +modellist image lines +.ih +PARAMETERS +.ls image +Image whose model is to be listed. +.le +.ls lines +Sample lines to be listed. +.le +.ls model = "gauss5" +Profile model to be used to create the model line. +The only model currently defined is \fIgauss5\fR. +.le +.ls columns = "*" +Image columns to be listed. +.le +.ls naverage = 20 +The number of image lines to be averaged to form the data values. +.le +.ls lower = -10 +Lower limit of model profiles measured in pixels from the spectra centers. +.le +.ls upper = 10 +Upper limit of model profiles measured in pixels from the spectra centers. +.le +.ih +DESCRIPTION +The model of the image for the selected sample \fIlines\fR +are used to generate model image lines. Only the model \fIgauss5\fR is +currently available. The output format is column, sample line, image pixel +value, and model pixel value. The image pixel data are formed by averaging +\fInaverage\fR lines about the sample lines. +.ih +EXAMPLES +To list the image and model pixel values for the first sample line after +fitting the \fIgauss5\fR model with \fBfitgauss5\fR: + + cl> modellist image 1 >outputlist + +The list file \fIoutputlist\fR can be used with the \fBlists\fR and +\fBplot\fR packages to graph the image and model lines or to compute +and graph residuals. +.ih +SEE ALSO +newimage +.endhelp diff --git a/noao/twodspec/multispec/doc/msextract.hlp b/noao/twodspec/multispec/doc/msextract.hlp new file mode 100644 index 00000000..fa361b38 --- /dev/null +++ b/noao/twodspec/multispec/doc/msextract.hlp @@ -0,0 +1,172 @@ +.help msextract Jul84 noao.twodspec.multispec +.ih +NAME +msextract -- Extract spectra from a multi-spectra image +.ih +USAGE +msextract image output +.ih +PARAMETERS +.ls image +Image to be extracted. +.le +.ls output +Filename for the three dimensional image to be created containing the +extracted spectra. +.le +.ls lower = -10 +Lower limit of the integral for integrated spectra or the first column of the +strip spectra. It is measured in pixels from the spectrum center +defined by the position function in the MULTISPEC database. +.le +.ls upper = 10 +Upper limit of the integral for integrated spectra or (approximately) the +last column of the strip spectra. It is measured in pixels from the +spectrum center defined by the position function in the MULTISPEC database. +.le +.ls spectra = "*" +Spectra to be extracted. +.le +.ls lines = "*" +Image lines to be extracted. +.le +.ls ex_model = no +Extract model spectra fit to the image spectra? +.le +.ls integrated = yes +Extract integrated spectra? +.le +.ls unblend = no +Correct for blending in the extracted spectra? +.le +.ls clean = yes +Replace bad pixels with model values? The following parameters are used: +.ls nreplace = 1000. +Maximum number of pixels to be replaced per image line when cleaning with +model \fIgauss5\fR or maximum number of pixels to be replaced per spectrum when +cleaning with model \fIsmooth\fR. +.le +.ls sigma_cut = 4. +Cleaning threshold in terms of sigma of the fit. +.le +.ls niterate = 1 +Maximum number of cleaning iterations per line when cleaning with model +\fIgauss5\fR. +.le +.le +.ls model = "smooth" +Choice of \fIgauss5\fR or \fIsmooth\fR. Minimum match abbreviation is +allowed. This parameter is required only if \fIex_model\fR = yes +or \fIclean\fR = yes. +.le +.ls naverage = 20 +Number of lines to be averaged in model \fIsmooth\fR. +.le +.ls fit_type = 2 +Model fitting algorithm for model \fIgauss5\fR. +.le +.ls interpolator = "spline3" +Type of image interpolation function to be used. +The choices are "nearest", "linear", "poly3", "poly5", and "spline3". +Minimum match abbreviation is allowed. +.le +.ls verbose = no +Print verbose output? +.le +.ih +DESCRIPTION +The MULTISPEC database describing the spectra positions and shapes +is used to guide the extraction of the spectra in the multi-spectra image. +The user selects the \fIspectra\fR and image +\fIlines\fR to be extracted and whether to extract integrated or strip spectra. +In addition options are available to extract model spectra, replace bad +pixels by model spectra values, and correct for blending of the spectra. +The \fIoutput_file\fR three dimensional +image consists of one band (the third dimension) per extracted spectrum, +the extracted lines (the second dimension) and either one column for +the integrated luminosity or the number of columns in the extracted strip. + +Integrated spectra (\fIintegrated\fR = yes) are extracted by summing +the pixel or model values over the specified limits \fIlower\fR and \fIupper\fR +measured relative to the spectra centers defined by the position functions in +the database. Partial pixel sums are used at the endpoints. + +Strip spectra (\fIintegrated\fR = no) are extracted by image interpolation +of the image line or model profiles to obtain a line of values for +each spectrum and for each image line. The length of the strip is the +smallest integer containing the interval between \fIlower\fR and \fIupper\fR. +The strips for each spectrum are aligned so that the first column is a distance +\fIlower\fR from the spectrum center as given by the position function in the +database. + +If \fIex_model\fR = yes, \fIunblend\fR = yes, or \fIclean\fR = yes model +spectra are fit to the spectra in the image. There are two models: +a five parameter Gaussian profile called \fIgauss5\fR and profiles obtained +by averaging \fInaverage\fR image lines surrounding the image line being +modeled called \fIsmooth\fR. The model is selected either when the parameter +\fIunblend\fR = yes or with the parameter \fImodel\fR. If \fIunblend\fR = yes +then the model is \fIgauss5\fR regardless of the value of \fImodel\fR. + +When \fIex_model\fR = yes the effect is to substitute model spectra for the +image spectra in the output extraction image. + +When \fIclean\fR = yes pixels with large residuals from the model are +detected and removed from the model fit. The selected model is +fit to the pixels which are not in the bad pixel list (not yet implemented) +and which have not been removed from the model fit. The sigma of the fit +is computed. Deviant pixels are detected by comparing them to the model +to determine if they differ by more than \fIsigma_cut\fR times the sigma. +The model fit is iterated, removing deviant pixels at each iteration, until +no more pixels are found deviant or \fInreplace\fR pixels have been found. +The pixels removed or in the bad pixel list are then replaced with +model values. (To clean an image with this algorithm see \fBnewimage\fR.) + +There are some technical differences in the model fitting and cleaning +algorithms for the two models. In model \fIsmooth\fR +the fit for the profile scale factors is done independently for each spectrum +and automatically corrected when a bad pixel is detected. This fitting process +is fast and rigorous. The parameter \fInreplace\fR in this model refers to +the maximum number of pixels replaced \fIper spectrum\fR. + +In model \fIgauss5\fR, however, the profile scale factors are fit +to the entire image line (hence its ability to fit blended spectra). +There are two fitting algorithms; a rigorous simultaneous fit +and an approximate method. The simultaneous fit is selected when +\fIfit_type\fR = 1. This step is relatively slow. The +alternative method of \fIfit_type\fR = 2 sets the scale factor for each +spectrum by taking the median scale, where scale = data / model profile, +for the three pixels nearest the center of the profile. The median +minimizes the chance of a large error due to a single bad pixel. This +scale may be greatly in error in the case of extreme blending but is also +quite fast; the extraction time is reduced by at least 40%. +The steps of profile fitting and deviant pixel detection are alternated +and the maximum number of iterations through these two steps is +set by \fIniterate\fR. The default of 1 means that the model fitting is not +repeated after detecting deviant pixels. + +When \fIunblend\fR = yes the \fIgauss5\fR model +is fitted to the image spectra (including possible cleaning). +The relative contributions to the total image pixel value from each of the +blended spectra are determined from the model and applied toward either the +integrated or strip spectra. If \fIex_model\fR = yes then this option has +no effect other than to force the selection of model \fIgauss5\fR. + +The option \fIverbose\fR is used to print the image lines being extracted +and the number of pixels replaced by the cleaning process. +.ih +EXAMPLES +To extract all the integrated spectra from all the image lines: + + cl> msextract image image.ms + +To extract model strip spectra: + + cl> msextract image image.ms ex_model=yes int=no + +To extract integrated spectra without any modeling: + + cl> msextract image image.ms clean=no +.ih +SEE ALSO +newimage +.endhelp diff --git a/noao/twodspec/multispec/doc/mslist.hlp b/noao/twodspec/multispec/doc/mslist.hlp new file mode 100644 index 00000000..461b52b4 --- /dev/null +++ b/noao/twodspec/multispec/doc/mslist.hlp @@ -0,0 +1,77 @@ +.help mslist Jul84 noao.twodspec.multispec +.ih +NAME +mslist -- List entries in a MULTISPEC database +.ih +USAGE +mslist image keyword lines spectra +.ih +PARAMETERS +.ls image +Image whose MULTISPEC database entries are to be listed. +.le +.ls keyword +Keyword for the database entry to be listed. The keywords are: +.ls header +List general header information. +.le +.ls comments +List the comments. +.le +.ls samples +List the sample image lines. +.le +.ls x0 +List the spectra positions for the specified sample lines and spectra. +.le +.ls i0 +List the model profile scales for the specified sample lines and spectra. +.le +.ls s0, s1, or s2 +List the gauss5 model shape parameter s0, s1, or s2 for the specified sample +lines and spectra. +.le +.ls gauss5 +List the gauss5 model parameters x0, i0, s0, s1, and s2 for the specified +sample lines and spectra. +.le +.ls x0 spline +List the spline evaluation of the spectra positions for the specified +image lines and spectra. +.le +.ls s0 spline, s1 spline, or s2 spline +List the spline evaluation of the gauss5 model shape parameters s0, s1, or s2 +for the specified image lines and spectra. +.le +.le +.ls lines +Lines to be listed. For the entries x0, i0, s0, s1, s2, and gauss5 the +lines refer only to the sample image lines. For the spline entries the +lines refer to the image lines at which the spline is to be evaluated. +.le +.ls spectra +Spectra to be listed. +.le +.ls titles = no +Print additional titles? +.le +.ih +DESCRIPTION +This task is a general MULTISPEC database listing tool. A keyword is selected +and the referenced data is listed. Some entries require the specification of +the desired sample or image lines and the desired spectra. +.ih +EXAMPLES +To list the spectra positions for spectrum 3 at all the sample lines: + + cl> mslist image x0 "*" 3 + +To list the model profile scale parameter for sample line 1: + + cl> mslist image i0 1 "*" + +To list the gauss5 model parameters for spectra 2 and 3 and sample lines 5 +and 7: + + cl> mslist image gauss5 "5,7" "2-3" titles+ +.endhelp diff --git a/noao/twodspec/multispec/doc/msplot.hlp b/noao/twodspec/multispec/doc/msplot.hlp new file mode 100644 index 00000000..f08eac1b --- /dev/null +++ b/noao/twodspec/multispec/doc/msplot.hlp @@ -0,0 +1,44 @@ +.help msplot Oct85 noao.twodspec.multispec +.ih +NAME +msplot -- Plot data and model image line +.ih +USAGE +msplot image line +.ih +PARAMETERS +.ls image +Image to be plotted. +.le +.ls line +The image line to be plotted. Actually the nearest sample line will be +plotted. +.le +.ls naverage = 20 +Number of image lines to average about the specified line. +.le +.ls lower = -10., upper = 10. +Limits of the model profiles relative to the center of each profile. +.le +.ls graphics = "stdgraph" +Graphics output device. +.le +.ls cursor = "" +Graphics cursor input. If a file is given then the cursor input is taken +from the file. If no file is given then the standard graphics cursor will +be used. +.le +.ih +DESCRIPTION +A line of image data and the profile model for the line is graphed. +The model is graphed with a dashed line. The graph may be then expanded, +manipulated, and printed with the standard cursor mode commands. +.ih +EXAMPLES +To plot the model fit for image sample for image line 400: + + cl> msplot sample 400 +.ih +SEE ALSO +modellist +.endhelp diff --git a/noao/twodspec/multispec/doc/msset.hlp b/noao/twodspec/multispec/doc/msset.hlp new file mode 100644 index 00000000..689e525a --- /dev/null +++ b/noao/twodspec/multispec/doc/msset.hlp @@ -0,0 +1,104 @@ +.help msset Jul84 noao.twodspec.multispec +.ih +NAME +msset -- Set entries in a MULTISPEC database +.ih +USAGE +msset image keyword value +.ih +PARAMETERS +.ls image +Image in which the MULTISPEC database entries are to be modified or initialized. +.le +.ls keyword +Keyword for the database entry to be set. The keywords are: +.ls nspectra +Set the number of spectra in the header. +.le +.ls comments +Add comments lines to the database comment block. +.le +.ls x0 +Set the spectra positions for the specified sample lines and spectra. +.le +.ls i0 +Set the model profile central intensities for the specified sample lines +and spectra. +.le +.ls s0, s1, or s2 +Set the gauss5 model shape parameter s0, s1, or s2 for the specified sample +lines and spectra. +.le +.le +.ls value +Value to be used for value input. +.le +.ls lines = "*" +Sample lines to be affected by value input. +.le +.ls spectra = "*" +Spectra to be affected by value input. +.le +.ls read_list = no +If yes use list input and if no use value input. +.le +.ls list = "" +List for list input. See the description below for the appropriate format. +.le +.ih +DESCRIPTION +The entries in a MULTISPEC database associated with the image +are modified or initialized. +The parameters \fIimage\fR, \fIkeyword\fR, and \fIread_list\fR +determine the database to be operated upon, the database entry to +be set, and the input type. There are two forms of input; +list input and value input. +The input type is selected by the boolean parameter +\fIread_list\fR. For list input the parameter \fIlist\fR +is used and for value input the parameter \fIvalue\fR and +possibly the parameters \fIlines\fR and \fIspectra\fR are used. +The required parameters and input formats for the different keywords +are outlined below. +.ls nspectra +For list input the list format is the number of spectra and +for value input the \fIvalue\fR parameter is the number of spectra. +.le +.ls comments +For list input the list format is lines of comments and for value +input \fIvalue\fR parameter is a comment string. +.le +.ls x0, i0, s0, s1, s2 +For list input the list format is sample line, spectrum number, and +parameter value +and for value input \fIlines\fR is a range string selecting the +sample lines to be affected, \fIspectra\fR is a range string selecting +the spectra to be affected, and \fIvalue\fR is the value to be set for all +the selected lines and spectra. +.le +.ih +EXAMPLES +To add several comments to the database by query: + +.nf + cl> msset image "comments" read_list+ + Input list> First comment here. + Input list> Second comment here. + Input list> <eof> +.fi + +where <eof> is the end of file character terminating the list. +To set the value of s0 to 1 for all the spectra in sample line 1: + + cl> msset image "s0" 1 + +To set the spectra positions from a list: + + cl> msset image "x0" read_list+ list=positionlist + +To add a single comment such as in a script: + + cl> msset image "comments" "Comment here." +.ih +SEE ALSO +findspectra mslist +.endhelp diff --git a/noao/twodspec/multispec/doc/multispec.ms b/noao/twodspec/multispec/doc/multispec.ms new file mode 100644 index 00000000..cc17352e --- /dev/null +++ b/noao/twodspec/multispec/doc/multispec.ms @@ -0,0 +1,532 @@ +.EQ +delim $$ +.EN +.TL +The Multi-Spectra Extraction Package (multispec) +.AU +Francisco Valdes +.AI +IRAF Group +.K2 +October 1984 +.NH +Introduction +.PP +This document provides an introduction and overview of the multi-spectra +extraction package \fBmultispec\fR. Detailed descriptions and usage +information for the tasks of the package are available in the manual +pages. The tasks in the package are: + +.TS +center; +n. +findpeaks \&- Find the peaks +fitfunction \&- Fit a function to the spectra parameter values +fitgauss5 \&- Fit spectra profiles with five parameter Gaussian model +modellist \&- List data and model pixel values +msextract \&- Extract spectra +mslist \&- List entries in a MULTISPEC database +msplot \&- Plot a line of image and model data +msset \&- Set entries in a MULTISPEC database +newextraction \&- Create a new MULTISPEC extraction database +newimage \&- Create a new multi-spectra image +.TE + +.PP +The \fBmultispec\fR package is a subpackage of the \fBtwodspec\fR package. +It provides tools to locate, model, clean, correct for blending, +and extract integrated or strip spectra from two dimensional, multi-spectra +images. These tools may be used directly or combined in scripts to +extract specific types of spectra or spectra from specific instruments. +Examples of the latter usage are the tasks in the image reduction package +\fBcryomap\fR. +.PP +The extraction of spectra consists of locating pixels along each +image line which intersect the spectra and recording either the sum of +the pixels, \fIintegrated spectra\fR (some times referred to as +one-dimensional spectra), or the set of pixels, +\fIstrip spectra\fR, for each line and for each spectrum as output. +The size and limits of the intersection region are specified by the +user relative to the centers of the spectra. +The locations of the spectra in each image line are determined separately +so that the spectra need not be aligned along the columns of the image nor +be perfectly straight. However, since the extraction is done by image line, +if the spectra are not aligned with the columns then the spectral resolution +will be decreased. If the spectra are aligned with the image lines then +the image should be rotated or transposed using \fBimtranspose\fR. +.PP +The \fBmultispec\fR extraction produces three dimensional images with +one image band (the third dimension) for each extracted spectrum +and one line (the second dimension) for each extracted image line. +For integrated spectra there is only one column +while for strip spectra, the number of columns is equal to the extraction +strip width. The strips are aligned to the same positions relative to the +spectra centers by image interpolation. If desired the output extractions can +be reformated in a variety of ways. +.PP +In addition to direct extraction of the image spectra the \fBmultispec\fR +package provides for modeling the spectrum profiles. The model +may be extracted instead of the image spectra as either integrated or +strip spectra. The model may be used to correct for blending of the spectra +and to detect and replace bad pixels. The cleaning replaces data pixels which +are discrepant from the model by the model values. +.PP +The modeling and cleaning features of the \fBmultispec\fR package can also +be used for creating new multi-spectra images. In other words a new +image is created containing cleaned or model spectra for selected +lines. +.PP +Section 2 gives an overview of the \fBmultispec\fR package and the extraction +process. The next section briefly describes the tasks in the package. +This is followed by a description of the extraction database. +The final section defines the model profiles used in the \fBmultispec\fR +package. +.NH +Overview of the Multispec Package and the Extraction Process +.PP +The \fBmultispec\fR package consists of general and flexible tools +for creating and manipulating databases which describe multi-spectra +images. The contents of the databases are described in a later section. +Each database is associated with a particular image and is referenced +through the image name. The first positional argument in all the +\fBmultispec\fR tasks is an image. In the current version of the package each +database exists as a separate binary file with a filename formed by adding +the extension '.db' to the image name. Note, however, that this +need not be the case in the future. +.PP +The organization of the package as a set of tools operating on a database +allows room for the package to evolve. Different algorithms may be +designed for different types of multi-spectra images by using combinations +of the existing tools and by adding new tools. The discussion below +points out areas where new tasks might be added as well as citing the +applicable existing tasks. +.PP +The extraction of spectra from a multi-spectra image consists of two +basic steps; determining the locations of the spectra in the image and +extracting the spectra. The positions of the spectra in a multi-spectra +image are determined at a set of "sample" image lines. These positions +are used to fit an interpolation function defining the spectrum positions +at all the image lines. This function is then used in the extraction of +the spectra. +.PP +The sample image lines are chosen by the user when the database is first +created by the task \fBnewextraction\fR. An exception to this is when +a template image is used (discussed below). However, in this case the +sample image lines are still those chosen by the user when the template +image database was created. The sample image lines may consist of +anywhere from one image line to all the image lines. The purpose +of the sample lines is to sample the image often enough to follow changes +in the positions and shapes of the spectra but to still minimize the +time spent in finding the spectra and the size of the database. The choice +of sample lines also depends on the algorithm used to determine the +positions of the spectra; a large number of sample +lines for a fast, approximate method and a smaller number of lines +for a complex and accurate method. For example, in order to deal with +very blended spectra the task \fBfitgauss5\fR provides a sophisticated +model fitting algorithm. This technique is computationally slow and, so, +the user should not choose too many sample lines. +.PP +After the database has been created the minimum information needed +for extraction is the spectrum positions at the sample lines. There +are many ways in which the positions may be determined. Some +possibilities are listed below. + +.IP (1) +Enter the spectrum positions from a list using \fBmsset\fR. The +list might be generated from a graphics/cursor task. +This is method is very time consuming when the number of spectra and +the number of images are large. +.IP (2) +Determine the spectrum positions automatically by finding the peaks in +each sample image line. The task \fBfindpeaks\fR performs this function. +.IP (3) +Determine the spectrum positions at just one sample image line +using either (1) or (2) and trace the spectra by a fast and refined +peak finding method. Such a task is desirable but is not a part of the +current package. +.IP (4) +Determine the spectrum positions at just one sample image line +using either (1) or (2) and trace the spectra by fitting model +spectrum profiles. The task \fBfitgauss5\fR does this using +the model gauss5 described in section 5. Additional model fitting +tasks can be added as needed. +.IP (5) +Use the positions determined for a previous image and, if necessary, +refine the positions. \fBFitgauss5\fR is used to +refine the spectrum positions at each sample line independently. + +.PP +Several position finding algorithms may be used in stages to achieve +the degree of accuracy required by the user. +Thus, the first position determinations may be relatively crude and +then, if needed, more sophisticated methods may be applied to refine the +positions. The task \fBfindpeaks\fR is a crude peak finder. The positions +are only determined to the nearest pixel. The task \fBfitgauss5\fR is +a sophisticated model fitting techique which is used after \fBfindpeaks\fR +first determines the approximate positions of the spectra. +.PP +The determination of the spectra locations may be performed independently +at each sample line as in (1) and (2) above or the spectra locations may +be traced starting from one sample line as in (3) and (4). The second method +is preferable. Generally, \fBfindpeaks\fR is used at only one sample line +to initially determine the number and approximate locations of the spectra. +\fBFitgauss5\fR then fits model gauss5 to the spectrum profiles and +the model solution is used at the next sample line as the starting +point for the next model fit. In this manner the positions of +the spectra are determined at the other sample image lines. +.PP +The results of the peak finding and profile fitting are improved +by using an average of many image lines about the sample image line rather +than just the sample image line by itself. Both \fBfindpeaks\fR and +\fBfitgauss5\fR have this ablility. +.PP +It is often the case that several multi-spectra images have essentially +the same format; i.e. the same image size, the same number of spectra, +and the same positions (either approximately or identically). +Commonly, one of the images is used for calibrations and has strong, +high signal-to-noise spectra while the other images have weaker spectra. +In this case it is not necessary to repeat the position determinations. +The spectrum positions in one of the images, generally the one with +the strong calibration spectra, are determined first. This image is +then used as a "template" to provide the initial position estimates for +the other images. If the positions are identical no further work is needed, +otherwise, the positions can be refined to correct for small changes in the +positions and shapes of the spectra. +.PP +The task \fBnewextraction\fR creates new databases. If a template image +is specified then a copy is made of the template image database. This means +that the number of spectra and the sample image lines remain the same. +If the spectrum positions are slightly different from the template image +then the task \fBfitgauss5\fR is used to determine the new positions. +.PP +The spectrum positons and possibly any model parameters are interpolated +from the sample lines to the remaining image lines by fitting a function +to values at the sample lines. In addition, the function fits may +leave out poorly determined points and also smooth the values at the +sample lines. The task \fBfitfunction\fR fits selected functions of +specified order to the selected spectra and sample image lines. +.PP +The extraction of the spectra from multi-spectra images is performed by +the task \fBmsextract\fR. The task extracts either integrated or strip +spectra, either data or model values, with or without blending corrections, +and with or without replacing bad pixels by model values. +The user specifies the limits of the extraction +strip as well as the spectra and image lines to be extracted. +.PP +For the simplest type of data extractions (basically strip extraction) +no modeling is required. Other types of extractions, such as model +extractions and/or with cleaning and blending corrections require some +degree of modeling. There are two models which may be used; +"smooth" and "gauss5". These models are described in section 5. +The model parameters for model gauss5 must be set by \fBfitgauss5\fR +before \fBmsextract\fR is used. Additional models may added for +extraction as well as for the spectrum position determinations. +.PP +The model based features of \fBmsextract\fR -- model extractions +and cleaning -- are available in the related task \fBnewimage\fR. +This task creates new images which consist of either model spectra +or cleaned data spectra. +.PP +The models in the \fBmultispec\fR package assume that the profiles +go to zero; i.e. there is no background light. Background light +may be removed using \fBbackground\fR. In the future a task will +be provided create a mask defining the locations of the spectra from +the database which can be used with general surface fitting tasks +to create a background surface to be subtracted from the image. +.PP +The final step in using the \fBmultispec\fR package is to convert the +extraction output to the desired format. This may include graphs, +card image formats, and files for the \fBonedspec\fR and \fBlongslit\fR +packages. Currently, the available formats are images and IIDS +card images. +.NH +The Tasks of the Multispec Package +.PP +Use of the \fBmultispec\fR package begins with \fBnewextraction\fR and +ends, usually, with \fBmsextract\fR. In between there are tasks which +update, refine or change the database and tasks which provide diagnositic +information. The informational tasks can be combined with tasks from +other packages to produce tabular or graphical output. The task +\fBmsplot\fR is an example. In this section a brief description of +each task is given. Further information about the tasks, including usage, +is available in the manual pages. +.SH +findpeaks +.IP +Selected sample image lines are examined to determine the number and +column positions of data peaks in the line. An average of a number of image +lines surrounding the sample lines is formed in which the local maxima +are located. Various criteria are applied to cull the list of local +maxima to the desired peaks. These criteria include a peak threshold, +a maximum peak-to-peak contrast, a minimum peak separation, and a +maximum number of peaks. This task is used to determine crude, initial +estimates for the spectrum positions. It could be used alone for +simple extractions. +.SH +fitfunction +.IP +This task has two roles. It's primary role is to define the +interpolation/extrapolation function for the spectra +positions between the sample lines. The fitting function can be +either purely interpolative or may also provide smoothing of the +parameters from the sample lines. The second role is to provide +smoothing of the model parameters along the dispersion and the +ability to replace bad values by the function fit to the remaining +parameters. In this second role the user may iterate between +smoothing and model fittng. The functions are always defined between +the first and last image lines. +.SH +fitgauss5 +.IP +The model profiles gauss5, described in section 5, are fit to the +selected spectra and sample lines. The parameters to be determined +and the fitting algorithm may also be selected. +The model parameters are recorded in the database. +The model may be tracked from a starting line to other sample image +lines or each sample line may be fitted independently. +This task is used to accurately determine the spectrum positions +and provide an extraction model for heavily blended spectra. +.SH +modellist +.IP +For the selected sample image lines and image columns data +and model values are listed. This task is used to check how well +the model fitting tasks (currently just \fBfitgauss5\fR) have fit +the sample image line. The task \fBmsplot\fR is used to produce +graphical output. +.SH +msextract +.IP +This task does the actual extraction of spectra. It requires that +the spectrum positions are defined by fitting functions in the +database. If model gauss5 is to be used then the database must +also contain the model parameters for the sample image lines. It +extracts integrated or strip spectra, using data or model values, +with or without blending corrections, and with or without cleaning +of bad pixels. +.SH +mslist +.IP +Of the diagnositic or informational tasks \fBmslist\fR is the most +general. The user selects the type of information from the database +which is desired and it is then printed. The types of information +include the database header, the database comments, the spectra +positions and model parameter values for the sample lines, and the +interpolation/smoothing function values for any desired set of +image lines. +.SH +msplot +.IP +This task extracts data and models values and plots them superposed. +This task is used as a diagnositic tool to inspect how well model fitting +represents the image spectra. +.SH +msset +.IP +This task is a general tool for modifying or setting some of the quantities +in the database. The quantity to be changed or set is +selected by a keyword and the values are input in two ways; +with a list structured parameter (such a file containing the list of +values or the standard input) or as a parameter value. This task +is the way a user may enter comments in the database or manually +set the number and positions of the spectra. It is also used to +set the initial values for the gauss5 model parameters s0, s1, and s2 +prior to using \fBfitgauss5\fR. +.SH +newextraction +.IP +This task has three important roles. First it creates the database +associated with the multi-spectra image. Second, it defines the sample +image lines to be used. The user can specify as many or as few sample lines +as desired. It should be kept in mind that the more sample lines used +the larger the database becomes and the longer the processing time when +modeling the spectra. Finally, \fBnewextraction\fR allows +a database from another image (called a template image) to initialize the +database for the new multi-spectra image. The template image is generally +a calibration image with strong, well-defined spectra. +Initializing a database with a template image saves time, reduces problems +with bad pixels, and is more accurate when an image with weak spectra is +to be extracted. +.SH +newimage +.IP +This task is similar to \fBmsextract\fR; it uses the same algorithms +and parameters. It differs in the type output. +Rather than producing extracted integrated or strip spectra this task +produces new image lines. It is particularly useful for extracting +model images to be compared against the original image or to +produce images which have been cleaned. +.NH +The Multispec Database +.PP +The tasks in the \fBmultispec\fR package create and manipulate a database. +The database contains a description of the multi-spectra image which +is modified, refined, examined, or otherwise used by the tasks in the package. +In the current version the database is a separate binary file with a filename +formed by appending ".db" to the image name described by the database. +.PP +The database contains four basic types of data; general information, +comments and history, position parameters, and model parameters. +The data in the database is examined with the task \fBmslist\fR. +The general information section, called the database header, contains the +the name of the image, the size of the image, and the number of spectra in +the image. Once the number of spectra in the image has +been entered in the database it is an error to attempt to change this +number. The database must be deleted and a new database created in order +to change the number of spectra. +.PP +The comment and history section of the database contains text +strings. Each task which modifies the contents of the database places +a dated history line in this section. The user may also add comments +with \fBmsset\fR. Currently this information is not passed on to +the extraction output. +.PP +There are three types of position information in the database. The +first is a set of sample image lines. The sample lines are set when +the database is created by \fBnewextraction\fR. The sample lines select +which image lines from the multi-spectra image are to be examined and used +during the extraction. Information from these sample lines, and only +these sample lines, is entered in the database. The sample lines +may be listed with \fBmslist\fR. +.PP +The second type of position information is the positions of the +spectra (centers) at each sample line. These positions are initially +set by either \fBfindpeaks\fR or, manually, by \fBmsset\fR. The +position information is refined by fitting model profiles. +.PP +The third type of position information is a function fit to the +positions from all the sample lines for each spectrum. +These function fits are produced by \fBfitfunction\fR. +The functions define the positions of the spectra at all the image +lines. The spectra positions at the sample lines or the function +evaluation for any image line may be listed with \fBmslist\fR. +.PP +The finally type of basic data contained in the database are +model parameter values. A model need not be used in the extraction +but if one is used then the parameters determining the model profiles +are recorded in the database. The specific parameters depend on the +model. Currently the only model is \fIgauss5\fR. The model and its +parameters are described in section 5. +.PP +As with the spectra positions the parameters are stored in the database +in two forms; as values for each spectrum at each sample image line +and as function fits to the values at the sample lines which interpolate +them to any image line. The sample line values are +set by the model fitting tasks and the function fits are set by +\fBfitfunction\fR. The parameter values at the sample lines or the +function evaluations for any image lines may be listed with \fBmslist\fR. +.NH +Multispec Spectrum Profile Models +.PP +The spectra profiles in the image are modeled for many reasons: +To provide accurate, subpixel position determinations, to extract model +spectra or model images, to detect and replace bad pixels, and +to estimate and correct for blending between the spectra. +There are currently two models used in the \fBmultispec\fR package, "gauss5" +and "smooth". +.NH 2 +Model Gauss5 +.PP +The gauss5 model profiles are Gaussian but with a scale which varies +smoothly between the center and the edge of the profile. There +are five parameters: + +.RS +.IP x0 +The column position in the image line of the center of the profile. +.IP i0 +The intensity scale of the profile. It corresponds to the intensity +of the center of the profile. +.IP s0 +The zeroth order, constant, term in the Gaussian scale. +.IP s1 +The even first order term in the Gaussian scale. +.IP s2 +The odd first order term in the Gaussian scale. +.RE + +.PP +The mathematical form of the the model is shown in equation (1): +.EQ (1) +roman profile (x)~=~i0 exp~left { -s( DELTA x )~DELTA x sup 2 right } +.EN +where +.EQ +DELTA x ~=~x~-~x0~, +.EN +.EQ +s( DELTA x)~=~s0~+~s1~|y| +~s2~y~, +.EN +and +.EQ +y~=~ DELTA x / ( DELTA x sup 2 + alpha ) sup half ~. +.EN +The profile is defined within the user specified limits \fIlower\fR and +\fIupper\fR measured relative to the the profile center and +$alpha~=~(upper-lower)/4$. The quantity $y$ lies in the range +-1 to 1 over the interval in which the profile is defined. The odd +and even terms, s1 and s2, allow for symmetric and antisymmetric profile +changes relative to a simple Gaussian profile. +.PP +The task \fBfitgauss5\fR fits the gauss5 model to the spectrum profiles in +the sample image lines to determine one or more of the model parameters for +each spectrum. The parameter values are stored in the database for the image. +In \fBmsextract\fR the model profiles for each +image line are obtained by interpolating the profile shapes from the sample +lines (with the model parameters in the database determined by +\fBfitgauss5\fR) and then fitting only the intensity scale "i0". +There are a number of technical details associated with the model fitting +in each of these tasks which are discussed in the manual pages. +.PP +The gauss5 model is used to accurately determine the positions of the +spectrum centers at the sample image lines. Fitting simultaneously +for the model parameters allows the spectra to be blended. +This is the chief advantage of this model. +This model is also used during extraction to correct for blending of +the spectra and to detect and replace bad pixels. +.NH 2 +Model Smooth +.PP +The spectrum profiles from the lines immediately preceeding +the image line in which the spectrum profile is to be fit are shifted +to a common center and averaged to form the model profile. +An intensity scale factor is then determined which best fits the model +profile to the image profile. This is done for each spectrum in the +image. The scale factors are determined by least squares with +possible bad pixel rejection. Rejected pixels are eliminated +when the image line is later used in forming new average model profiles. +.PP +The advantages of this model are that the image spectrum profiles may +have any shape and the least squares fitting with bad pixel rejection +is fast and rigorous. By passing through the image lines sequentially +the image lines need be accessed only once and the profile averages +can be quickly updated for the next image line. +.PP +The disadvantages of this model are that the spectrum profiles cannot +be blended and the model does not measure profile positions. +This means that the spectrum profile positions must be +known. This model is suitable for model extractions and cleaning of +bad pixels in unblended multi-spectra images. It is available in +the task \fBmsextract\fR. +.bp +.SH +Glossary +.LP +\fBmultispec\fR +.IP +Acronym for Multi-Spectra Extraction as in \fBmultispec\fR Package. +.LP +integrated spectra +.IP +The spectra are extracted by integrating the pixel values across the spectrum +to produce a single aperture luminosity value. +.LP +sample image line +.IP +The spectra positions and model profile shapes are determined at a set +of image lines selected when the database is created. +.LP +strip spectra +.IP +The spectra are extracted as a strip of fixed with the spectra shifted by +image interpolation to a common center. diff --git a/noao/twodspec/multispec/doc/newextract.hlp b/noao/twodspec/multispec/doc/newextract.hlp new file mode 100644 index 00000000..37123f28 --- /dev/null +++ b/noao/twodspec/multispec/doc/newextract.hlp @@ -0,0 +1,61 @@ +.help newextraction Jul84 noao.twodspec.multispec +.ih +NAME +newextraction -- Initialize a new MULTISPEC extraction +.ih +USAGE +newextraction image template +.ih +PARAMETERS +.ls image +Image to be extracted. +.le +.ls template +The previously created database for the template image is used to initialize +the new database. If the null string is given then the database is not +initialized. +.le +.ls sample_lines = "10x50" +Sample image lines in which the spectra positions are to be determined and, +optionally, modeled. This parameter is not used if a template image is given. +.le +.ih +DESCRIPTION +To extract the spectra from a multi-spectra image a database must be created +and associated with the image. This task creates the database with a name +formed by adding the extension '.db' and initializes some of the database +entries. + +The sample lines are used to track the spectra positions and, if an analytic +profile model is to be fit to the spectra, to map profile shape changes. +The image lines only need be sampled enough to track \fInon-linear\fR position +distortions and significant profile shape changes since interpolation +is used between the sample lines. Though specifying just one sample +line is allowed using at least two sample lines is recommended to allow for +any slope in the position of the spectra. Specifying all the image lines +will greatly increase the processing time and is never justified. + +Using a previous database to initialize the new database is useful if the +new image is only slightly different in the positions and profiles of the +spectra. In some cases extraction may proceed immediately without any +further position determination and modeling. Further modeling +and spectra position determinations will refine the previously determined +parameters with an increase in execution time. Using a template image is +particularly important if the first image extracted has strong spectra +and subsequent images have much weaker spectra since the automatic spectra +position location and profile modeling may yield poor results for very weak +spectra. +.ih +EXAMPLES +To initialize a MULTISPEC database for extracting the spectra in +the image \fIimage1\fR: + + cl> newextraction image1 "" + +To create a new MULTISPEC database for extracting the spectra in +the image \fIimage2\fR using \fIimage1\fR as a template image: + +.nf + cl> newextraction image2 image1 +.fi +.endhelp diff --git a/noao/twodspec/multispec/doc/newimage.hlp b/noao/twodspec/multispec/doc/newimage.hlp new file mode 100644 index 00000000..1ef7fbe0 --- /dev/null +++ b/noao/twodspec/multispec/doc/newimage.hlp @@ -0,0 +1,130 @@ +.help newimage Jul84 noao.twodspec.multispec +.ih +NAME +newimage -- Create a new multi-spectra image +.ih +USAGE +newimage image output +.ih +PARAMETERS +.ls image +Image to be used to create the new image. +.le +.ls output +Filename for the new multi-spectra image. +.le +.ls lower = -10 +Lower limit for model profiles. It is measured in pixels from the +spectra centers defined by the position functions in the database. +.le +.ls upper = -10 +Upper limit for model profiles. It is measured in pixels from the +spectra centers defined by the position functions in the database. +.le +.ls lines = "*" +Image lines of the multi-spectra image to be in the new multi-spectra image. +.le +.ls ex_model = no +Create a model image? +.le +.ls clean = yes +Replace bad pixels with model values? The following parameters are used: +.ls nreplace = 1000. +Maximum number of pixels to be replaced per image line when cleaning with +model \fIgauss5\fR or maximum number of pixels to be replaced per spectrum when +cleaning with model \fIsmooth\fR. +.le +.ls sigma_cut = 4. +The cleaning threshold in terms of the predicted pixel sigma. +.le +.ls niterate = 1 +Maximum number of cleaning iterations per line when cleaning with model +\fIgauss5\fR. +.le +.le +.ls model = "smooth" +Choice of \fIgauss5\fR or \fIsmooth\fR. Minimum match abbreviation is +allowed. This parameter is required only if \fIex_model\fR = yes +or \fIclean\fR = yes. +.le +.ls fit_type = 2 +Model fitting algorithm for model \fIgauss5\fR. +.le +.ls naverage = 20 +Number of lines to be averaged in model \fIsmooth\fR. +.le +.ls interpolator = "spline3" +Type of image interpolation function to be used. +The choices are "nearest", "linear", "poly3", "poly5", and "spline3". +Minimum match abbreviation is allowed. +.le +.ls verbose = no +Print verbose output? +.le +.ih +DESCRIPTION +A new multi-spectra image is created using the description of the +multi-spectra image in the MULTISPEC database associated with \fIimage\fR. +The user selects the image \fIlines\fR from the original image to be in +the new image. The options allow the creation of model images or images in +which the bad or deviant pixels are replaced by model profile values. + +If \fIex_model\fR = yes or \fIclean\fR = yes model +spectra are fit to the spectra in the image. There are two models: +a five parameter Gaussian profile called \fIgauss5\fR and profiles obtained +by averaging \fInaverage\fR image lines surrounding the image line being +modeled called \fIsmooth\fR. The model is selected with the parameter +\fImodel\fR. + +When \fIex_model\fR = yes an image containing model spectra is produced. + +When \fIclean\fR = yes pixels with large residuals from the model are +detected and removed from the model fit. The selected model is +fit to the pixels which are not in the bad pixel list (not yet implemented) +and which have not been removed from the model fit. The sigma of the fit +is computed. Deviant pixels are detected by comparing them to the model +to determine if they differ by more than \fIsigma_cut\fR times the sigma. +The model fit is iterated, removing deviant pixels at each iteration, until +no more pixels are found deviant or \fInreplace\fR pixels have been found. +The pixels removed or in the bad pixel list are then replaced with +model values. (To clean and extract the spectra with this algorithm see +\fBmsextract\fR.) + +There are some technical differences in the model fitting and cleaning +algorithms for the two models. In model \fIsmooth\fR +the fit for the profile scale factors is done independently for each spectrum +and automatically corrected when a bad pixel is detected. This fitting process +is fast and rigorous. The parameter \fInreplace\fR in this model refers to +the maximum number of pixels replaced \fIper spectrum\fR. + +In model \fIgauss5\fR, however, the profile scale factors are fit +to the entire image line (hence its ability to fit blended spectra). +There are two fitting algorithms; a rigorous simultaneous fit +and an approximate method. The simultaneous fit is selected when +\fIfit_type\fR = 1. This step is relatively slow. The +alternative method of \fIfit_type\fR = 2 sets the scale factor for each +spectrum by taking the median scale, where scale = data / model profile, +for the three pixels nearest the center of the profile. The median +minimizes the chance of a large error due to a single bad pixel. This +scale may be greatly in error in the case of extreme blending but is also +quite fast; the extraction time is reduced by at least 40%. +The steps of profile fitting and deviant pixel detection are alternated +and the maximum number of iterations through these two steps is +set by \fIniterate\fR. The default of 1 means that the model fitting is not +repeated after detecting deviant pixels. + +The option \fIverbose\fR can be used to print the image lines being extracted +and any pixels replaced by the cleaning process. +.ih +EXAMPLES +To create a cleaned version of the image using model \fIsmooth\fR for cleaning: + + cl> newimage image newimage + +To create an model image using model \fIgauss5\fR: + + cl> newimage image newimage ex_model=yes model="gauss5" +.ih +SEE ALSO +msextract +.endhelp |