; XGAUS ;--------------------------------------------------------------- ;! Fits 1-dimensional Gaussians to images: restartable ;# TASK ANALYSIS PLOT SPECTRAL ONED MODELING TV-APPL ;----------------------------------------------------------------------- ;; Copyright (C) 1995, 2005, 2008, 2013-2014 ;; Associated Universities, Inc. Washington DC, USA. ;; ;; This program is free software; you can redistribute it and/or ;; modify it under the terms of the GNU General Public License as ;; published by the Free Software Foundation; either version 2 of ;; the License, or (at your option) any later version. ;; ;; This program is distributed in the hope that it will be useful, ;; but WITHOUT ANY WARRANTY; without even the implied warranty of ;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ;; GNU General Public License for more details. ;; ;; You should have received a copy of the GNU General Public ;; License along with this program; if not, write to the Free ;; Software Foundation, Inc., 675 Massachusetts Ave, Cambridge, ;; MA 02139, USA. ;; ;; Correspondence concerning AIPS should be addressed as follows: ;; Internet email: aipsmail@nrao.edu. ;; Postal address: AIPS Project Office ;; National Radio Astronomy Observatory ;; 520 Edgemont Road ;; Charlottesville, VA 22903-2475 USA ;----------------------------------------------------------------------- XGAUS LLLLLLLLLLLLUUUUUUUUUUUU CCCCCCCCCCCCCCCCCCCCCCCCCCCCC XGAUS Fits 1-dimensional Gaussians to images INNAME Input image name (name) INCLASS Input image name (class) INSEQ 0.0 9999.0 Input image name (seq. #) INDISK 0.0 9.0 Input image disk unit # INVERS 0.0 XG table version number OUTNAME Output image name (name) OUTCLASS Output image name (class) OUTSEQ -1.0 9999.0 Output image name (seq. #) OUTDISK 0.0 9.0 Output image disk unit #. BLC Bottom left corner of input TRC Top right corner of input YINC 0.0 128.0 Do every YINCth row in the first pass ZINC 0.0 128.0 Do every ZINCth plane in the first pass FLUX 0.0 Flux cutoff (see HELP) BCHAN -1.0 512.0 # channels at left for baseline (-1 => no baseline) ECHAN -1.0 512.0 # channels at right for baseline (0 => none) DOOUTPUT -1.0 3.0 Write residual image: 1,3 write parameter img: 2,3 DOTV -1.0 2.0 Plot data on TV; 2 => use TV not terminal for questions DORESID -1.0 1.0 Plot residuals on TV if DOTV PIXRANGE Min,Max of image intensity Max <= Min => entire range LTYPE -410.0 410.0 Type of labeling: 1 border, 2 no ticks, 3 standard, 4 rel to center, 5 rel to subim cen 6 map pixels PIXVAL 0.0 Display if peak < PIXVAL NITER 0.0 4000.0 Max # fit iterations 0 -> 100. NGAUSS 1.0 8.0 Number of Gaussians (<= 8) Guess of model parameters RMSLIMIT Switch from automatic to interactive if rms > RMSLIMIT BADDISK Disk to avoid for scratch ---------------------------------------------------------------- XGAUS Task: Fits one-dimensional Gaussians to each row of an image. (Normally XGAUS will be used on images with frequency or velocity as the first axis, but it will proceed with others as well.) The task fits up to 8 Gaussians plus a linear baseline to each row. Optionally, it writes up to 4 + 8*NGAUSS n-1 dimensional images containing the fit parameters and their uncertainties and/or the residual image. With a careful choice of input parameters, the task is capable of running in a batch-like mode. However, particularly for NGAUSS > 1, an interactive mode (using the TV graphics planes) is recommended. Interactivity is requested by setting DOTV greater than zero. Note that, since the TV is not allowed in BATCH jobs, XGAUS turns off all interactive options and TV displays when it is run in batch. The TV will let color differentiate what is plotted: DOTV > 0 (plot data, allow initial guesses to be set AND see what was fit). The initial guess uses channels 1-BCHAN and Nchan to Nchan-ECHAN+1 to set the DC level and slope. It is plotted as X signs. The data are plotted with stepped connecting lines, the fit as a smooth connecting line, and the residual as a dashed line. The guess is in graphics channel 2 (usually green), the data in 1 (usually yellow), the fit in 4 (usually cyan), and the residual in 3 (usually pink). In 31DEC13, XGAUS has been transformed into the task it should have been. When run in interactive mode, no output files are created except for a large table file used to hold the results of the fitting (an XG file). The table is initialized with pixel positions and fluxes before any fitting is done. Then the fitting is done sequentially through the input cube in two passes. The first pass does every YINC'th row and ZINC'th plane, while the second pass does every pixel not already done. After that, in interactive mode, a menu-driven routine is brought up to allow you to review the fits including re-doing poor ones, rearranging the numbering of the components, and other editing of the results. Finally output files are written if requested or if the task has been run in a non-interactive mode. XGAUS may be re-started using the XG file either resuming the sequential fitting (if it was not completed) or dropping into the edit stage directly. You may create the XG file with expansive ranges for BLC(2) - TRC(2), BLC(3) - TRC(3) and then quit as soon as the first spectrum to be fit is displayed. You may then re-start using rather smaller ranges in Y and Z covering particular regions of similar line shapes. The output images will cover the entire range initially specified. When doing this, it is a good idea not to change BLC and TRC (1,4,5,6,7). Note too that the various sessions are allowed to use different values of NGAUS. The table routines always handle 8 Gaussians, but it is much easier to provide guesses for a small number of Gaussians whose values will be then stored in the low numbered Gaussians. In the edit stage, you could specify a large number of Gaussians and move some of these low numbered solutions into higher numbers. Details of the interactive operation and options are described in the EXPLAIN file. Adverbs: INNAME.....Input image name (name). Standard defaults. INCLASS....Input image name (class). Standard defaults. INSEQ......Input image name (seq. #). 0 => highest. INDISK.....Disk drive # of input image. 0 => any. INVERS.....Input XG table version. 0 => a new one. Set this adverb if you wish to re-start a fitting session. OUTNAME....Output image name (name). Standard defaults. OUTCLASS...Output image name (class). Standard defaults. Used only for the residual image. The OUTCLASSes for the fit parameters are SLOPE, CONST, AMPLn, CENTRn, and WIDTHn and for the uncertainties are DSLOPE, DCONST, DAMPLn, DCENTn, and DWIDTn. The "flux" of each component (ampl * width * 1.06) is written with OUTCLASSes of FLUXn and DFLUXn. OUTSEQ.....Output image name (seq. #). 0 => highest unique. OUTDISK....Disk drive # of output image. 0 => highest number with sufficient space. BLC........Bottom right corner in input image of desired subimage. Default is entire image. If INVERS is specified, this BLC - TRC region in axes 2 and 3 must match or be interior to the region specified when the XG table was created. Processing is restricted to the current BLC to TRC. TRC........Top right corner in input image of desired subimage. Default is entire image. YINC.......Do only every YINC'th row (beginning at BLC(2)) in the first pass. Do the other rows in the second pass. ZINC.......Do only every ZINC'th plane (beginning at BLC(3)) in the first pass. Do the other planes in the second pass. FLUX.......A flux cutoff in the same units as the input image (i.e. Jy/beam). If a row does not have three consecutive points above this level, no component is fit to the row. It is also used to limit the points which determine the initial guess. If you do a re-start with a new lower value of FLUX, XGAUS will go into a special sequential mode to fit all rows with values between the previous and new FLUX levels. BCHAN......Number of points at the left edge of the input row used to determine the initial guess for the linear baseline. -1 => do not fit a baseline. 0 => no channels at the left, i.e. the initial guess for the baseline slope will be zero. ECHAN......Number of points at the right edge of the input row used to determine the initial guess for baseline. 0 => no channels at the right, i.e. the initial guess for the baseline slope will be zero. DOOUTPUT...= 1 or 3 requests that the residual (data-model) image be written out as a cataloged image using OUTCLASS. = 2 or 3 requests that the model images be written out as cataloged images. Note that DOOUTPUT can be changed interactively during the imaging portion of the task and the value of the adverb at exit is all that matters. DOTV.......True (> 0) implies plot the data, the initial guess, and the fit model on the TV. It must be set > 0. if you wish the option to correct initial guesses with the TV cursor. If DOTV > 0 and NGAUSS > 1, the program pauses after plotting its initial guess and requests user input from a TV menu (DOTV=2) or the terminal. Click DO_FIT (hit Enter key on terminal) to proceed with the program's initial guess. Enter BAD to have the solution for the row blanked or QUIT to exit the task directly. Click RE-GUESS (E or e on the terminal) to have the program request your guidance as to a better initial guess. During the entering of a new initial guess with the TV cursor, you may request that a particular component be null (held at 0 amplitude) by indicating that its peak or half-width point is outside the plot area. DORESID....True (> 0) requests a plot of the final residual values in each row. If DOTV is true, the program will pause after each row and request input from a TV menu (DOTV=2) or the terminal. At that point, click GOOD (Enter key on the terminal) to keep the answer (for now). You may enter BAD to have the solution blanked. The program will warn you if any of the fit parameters are way out of normal limits. For NGAUSS > 1, one may also enter RETR to loop back and try again with a manually entered initial guess. PIXRANGE...Min,Max of Image intensity. 0,0 => min,max of EACH row (separately). LTYPE......Labelling type, see HELP LTYPE for details: 1 = border, 2 = no ticks, 3 or 7 = standard, 4 or 8 = relative to ref. pixel, 5 or 9 = relative to subimage (BLC, TRC) center, 6 or 10 = pixels. 7-10 all labels other than tick numbers and axis type are omitted. Less than 0 is the same except that the plot file version number and create time are omitted. Add n * 100 to alter the metric scaling. PIXVAL.....Plot row only if the peak value < PIXVAL. Limits the plots to those of uncertain intensities. 0 -> 10^9. NITER......Maximum function evaluations during the fit of each row. (< 10 -> 100 for NGAUSS > 1, < 100 -> 150 for 1 Gaussian) NGAUSS.....Number of Gaussians to fit (between 1 and 8). RMSLIMIT...When the fitting has been told to continue without the TV, switch interactive mode back on whenever the Q and/or U rms exceeds RMSLIMIT. 0 -> 1000000. BADDISK....Disk drives to avoid for scratch files. ---------------------------------------------------------------- IMAGR: One-dimensional Gaussian fitting of image cubes Documenter: E. W. Greisen NRAO Related Programs: CUBIT, GAL, IMFIT, JMFIT, SAD, SLFIT Beginning with the 31DEC13 version of AIPS, XGAUS has been overhauled into the task it should have been. The interactive Tektronix option was removed since Tek emulation in xterm windows has proved unreliable and since full TV menus were added to the later stages of the task. At present, XGAUS assumes that there are at most 3 axes with more than one pixel in the input image. XGAUS begins by making an XG table containing one row for every pixel in the output plane from BLC(2,3) to TRC(2,3). Each row of the table is initialized with the peak value in each row of the input image, where the "peak" is the largest value averaged over 3 adjacent pixels. Then XGAUS begins fitting those rows having a peak value above FLUX, doing every YINC pixel in the Y axis and ZINC pixel in the Z axis. Initial guesses for each fit are found from the row data when NGAUSS = 1 or from a previous fit for NGAUSS > 1. This is one of the reasons the interactive mode is recommended strongly for NGAUSS > 1. After each input image row is fit, the results are added to the XG table immediately. This allows the interactive user to quit at any time and then restart the process later (set INVERS to point at the pre-existing XG table which you intend to modify). After the YINC by ZINC pass through the data, a second pass through the data is done on every pixel. Those pixels which have already been fit are skipped, but the results from them contribute to the next initial guess. Finally, when all pixels have been fit, XGAUS enters an interactive routine designed to improve the results before any output images are written. However, before describing that function, we need to provide details of the interactive fitting of each row of the input image. This process is followed whenever a row is fit in the first YINC/ZINC function, in the second every pixel function, and in the third result editing function. For each row, the first step in the fitting process is to determine an initial guess non-interactively. If XGAUS is currently in interactive mode, the next steps are: 1. If DOTV > 0 (STRONGLY recommended), the input data are plotted on the TV and the initial guess is added to the plot of the data. 2. If DOTV > 0 and either NGAUSS > 1 or a fitting has been forced, XGAUS asks in the AIPS window for instructions using a TV menu (DOTV=2) or the terminal input window. If your answer begins with E or e (RE-GUESS on the TV), XGAUS will ask you to enter a new initial guess (step 3). If the answer begins with B or b (BAD on the TV), XGAUS will mark this pixel as bad and, if the answer begins Q or q (QUIT on the TV), XGAUS will simply exit, allowing you to restart it at a later time. If the answer is anything else, e.g. a simple carriage return ("Enter") or GOOD on the TV, XGAUS goes on to step 4 below. 3. To enter a new initial guess, you will be prompted to position the TV cursor at the center and height of each Gaussian (note that the X and Y positions matter) and at the half-width of each Gaussian (note that only the X position matters). To mark each point, hit any button (A, B, C, or D). XGAUS then returns to step 1 above to plot the new guess and ask again. 4. Once an acceptable initial guess has been found, XGAUS proceeds to call a non-linear least squares fitting routine to determine the Gaussian parameters that appear to fit the data best. the answers are then checked to see if they are "reasonable" - negative components, components centered outside the input data, and components < 1.5 cells wide (FWHM) are thought to be unreasonable. 5. If the result is unreasonable and the current mode is not interactive, but it started as interactive, the interactive mode is turned back on (with a message and a display of the answers in the message terminal window) and XGAUS returns to step 1. If the task did not start as interactive and the result is unreasonable, the solution for that input row is marked bad and the task goes on to the next row. 6. If the current mode is interactive, the residuals are added to the plot (if DORESID > 0) and the final model is also added to the plot (if DOMODL > 0). In the input terminal, you are then told if the answers are "unreasonable" and are shown the answers, reasonable or not. You are then offered a variety of choices. If your response begins with: a) B or b - the solution is marked as bad and XGAUS goes on to the next input row (BAD on the TV) b) Q or q - XGAUS exits cleanly, leaving the XG table to be worked on again at a later time (QUIT on the TV). c) R or r - XGAUS returns to step 1 to try again with the current number of Gaussians (<= NGAUSS) (RE-GUESS on the TV). c) E or e - XGAUS returns to step 1 to try again with the current number of Gaussians (<= NGAUSS) (RE-GUESS on the TV). e) 1, 2, 3, 4, 5, 6, 7, 8 - XGAUS returns to step 1 to try again with the specified number of Gaussians (note that only numbers <= NGAUSS are respected) (1, 2, 3, 4, 5, 6, 7, and 8 on the TV as needed) f) H or h - XGAUS prompts you to enter the Gaussian parameters for each component (HAND on the TV). You enter the peak in input image units, the center in pixels wrt the reference pixel, and width in pixels. All 3 numbers must appear in one line. You may also put in flags to control whether your hand-entered numbers may be changed if the spectrum is re-fit. Values > 0 mean that the parameters will be fit, and all omited flags are taken as 1. It will prompt for all NGAUSS components, but will change the prompt if the current number of Gaussians is < NGAUSS. If you enter something besides 0 in that case, the current number of Gaussians is increased appropriately. XGAUS then returns to the start of this step 6 to allow you to see if you made a good guess. g) D or d - After a HAND operation, a DO FIT option is available. It loops back using the hand-entered parameters as the initial guess to the fitting routine. h) T or t - XGAUS turns off the interactive mode and does solutions in a batch-like fashion until an unreasonable solution is found or it finishes the current function. (TVOFF on the TV) i) Anything else - XGAUS saves the answer (with uncertainties) and goes on to the next input row (GOOD on the TV) Note, on retries (R, E, 1, 2, 3, and 4), step 2 will skip the questions and go on to step 3 directly (as if you had answered E). Recommended adverb values are DOTV = 2 and probably DORESID = 1. This allows full information and interactivity and uses different graphics planes for the different plots. DOTV = 2 will save moving back and forth between the TV and the terminal all the time. After all pixels that are strong enough have been given a solution, XGAUS enters a menu-driven function. The menu has in the left column: ------------------- | EXIT | Exit XGAUS, writing output images if DOOUTPUT is now > 0. | SET MIN S/N | Set minimum amplitude S/N(s) for okay solutions | SET MAX RES | Set maximum residual for okay solutions | SET PEAK RANGE | Set peak value range(s) for okay solutions | SET OFFX RANGE | Set offset range(s) for okay solutions | SET WIDTH RANGE | Set width range(s) for okay solutions | SET MAX ERR WID | Set maximum error(s) in width for okay solutions | REDO ALL | Re-do all solutions which are not okay | FLAG ALL | Mark bad all solutions which are not okay | OFF ZOOM | Turn of TV zoom | OFF TRANSFER | Turn off black & white and color TV enhancements | SET DOOUTPUT | Increment DOOUTPUT in loop 0-3 - with 1 and 3 causing residual images and 2 and 3 causing parameter images to be written on EXIT | ADD TO LIST | Type in output pixel coordinates to add to list | SHOW LIST | Display coordinates in list | REDO LIST | Re-do solutions for all pixels in list | FLAG LIST | Flag solutions for all pixels in list | SWAP LIST 1-2 | Swap solutions for components 1 and 2 for all pixels in list | SWAP LIST 1-3 | Swap solutions for components 1 and 3 for all pixels in list | SWAP LIST 2-3 | Swap solutions for components 2 and 3 for all pixels in list | SWAP LIST 1-4 | Swap solutions for components 1 and 4 for all pixels in list | SWAP LIST 2-4 | Swap solutions for components 2 and 4 for all pixels in list | SWAP LIST 3-4 | Swap solutions for components 3 and 4 for all pixels in list ------------------- Only appropriate SWAP LIST options are actually displayed. Three editing concepts are present in this menu. The first is to establish what parameter values (S/N in peak, peak residual, peak values, center values, width values, and error in widths) for each component constitute an "okay" solution. The selected limits are displayed above the menu. Then you can choose to try to REDO ALL solutions which are not okay or simply FLAG ALL, marking them as bad. The second concept is to create a list of up to 1000 pixels which are thought to have some problem. You can enter these pixels by hand (or remove pixels from the list by hand) with ADD TO LIST. A faster way to add to the list is to show the image of some parameter and then select interesting pixels while doing CURVALUE (see below). The list of pixels can be re-done (REDO LIST) or marked as bad (FLAG LIST). The list is cleared by these operations so you can make a new list as needed. The third editing method available here is to take the pixel list and swap the solutions at those pixels between component N with component M. The right hand column of options include: ------------------ | SHOW IMAGE A1 | Enter image interaction with peak value of component 1 | SHOW IMAGE C1 | Enter image interaction with center pixel of component 1 | SHOW IMAGE W1 | Enter image interaction with width of component 1 | SHOW IMAGE F1 | Enter image interaction with "flux" of component 1 | SHOW IMAGE EA1 | Enter image interaction with uncertainty in peak value of component 1 | SHOW IMAGE EC1 | Enter image interaction with uncertainty in center pixel of component 1 | SHOW IMAGE EW1 | Enter image interaction with uncertainty in width of component 1 | SHOW IMAGE EF1 | Enter image interaction with uncertainty in "flux" of component 1 ------------------ For NGAUSS > 1, appropriate additional choices are offered. For NGAUSS > 4, additional columns are needed in the menu. When you select one of the above options, the above menus are turned off, the image plane maintained in memory is displayed on the TV screen, and a new menu is displayed. If offers: --------------- | RETURN | Return to the above menus, image stays displayed | LOAD AS SQ | Re-load image with square root transfer function | LOAD AS LG | Re-load image with log transfer function | LOAD AS L2 | Re-load image with extreme log transfer function | LOAD AS LN | Re-load image with linear transfer function | OFF TRANSF | Turn off enhancement done with TVTRANSF | OFF COLOR | Turn off color enhancements | TVTRANSF | Black & white image enhancement | TVPSEUDO | Color enhancement of numerous sorts | TVPHLAME | Color enhancement of flame type, multiple colors | TVZOOM | Interactive zooming and centering of image | CURVALUE | Display value under cursor, mark pixels for list | SWAP 1-2 | Swap solutions for components 1 and 2 interactively | SWAP 1-3 | Swap solutions for components 1 and 3 interactively | SWAP 2-3 | Swap solutions for components 2 and 3 interactively | SWAP 1-4 | Swap solutions for components 1 and 4 interactively | SWAP 2-4 | Swap solutions for components 2 and 4 interactively | SWAP 3-4 | Swap solutions for components 3 and 4 interactively | NEXT WINDOW | Move to next window into large images --------------- For NGAUSS < 4, appropriate SWAP options are suppressed. Only one of the LOAD AS xx options is offered - SQ when the current function type is LN, LG when the current function type is SQ, L2 when the current function type is LG, and LN when the current type is L2. Selecting this option, reloads the image with the newly selection function type and changes the menu option accordingly. The options OFF TRANSF through TVZOOM are essentially the same as the corresponding verbs in AIPS. CURVALUE is like the verb of that name, but, when you press buttons A or B, the pixel under the cursor is added to the editing list. The SWAP n-m options invoke the interactive polygon setting operation used by the AIPS verb TVSTAT. This starts in the "set polygon n" mode where button A selects a vertex in the polygon, B sets the last vertex in the polygon and gets ready to set polygon n+1, C sets the last vertex and enters a vertex editing mode, and D sets the last vertex and exits the polygon setting. In the vertex editing mode, move the cursor to a vertex to be moved and press button A or B to reset its position. After moving it to the desired point, press button A or B to fix that point and restart the vertex editing mode. Press button C to fix that point and then start creating polygon n+1. Press button D to fix that point and exit the polygon setting. Once button D has been set, the selected areas of pixels have their solutions swapped. The images are updated and reloaded automatically. Very large images may not be able to fit on your TV. When the image exceeds the size of the TV, the top line will show the component number followed by a subimage number in parentheses and the NEXT WINDOW option will appear. The first sub-images displayed is called number 0 and shows the full image every n'th pixel in X and Y. Sub-image 1 begins at the lower left, moves right, then back to the left and up, and so forth until the top right is reached. Every pixel is displayed in these sub-images. Use the NEXT WINDOW option to step through the sub-images in a circular fashion.