7.5 Fitting of images

There are three tasks and two verbs which estimate the position and intensity of a component on a two-dimensional image. The simplest and fastest methods are the verb IMCENTER and MAXFIT. The latter fits a two-dimensional parabola to the maximum within a few pixels of an image position, and gives the peak and its position. The tasks IMFIT and JMFIT are similar and fit an image subsection with up to four Gaussian components with error estimates. Tasks SAD and TVSAD attempt to automate the process of finding and fitting Gaussian components in an image. Additionally, in one dimension, the task SLFIT fits Gaussian components to slice data and the task XGAUS fits Gaussian components to each row of an image.

7.5.1 Centroid fits (IMCENTER) and parabolic fit to maximum (MAXFIT)

You may determine a centroid for a region in an image with the verb IMCENTER. Set the name parameters for the desired image and then define the region with adverbs BLC and TRC, perhaps using TVLOAD; TVWIN  C R. Set FLUX if you wish to limit the computations to pixel values greater than FLUX. Then

MAXFIT’s speed makes it useful for simple regions. Type:

The inputs should be self-explanatory. The IMSIZE parameter can be important in crowded fields. MAXFIT can be used conveniently by first displaying the image on the TV and then typing:

First the cursor will appear on the TV. Move it close to a maximum, press the left mouse button, and hit button A, B, C, or D. The fit will appear in your  window. Adverb values PIXXY, PIXVAL, COORDINA, and ERROR will be set appropriately. Adverb FSHIFT will also be set to suggest the shift needed to place the source directly on a pixel.

7.5.2 Two-dimensional Gaussian fitting (IMFIT, JMFIT)

A more sophisticated least-squares fit of an image is obtained with IMFIT, which fits an image with up to four Gaussian components and attempts to derive error estimates. A linear or curved, two-dimensional “baseline” may also be fitted. A sample set-up is as follows:

To improve accuracy, include as small an area as possible in the fit. In some cases, it is useful to hold some of the parameters constant, particularly when fitting a complex clump of emission with several components. The parameters can interact. Error estimates are given for each component based on the image rms in the image header or determined by robust means over the full input image plane. IMFIT will sometimes fail to converge in complicated regions. When this happens, you might try using the task JMFIT, which is very similar in function, but uses a different mathematical method to minimize the rms. Comparison of the results of IMFIT and JMFIT will sometimes be instructive. The tasks will correct the results for the effects of the primary beam and bandwidth smearing if you wish. These tasks return adverbs FMAX, FPOS, and FWIDTH with the answers to the fit, DOMAX, DOPOS, and DOWIDTH with the uncertainties in the fits, and FSHIFT with the shift needed to place the first component directly on a pixel. It is wise to treat the results of MAXFIT, IMFIT and JMFIT with caution. The estimates of the errors, in particular, are based on theory and on trials of deconvolution over a range of widths.

The verb MFITSET may be used to set many of the input adverbs for IMFIT and JMFIT from the image displayed on the TV. The verb sets the image name adverbs, then directs you to set the fitting window (like TVWINDOW), and then to point at the peak of Gaussian component 1 (sets GMAX(i) and GPOS(*,i)), then its major axis half-power point (sets GWIDTH(1,i) and GWIDTH(3,i)), then its minor axis half-power point (sets GWIDTH(2,i)). Push buttons A, B, or C to set each of these points or push button D to stop the process and set NGAUSS. Instructions will appear on your terminal.

Alternatively, use RUN INPFIT  C R (see §12.2.1) to obtain a procedure which will help to supply input parameters to IMFIT. This RUN file loads a procedure called INPFIT into AIPS. To invoke it, load the image which you want to fit onto the TV with TVALL and type INPFIT ( 3 )  C R to specify three components. The procedure will prompt you to set the desired sub-image window with the TV cursor (it uses verb TVWINDOW) and then to point the TV cursor at the peaks of each of the Gaussians, click the left mouse button when the cursor is correctly placed, and push button A, B, C, or D. The inputs GMAX, GPOS, BLC, and TRC are set in this way.

7.5.3 Source recognition and fitting (SAD, TVSAD)

The task SAD (§10.4.4) attempts to find all sources in a sub-image whose peaks are brighter than a given level. It searches the sub-image specified by BLC and TRC for all points above this level and merges such points in contiguous “islands.” For each island, initial estimates of the strength, size, and number of components are generated. Then the fitting algorithm used in JMFIT is called to determine the least square Gaussian fit. Solutions which fail to meet certain criteria can be retried as two components and, if they still fail, rejected. SAD is a task with many adverbs, a full description of which would be beyond the scope of this ookook. Enter EXPLAIN SAD  C R for a full description of this task and its parameters. The effects of bandwidth smearing and the primary beam may be corrected. SAD produces a Model-Fit extension file which may be converted to a stars file (§6.3.2) with MF2ST. The MF file may be printed with MFPRT in formats suitable for STARS and in formats which may be used, with task BOXES, to prepare Clean boxes for input to the imaging tasks. SAD can now directly attach a “stars” extension file to the input image containing the fit Gaussians. Plot tasks can then overlay them on images.

A version of SAD called TVSAD has been written along with an  Memo¹ was written describing its use in detail. Like SAD it finds source islands and makes an initial guess. At this point it displays an interpolated image of the area surrounding the island and marks it with the island boundary and the initial guesses (plotted as ellipses at the half-power point). You may choose from a TV menu to enhance the display, to change the island boundary, to change one of the initial guesses, to change the number of Gaussians guessed and set their values (like verb MFITSET), to have the fit attempted, or to reject this island and go on to the next. After the fit, a residual image is displayed with the island boundary and fit Gaussians marked. At this stage you may accept the fit, go back to the previous menu to try again, or reject this island and go on the next island. You may turn off the TV interaction and have the task run automatically. It will turn the TV interaction back on if it finds a fit that does not meet the acceptance criteria. Finally, TVSAD writes out the residual image, printer displays, and extension files exactly like SAD.

7.5.4 Gaussian fits to slices (SLFIT)

You can generate a one-dimensional slice (profile) through any plane (characterized by the first two coordinates) of an image file using the  task SLICE. The output file is appended to the image file as an SL extension file. Slices are computed along lines in the two-dimensional image joining any valid pair of points selected by BLC and TRC. A slice along the third axis, at a single pixel (not necessarily integers) in the first two axes, may also be prepared. Slices directly along any single axis may be saved without interpolation if desired. If the slice is directly along an FQID axis, SLICE will output a slice which is linearly sampled in frequency. Tasks ISPEC, BLSUM, and TVSPC may also save slice extension files. These are somewhat different in that they can represent a sum over an area in the sky as a function of frequency. Nonetheless, the slice file display and fitting software is of considerable use with these as well. The set of software dealing with slice file analysis and display can be obtained on your terminal by typing ABOUT ONED  C R. The list is also given in Chapter 13.

To generate a slice:

Use INDISK and GETNAME to select the input image. The beginning (BLC) and ending (TRC) points for the slice can be specified conveniently using the TV cursor if the image to be sliced is first displayed on the TV with TVLOD or TVALL. To set these points with the TV, type:

then set the TV cursor to the desired beginning point for the slice, press the left mouse button, and repeat for the ending point for the slice. Note that, for slices, BLC need not be below or to the left of TRC. Finally:

Slice files may be output as ASCII text files using the OUTTEXT adverb. They contain detailed pixel and coordinate information as well as the interpolated slice values. They may also be printed later by SLPRT including any one of the models fit to it. Slice files are archived in your disk catalog as SL extensions to the image file from which they were derived. Running SLICE again with new parameters does not overwrite the slice file, but makes another with a higher “version” number. To review and/or delete slice files, follow the instructions for EXTLIST and EXTDEST of plot files in §6.3 above, but use INEXT ’SL’  C R in place of INEXT ’PL’  C R.

When SLICE has terminated, the file may be plotted in the TV display on your workstation using:

The default scales will plot all slice points on a vertical scale from the slice minimum to the slice maximum. You can alter the part of the slice that is plotted and the vertical scale by specifying, for example:

Note: several slices may be put on one TV plot. Use TVASLICE  C R for the additional ones. Multiple colors may be achieved by using different graphics channels (GRCHAN).

Slice files may be converted into plot files by:

The resulting plot files may then be output by:

The task SLFIT fits Gaussian components to one-dimensional data in slice files. Assuming that the usual GETNAME step has been done, a typical session would go like:

This verb will prompt you to POSITION CURSOR AT CENTER & HEIGHT OF GAUSSIAN COMP 1. Move the cursor to the requested position and hit any button. Then you are asked to POSITION CURSOR AT HALFWIDTH OF GAUSSIAN COMP 1. Move the cursor to the half-intensity point of the component and click any button. Continue until all components have been entered. (Note: these operations are also available on the Tek device with verbs beginning with TK. We recommend the TV versions since cursor reading in X-Windows emulations of Tek devices appears to be unreliable.) Then type:

If everything looks ok, then:

When the task gets an answer, the solution will be displayed as  messages, recorded in the message file, and recorded in the slice file itself. To get a hard copy of the results:

and, to display the results in the TV window, enter:

The TV plot may be easier to examine if you use different GRCHAN for the different types of plot. To get a higher quality plot of the results, an example of which is shown in §6.3.2.1, type:

7.5.5 Spectral parameter fitting (TVSPC, XGAUS, ZEMAN, RMFIT, FARS)

Task TVSPC assists users in getting to know their data cubes. It uses a 2-dimensional image plane of some parameter that is useful in visualizing the full field. A moment-zero image is one obvious choice. The user selects a pixel on this image using the TV display and a spectrum at that coordinate is displayed. The task is highly interactive allowing a quick choice of an interesting pixel. The spectrum, or a portion of the spectrum, at that point may be fit by up to 4 Gaussians plus baseline and may be saved as a SLice file. The corresponding spectrum from a second transposed cube may also be displayed, allowing for example the I and V or the Q and U polarizations to be displayed simultaneously. Optionally, a third spectral cube may be displayed with the spectral plane selected interactively from the spectral plots. See the associated  Memo² . In earlier releases, use XPLOT with DOTV 1 to familiarize yourself with your image cube.

A task called TVIEW offers interactive TV displays of an image cube in all 3 transpositions plus a reference image or spectrum. The task is similar to TVSPC but is designed solely to help the user become better acquainted with the data cube. The task is described in an  Memo³ .

XGAUS is an interactive task which can fit up to eight Gaussians and a linear baseline to each row of an image. In 31DEC13, it was redesigned so that it begins by setting up a table of the peak image intensities and fit results, so that it may be re-started multiple times with the same table, and so that it provides multiple means to edit the results before it finally writes its results as a set of n - 1 dimensional image files. Although XGAUS was designed for use primarily on transposed spectral-line cubes (see §8.5.2), it has a wide variety of other applications. The interaction is optional and uses the TV window on your workstation. The data, initial guess, model fit (total and components), and the residual for each row may be plotted on the TV screen. If the number of Gaussians being fit is larger than one, you may choose for each row to enter a revised initial guess using the cursor in the TV window. This process is similar to that of TVSET described above (§7.5.4). If things go well, you may turn off the interactive model setting; the task will turn it back on if there is a bad fit at some pixel. After all pixels have been fit, the task enters an image and edit mode. Images are computed for all Gaussian parameters and their uncertainties. A menu is shown offering displays of the images and several methods to revisit the solutions for selected pixels. This includes swapping the solutions between components and re-doing a list of individual pixels or all pixels for which the solutions or uncertainties exceed specified ranges. When you are satisfied with the result, you may instruct the task to write out images of the Gaussian parameters and their uncertainties. The use of XGAUS as well as ZEMAN, RMFIT, and the plotting and modeling tasks (below) is discussed in detail in  Memo 118.⁴

Task ZEMAN was written to fit for Zeeman splitting. The basic function fit to the V Stokes spectrum is a constant time the total intensity plus another constant times the derivative of the total intensity with frequency. The task is similar to XGAUS in that it first sets up a table of image intensities and fit results, then fits all pixels above specified intensities, and then offers the image display and edit modes. Unlike previous Zeeman programs (in other packages), ZEMAN also offers the option to use the results of XGAUS on the total intensity cube to fit the V Stokes parameter with a constant times the total intensity plus a separate constant for the derivative of each Gaussian component with frequency. This allows for multiple magnetic fields at the same pixel, each corresponding to a separate Gaussian spectral component.

Tasks AGAUS and ZAMAN have appeared to perform the same functions but on absorption-line spectral cubes. The fit is done to the observations, but the parameter that is Gaussian is the optical depth. These tasks and the essential mathematics are described in  Memo 122.⁵

Task XG2PL plots the solution for individual spectra from XGAUS. The corresponding ZEMAN fit may also be plotted in a separate panel. Absorption line and optical depth models may also be plotted.

Polarization spectra are often fit for rotation measure. The old task RM reads a transposed cube of the polarization angle made with multiple COMBs, then FQUBE or MCUBE followed by TRANS. It fits the rotation measure and the intrinsic magnetic field direction using one of two possible methods to resolve lobe ambiguities.

The task FARS reads two similar cubes of Q and U polarization images and performs a “rotation-measure synthesis.” The output of FARS is an image cube of the real part and another image cube of the imaginary part of the Fourier transform on the λ² axis, with or without one-dimensional Cleaning. Slice files of the complex RM transfer function (“dirty beam”) are also written. Task AFARS reads these cubes to produce images of the maximum rotation measure and of the amplitude or phase at the maximum. Task RFARS applies the rotation measure image to the FARS input cubes to produce de-rotated cubes of Q and U. These may be SQASH’ed to make an image of the average un-rotated polarization. The RUN file/procedure DOFARS assists with the whole process, doing matrix transpositions, finding channel weights, running FARS, transposing the outputs, and cleaning up.

Rotation-measure synthesis remains an experimental procedure at present. To test the handling of spectra in Q and U, the task TARS was written. It reads a text file giving the spectrum of Q and U at some single coordinate and performs the RM synthesis using methods from FARS plus an additional experimental form of complex Clean. Up to 20 model components may be added to the Q/U spectral data at specified rotation measure, amplitude, and phase. Task QUXTR was written to extract a spectrum from and pair of Q and U polarization cubes for input to TARS. The text-file output of TARS may then be plotted with TARPL, which allows multiple RM spectra from the same TARS output file to plotted on top of each other. In 31DEC20, new task MTARS prepares model Q/U spectra with up to 9999 model components and FTARS performs an interactive fitting of the input to TARS using the TARS output as a guide. This is RMFIT (see below) for a single spectrum. TARPL can now also lot the input to TARS.

It has been found that rotation-measure synthesis is severely limited in its ability to separate multiple RM components. Therefore, another XGAUS-like task called RMFIT was written. It uses the FARS rotation-measure output cubes to provide initial guesses for fits to the Q and U polarization spectral cubes. Like XGAUS, it first builds a table of image intensities to hold the fit results and then displays the FARS spectrum at each pixel to take an initial guess with which to fit the Q and U spectra. They are displayed and, if the fit is accepted, the task goes on to the next pixel. RMFIT is restartable and has the same sort of image/edit stage in which images are displayed and the fits may be revisited by means similar to those in XGAUS. An option to fit a spectral index in Q and U is offered, but, in general, the input cubes should have been corrected by the total intensity spectral index and the channel-dependent primary beam before being used in RMFIT. An additional spectral index in Q and U may then indicate “thick” rotation measure screens rather than a true source spectral index. Other models of thickness may be fit instead. The task may be re-started repeatedly and output images written only when desired. Despite rather simple fitting methods, RMFIT has shown the ability to separate closely spaced components not separable by rotation measure synthesis. Task RM2PL makes plot files from the solutions for individual Q/U or Polarization/angle spectra.