`process_image`: processing an image

A standard analysis is performed using the process_image task. This task reads in the input image, calculates background rms and mean images, finds islands of emission, fits Gaussians to the islands, and groups the Gaussians into sources. Furthermore, the process_image task encompases a number of modules that allow decomposing an image into shapelets, calculating source spectral indices, deriving source polarization properties, and correcting for PSF variations across the image.

When process_image is executed, PyBDSF performs the following steps in order:

Reads in the image and collapses specific frequency channels with weights (see Multichannel options) and produces a ‘continuum’ image (the ch0 image) for all polarisations with which source detection is done.
Calculates basic statistics of the image and sensible values of the processing parameters. First, the number of beams per source is calculated (see Determining whether an image is confused for details), using a sensible estimate of box size and step size (which can be set using the rms_box parameter). Next, the thresholds are set. They can either be hard thresholded (by the user or set as 5-sigma for pixel threshold and 3-sigma for island boundaries by default) or can be calculated using the False Detection Rate (FDR) method using a user defined value for \(\alpha\) (the fdr_alpha parameter). If the user does not specify whether hard thresholding or FDR thresholding should be applied, one or the other is chosen internally based on the ratio of expected false pixels to true pixels.
Calculates the rms and mean images. The 3-sigma clipped rms and mean are calculated inside boxes of defined by the rms_box parameter. Optionally, these images can be calculated using adaptive scaling of this box, so that a smaller box (defined the the rms_box_bright parameter) is used near bright sources (where strong artifacts are more likely). Intermediate values are calculated using bicubic spline interpolation by default (the order of the spline interpolation can be set with the spline_rank parameter). Depending on the resulting statistics (see Determining whether an image is confused for details), we either adopt the rms image or a constant rms in the following analysis.
Identifies islands of contiguous emission. First all pixels greater than the pixel threshold are identified. Next, starting from each of these pixels, all contiguous pixels (defined by 8-connectivity, i.e., the surrounding eight pixels) higher than the island boundary threshold are identified as belonging to one island, accounting properly for overlaps of islands.
Fits multiple Gaussians to each island. The number of multiple Gaussians to be fit can be determined by three different methods (using the ini_gausfit parameter). With initial guesses corresponding to these peaks, Gaussians are simultaneously fit to the island using the Levenberg-Marqhardt algorithm. Sensible criteria for bad solutions are defined (see Flagging options). If multiple Gaussians are fit and one of them is a bad solution then the number of Gaussians is decreased by one and fit again, until all solutions in the island are good (or zero in number, in which case it’s flagged). After the final fit to the island, the deconvolved size is computed assuming the theoretical beam, and the statistics in the source area and in the island are computed and stored. Errors on each of the fitted parameters are computed using the formulae in Condon (1997) [1] and Condon et al. (1998) [2].
Groups nearby Gaussians within an island into sources. See Grouping of Gaussians into sources for details. Fluxes for the grouped Gaussians are summed to obtain the total flux of the source (the uncertainty is calculated by summing the Gaussian uncertainties in quadrature). The source position is set to be its centroid (the position of the maximum of the source is also calculated and output). The total source size is measured using moment analysis (see http://en.wikipedia.org/wiki/Image_moment for a nice overview of moment analysis). Errors on the source position and size are estimated using Condon (1997) and Condon et al. (1998), but additional uncertainties due to uncertainties in the constituent Gaussians may be estimated using a Monte Carlo technique.
Continues with further processing, if the user has specified that one or more of the following modules be used:
- Shapelet decomposition (see Shapelet decomposition module for details)
- Wavelet decomposition (see À trous wavelet decomposition module for details)
- Estimation of PSF variation (see PSF variation module for details)
- Calculation of polarization properties (see Polarization module for details)
- Calculation of spectral indices (see Spectral index module for details)

General reduction parameters

Type inp process_image to list the main reduction parameters:

PROCESS_IMAGE: Find and measure sources in an image.
================================================================================
filename ................. '': Input image file name
adaptive_rms_box ..... False : Use adaptive rms_box when determining rms and
                               mean maps
advanced_opts ........ False : Show advanced options
atrous_do ............ False : Decompose Gaussian residual image into multiple
                               scales
beam .................. None : FWHM of restoring beam. Specify as (maj, min, pos
                               ang E of N) in degrees. E.g., beam = (0.06, 0.02,
                               13.3). None => get from header
flagging_opts ........ False : Show options for Gaussian flagging
frequency ............. None : Frequency in Hz of input image. E.g., frequency =
                               74e6. None => get from header. For more than one
                               channel, use the frequency_sp parameter.
interactive .......... False : Use interactive mode
mean_map .......... 'default': Background mean map: 'default' => calc whether to
                               use or not, 'zero' => 0, 'const' => clipped mean,
                               'map' => use 2-D map.
multichan_opts ....... False : Show options for multi-channel images
output_opts .......... False : Show output options
polarisation_do ...... False : Find polarisation properties
psf_vary_do .......... False : Calculate PSF variation across image
rms_box ............... None : Box size, step size for rms/mean map calculation.
                               Specify as (box, step) in pixels. E.g., rms_box =
                               (40, 10) => box of 40x40 pixels, step of 10
                               pixels. None => calculate inside program
rms_map ............... None : Background rms map: True => use 2-D rms map;
                               False => use constant rms; None => calculate
                               inside program
shapelet_do .......... False : Decompose islands into shapelets
spectralindex_do ..... False : Calculate spectral indices (for multi-channel
                               image)
thresh ................ None : Type of thresholding: None => calculate inside
                               program, 'fdr' => use false detection rate
                               algorithm, 'hard' => use sigma clipping
thresh_isl ............. 3.0 : Threshold for the island boundary in number of
                               sigma above the mean. Determines extent of
                               island used for fitting
thresh_pix ............. 5.0 : Source detection threshold: threshold for the
                               island peak in number of sigma above the mean. If
                               false detection rate thresholding is used, this
                               value is ignored and thresh_pix is calculated
                               inside the program

Each of the parameters is described in detail below.

filename

This parameter is a string (no default) that sets the input image file name. The input image can be a FITS or CASA 2-, 3-, or 4-D cube.

adaptive_rms_box

This parameter is a Boolean (default is False). If True, an adaptive box is used when calculating the rms and mean maps. See Adaptive box options for details of the options.

advanced_opts

This parameter is a Boolean (default is False). If True, the advanced options are shown. See Advanced options for details of the advanced options.

atrous_do

This parameter is a Boolean (default is False). If True, wavelet decomposition will be performed. See À trous wavelet decomposition module for details of the options.

beam

This parameter is a tuple (default is None) that defines the FWHM of restoring beam. Specify as (maj, min, pos ang E of N) in degrees. E.g., beam = (0.06, 0.02, 13.3). For more than one channel, use the beam_spectrum parameter. If the beam is not given by the user, then it is looked for in the image header. If not found, then an error is raised. PyBDSF will not work without knowledge of the restoring beam.

flagging_opts

This parameter is a Boolean (default is False). If True, the Gaussian flagging options will be listed. See Flagging options for details of the options.

frequency

This parameter is a float (default is None) that defines the frequency in Hz of the input image. E.g., frequency = 74e6. For more than one channel, use the frequency_sp parameter. If the frequency is not given by the user, then it is looked for in the image header. If not found, then an error is raised. PyBDSF will not work without knowledge of the frequency.

interactive

This parameter is a Boolean (default is False). If True, interactive mode is used. In interactive mode, plots are displayed at various stages of the processing so that the user may check the progress of the fit.

First, plots of the rms and mean background images are displayed along with the islands found, before fitting of Gaussians takes place. The user should verify that the islands and maps are reasonable before preceding.

Next, if atrous_do = True, the fits to each wavelet scale are shown. The wavelet fitting may be truncated at the current scale if desired.

Lastly, the final results are shown.

mean_map

This parameter is a string (default is 'default') that determines how the background mean map is computed and how it is used further.

If 'const', then the value of the clipped mean of the entire image (set by the kappa_clip option) is used as the background mean map.

If 'zero', then a value of zero is used.

If 'map', then the 2-dimensional mean map is computed and used. The resulting mean map is largely determined by the value of the rms_box parameter (see the rms_box parameter for more information).

If 'default', then PyBDSF will attempt to determine automatically whether to use a 2-dimensional map or a constant one as follows. First, the image is assumed to be confused if bmpersrc_th < 25 or the ratio of the clipped mean to rms (clipped mean/clipped rms) is > 0.1, else the image is not confused. Next, the mean map is checked to see if its spatial variation is significant. If so, then a 2-D map is used and, if not, then the mean map is set to either 0.0 or a constant depending on whether the image is thought to be confused or not.

Generally, 'default' works well. However, if there is significant extended emission in the image, it is often necessary to force the use of a constant mean map using either 'const' or 'mean'.

multichan_opts

This parameter is a Boolean (default is False). If True, the multichannel options will be listed. See Multichannel options for details of the options.

output_opts

This parameter is a Boolean (default is False). If True, the output options will be listed. See Output options for details of the options.

polarisation_do

This parameter is a Boolean (default is False). If True, polarization properties will be calculated for the sources. See Polarization module for details of the options.

psf_vary_do

This parameter is a Boolean (default is False). If True, the spatial variation of the PSF will be estimated and its effects corrected. See PSF variation module for details of the options.

rms_box

This parameter is a tuple (default is None) of two integers and is probably the most important input parameter for PyBDSF. The first integer, boxsize, is the size of the 2-D sliding box for calculating the rms and mean over the entire image. The second, stepsize, is the number of pixels by which this box is moved for the next measurement. If None, then suitable values are calculated internally.

In general, it is best to choose a box size that corresponds to the typical scale of artifacts in the image, such as those that are common around bright sources. Too small of a box size will effectively raise the local rms near a source so much that a source may not be fit at all; too large a box size can result in underestimates of the rms due to oversmoothing. A step size of 1/3 to 1/4 of the box size usually works well.

Note

The spline_rank parameter also affects the rms and mean maps. If you find ringing artifacts in the rms or mean maps near bright sources, try adjusting this parameter.

rms_map

This parameter is a Boolean (default is None). If True, then the 2-D background rms image is computed and used. If False, then a constant value is assumed (use rms_value to force the rms to a specific value). If None, then the 2-D rms image is calculated, and if the variation is statistically significant then it is taken, else a constant value is assumed. The rms image used for each channel in computing the spectral index follows what was done for the channel-collapsed image.

Generally, the default value works well. However, if there is significant extended emission in the image, it is often necessary to force the use of a constant rms map by setting rms_map = False.

shapelet_do

This parameter is a Boolean (default is False). If True, shapelet decomposition of the islands will be performed. See Shapelet decomposition module for details of the options.

spectralindex_do

This parameter is a Boolean (default is False). If True, spectral indices will be calculated for the sources. See Spectral index module for details of the options.

thresh

This parameter is a string (default is None). If thresh = 'hard', then a hard threshold is assumed, given by thresh_pix. If thresh = 'fdr', then the False Detection Rate algorithm of Hopkins et al. (2002) is used to calculate the value of thresh_pix. If thresh = None, then the false detection probability is first calculated, and if the number of false source pixels is more than fdr_ratio times the estimated number of true source pixels, then the 'fdr' threshold option is chosen, else the 'hard' threshold option is chosen.

thresh_isl

This parameter is a float (default is 3.0) that determines the region to which fitting is done. A higher value will produce smaller islands, and hence smaller regions that are considered in the fits. A lower value will produce larger islands. Use the thresh_pix parameter to set the detection threshold for sources. Generally, thresh_isl should be lower than thresh_pix.

Only regions above the absolute threshold will be used. The absolute threshold is calculated as abs_thr = mean + thresh_isl * rms. Use the mean_map and rms_map parameters to control the way the mean and rms are determined.

thresh_pix

This parameter is a float (default is 5.0) that sets the source detection threshold in number of sigma above the mean. If false detection rate thresholding is used, this value is ignored and thresh_pix is calculated inside the program

This parameter sets the overall detection threshold for islands (i.e. thresh_pix = 5 will find all sources with peak flux densities per beam of 5-sigma or greater). Use the thresh_isl parameter to control how much of each island is used in fitting. Generally, thresh_pix should be larger than thresh_isl.

Only islands with peaks above the absolute threshold will be used. The absolute threshold is calculated as abs_thr = mean + thresh_pix * rms. Use the mean_map and rms_map parameters to control the way the mean and rms are determined.

Adaptive box options

If adaptive_rms_box = True, the rms_box is reduced in size near bright sources and enlarged far from them. This scaling attempts to account for possible strong artifacts around bright sources while still acheiving accurate background rms and mean values when extended sources are present. This option is generally slower than non-adaptive scaling.

Use the rms_box parameter to set the large-scale box and the rms_box_bright parameter to set the small-scale box. The threshold for bright sources can be set with the adaptive_thresh parameter:

adaptive_rms_box ...... True : Use adaptive rms_box when determining rms and mean maps
  adaptive_thresh ..... None : Sources with pixels above adaptive_thresh*
                               clipped_rms will be considered as bright sources (i.e.,
                               with potential artifacts). None => calculate inside
                               program
  rms_box_bright ...... None : Box size, step size for rms/mean map
                               calculation near bright sources. Specify as (box, step)
                               in pixels. None => calculate inside program

adaptive_thresh

This parameter is a float (default is None) that sets the SNR above which sources may be affected by strong artifacts Sources that meet the SNR threshold will use the small-scale box (set by the rms_box_bright parameter) if their sizes at a threshold of 10.0 is less than 25 beam areas.

If None, the threshold is varied from 500 to 50 to attempt to obtain at least 5 candidate bright sources.

rms_box_bright

This parameter is a tuple (default is None) of two integers that sets the box and step sizes to use near bright sources (determined by the adaptive_thresh parameter). The large-scale box size is set with the rms_box parameter.

Advanced options

If advanced_opts = True, a number of additional options are listed. The advanced options do not usually need to be altered from the default values, but can be useful, for example, for fine tuning a fit or for quickly fitting a small region of a much larger image.

The advanced options are:

advanced_opts ......... True : Show advanced options
  aperture ............ None : Radius of aperture in pixels inside which aperture
                               fluxes are measured for each source. None => no aperture
                               fluxes measured
  aperture_posn .. 'centroid': Position the aperture (if aperture is not None) on: 'centroid' or
                               'peak' of the source.
  blank_limit ......... None : Limit in Jy/beam below which pixels are blanked. None => no such
                               blanking is done
  bmpersrc_th ......... None : Theoretical estimate of number of beams per
                               source. None => calculate inside program
  check_outsideuniv .. False : Check for pixels outside the universe
  detection_image ........ '': Detection image file name used only for
                               detecting islands of emission. Source
                               measurement is still done on the main image
  do_cache ........... False : Cache internally derived images to disk
  do_mc_errors ....... False : Estimate uncertainties for 'M'-type sources
                               using Monte Carlo method
  fdr_alpha ........... 0.05 : Alpha for FDR algorithm for thresholds
  fdr_ratio ............ 0.1 : For thresh = None; if #false_pix / #source_pix <
                               fdr_ratio, thresh = 'hard' else thresh = 'fdr'
  fittedimage_clip ..... 0.1 : Sigma for clipping Gaussians while creating fitted
                               image
  fix_to_beam ........ False : Fix major and minor axes and PA of Gaussians to beam?
  group_by_isl ....... False : Group all Gaussians in each island into a single
                               source
  group_method .. 'intensity': Group Gaussians into sources using 'intensity' map or
                               'curvature' map
  group_tol ............ 1.0 : Tolerance for grouping of Gaussians into sources:
                               larger values will result in larger sources
  ini_gausfit ..... 'default': Initial guess for Gaussian parameters: 'default',
                               'fbdsm', or 'nobeam'
  ini_method .... 'intensity': Method by which inital guess for fitting of Gaussians is chosen:
                               'intensity' or 'curvature'
  kappa_clip ........... 3.0 : Kappa for clipped mean and rms
  maxpix_isl .......... None : Maximum number of pixels with emission per island.
                               None -> no limit
  minpix_isl .......... None : Minimum number of pixels with emission per island.
                               None -> calculate inside program
  ncores .............. None : Number of cores to use during fitting, None => use
                               all
  peak_fit ............ True : Find and fit peaks of large islands before fitting
                               entire island
  peak_maxsize ........ 30.0 : If island size in beam area is more than this,
                               attempt to fit peaks separately (if
                               peak_fit=True). Min value is 30
  rmsmean_map_filename  None : Filenames of FITS files to use as the mean and rms maps,
                               given as a list [<mean_map.fits>, <rms_map.fits>]
  rmsmean_map_filename_det  None : Filenames of FITS files to use as the mean and rms maps
                               when a detection image is specified, given as a list
                               [<mean_map.fits>, <rms_map.fits>]
  rms_value ........... None : Value of constant rms in Jy/beam to use if rms_map
                               = False. None => calculate inside program
  spline_rank ............ 3 : Rank of the interpolating function for rms/mean
                               map
  split_isl ........... True : Split island if it is too large, has a large
                               convex deficiency and it opens well. If it doesn't
                               open well, then isl.mean = isl.clipped_mean, and
                               is taken for fitting. Splitting, if needed, is
                               always done for wavelet images
  splitisl_maxsize .... 50.0 : If island size in beam area is more than this,
                               consider splitting island. Min value is 50
  src_ra_dec .......... None : List of source positions at which fitting is done.  E.g.,
                               src_ra_dec = [(197.1932, 47.9188), (196.5573, 42.4852)].
  src_radius_pix ...... None : Radius of the island (if src_ra_dec is not None) in pixels. None
                               => radius is set to the FWHM of the beam major axis.
  stop_at ............. None : Stops after: 'isl' = island finding step or 'read'
                               = image reading step
  trim_box ............ None : Do source detection on only a part of the image.
                               Specify as (xmin, xmax, ymin, ymax) in pixels.
                               E.g., trim_box = (120, 840, 15, 895). None => use
                               entire image

aperture

This parameter is a float (default is None) that sets the radius (in pixels) inside which the aperture flux is measured for each source. The aperture is centered on the either the centroid or the peak of the source (depending on the value of the aperture_posn option). Errors are calculated from the mean of the rms map inside the aperture.

aperture_posn

This parameter is a string (default is 'centroid') that sets the how the aperture is positioned relative to the source. If ‘centroid’, the aperture is centered on the source centroid. If ‘peak’, the aperture is centered on the source peak. If aperture=None (i.e., no aperture radius is specified), this parameter is ignored.

blank_limit

This parameter is a float (default is None) that sets the limit in Jy/beam below which pixels are blanked. All pixels in the ch0 image with a value less than the specified limit and with at least 4 neighboring pixels with values also less than this limit are blanked. If None, any such pixels are left unblanked (and hence will affect the rms and mean maps, etc.). Pixels with a value of NaN are always blanked.

bmpersrc_th

This parameter is a float (default is None) that sets the theoretical estimate of number of beams per source. If None, the value is calculated as N/[n*(alpha-1)], where N is the total number of pixels in the image, n is the number of pixels in the image whose value is greater than 5 times the clipped rms, and alpha is the slope of the differential source counts distribution, assumed to be 2.5.

The value of bmpersrc_th is used to estimate the average separation in pixels between two sources, which in turn is used to estimate the boxsize for calculating the background rms and mean images. In addition, if the value is below 25 (or the ratio of clipped mean to clipped rms of the image is greater than 0.1), the image is assumed to be confused and hence the background mean is put to zero.

check_outsideuniv

This parameter is a Boolean (default is False). If True, then the coordinate of each pixel is examined to check if it is outside the universe, which may happen when, e.g., an all sky image is made with SIN projection (commonly done at LOFAR earlier). When found, these pixels are blanked (since imaging software do not do this on their own). Note that this process takes a lot of time, as every pixel is checked in case weird geometries and projections are used.

detection_image

This parameter is a string (default is '') that sets the detection image file name used only for detecting islands of emission. Source measurement is still done on the main image. The detection image can be a FITS or CASA 2-, 3-, or 4-D cube and must have the same size and WCS parameters as the main image.

do_cache

This parameter is a Boolean (default is False) that controls whether internally derived images are stored in memory or are cached to disk. Caching can reduce the amount of memory used, and is therefore useful when analyzing large images.

do_mc_errors

This parameter is a Boolean (default is False). If True, uncertainties on the sizes and positions of ‘M’-type sources due to uncertainties in the constituent Gaussians are estimated using a Monte Carlo technique. These uncertainties are added in quadrature with those calculated using Condon (1997). If False, these uncertainties are ignored, and errors are calculated using Condon (1997) only.

Enabling this option will result in longer run times if many ‘M’-type sources are present, but should give better estimates of the uncertainites, particularly for complex sources composed of many Gaussians.

fdr_alpha

This parameter is a float (default is 0.05) that sets the value of alpha for the FDR algorithm for thresholding. If thresh is 'fdr', then the estimate of fdr_alpha (see Hopkins et al. 2002 [3] for details) is stored in this parameter.

fdr_ratio

This parameter is a float (default is 0.1). When thresh = None, if #false_pix / #source_pix < fdr_ratio, thresh = 'hard' otherwise thresh = 'fdr'.

fittedimage_clip

This parameter is a float (default is 0.1). When the residual image is being made after Gaussian decomposition, the model images for each fitted Gaussian are constructed up to a size 2b, such that the amplitude of the Gaussian falls to a value of fitted_image_clip times the local rms, b pixels from the peak.

fix_to_beam

This parameter is a Boolean (default is False). If True, then during fitting the major and minor axes and PA of the Gaussians are fixed to the beam. Only the amplitude and position are fit. If False, all parameters are fit. Note that when this option is activated, as a consequence of using fewer free parameters, the estimated errors on the peak and total flux densities are a factor of sqrt(2) lower compared to the case in which all parameters are fit (see Condon 1997). Additionally, the reported errors on the major and minor axes and the PA are zero.

group_by_isl

This parameter is a Boolean (default is False). If True, all Gaussians in the island belong to a single source. If False, grouping is controlled by the group_tol parameter.

group_method

This parameter is a string (default is 'intensity'). Gaussians are deemed to be a part of the same source if: (1) no pixel on the line joining the centers of any pair of Gaussians has a (Gaussian-reconstructed) value less than the island threshold, and (2) the centers are separated by a distance less than half the sum of their FWHMs along the line joining them. If 'curvature', the above comparisons are done on the curature map (see Hancock et al. 2012). If 'intensity', the comparisons are done on the intensity map.

group_tol

This parameter is a float (default is 1.0) that sets the tolerance for grouping of Gaussians into sources: larger values will result in larger sources. Sources are created by grouping nearby Gaussians as follows: (1) If the difference between the minimum value between two Gaussians and the lower of the peak flux densities of the Gaussians in an island is less than group_tol * thresh_isl * rms_clip, and (2) if the centers are separated by a distance less than 0.5 * group_tol of the sum of their FWHMs along the PA of the line joining them, they belong to the same island.

ini_gausfit

This parameter is a string (default is 'default'). These are three different ways of estimating the initial guess for fitting of Gaussians to an island of emission. If 'default', the maximum number of Gaussians is estimated from the number of peaks in the island. An initial guess is made for the parameters of these Gaussians before final fitting is done. This method should produce the best results when there are no large sources present. If 'simple', the maximum number of Gaussians per island is set to 25, and no initial guess for the Gaussian parameters is made. Lastly, the 'nobeam' method is similar to the 'default' method, but no information about the beam is used. This method is best used when source sizes are expected to be very different from the beam and is generally slower than the other methods. For wavelet images, the value used for the original image is used for wavelet order j <= 3 and 'nobeam' for higher orders.

ini_method

This parameter is a string (default is 'intensity'). If 'intensity', the inital guess described in the help for the ini_gausfit parameter is calculated using the intensity (ch0) image. If 'curvature', it is done using the curvature map (see Hancock et al. 2012).

kappa_clip

This parameter is a float (default is 3.0) that is the factor used for Kappa-alpha clipping, as in AIPS. For an image with few source pixels added on to (Gaussian) noise pixels, the dispersion of the underlying noise will need to be determined. This is done iteratively, whereby the actual dispersion is first computed. Then, all pixels whose value exceeds kappa clip times this rms are excluded and the rms is computed again. This process is repeated until no more pixels are excluded. For well behaved noise statistics, this process will converge to the true noise rms with a value for this parameter ~3-5. A large fraction of source pixels, less number of pixels in total, or significant non-Gaussianity of the underlying noise will all lead to non-convergence.

maxpix_isl

This parameter is an integer (default is None) that sets the maximum number of pixels in an island for the island to be included. If None, there is no limit.

minpix_isl

This parameter is an integer (default is None) that sets the minimum number of pixels in an island for the island to be included. If None, the number of pixels is set to 1/3 of the area of an unresolved source using the beam and pixel size information in the image header. It is set to 6 pixels for all wavelet images.

ncores

This parameter is an integer (default is None) that sets the number of cores to use during fitting. If None, all available cores are used (one core is reserved for plotting).

peak_fit

This parameter is a Boolean (default is True). When True, PyBDSF will identify and fit peaks of emission in large islands iteratively (the size of islands for which peak fitting is done is controlled with the peak_maxsize option), using a maximum of 10 Gaussians per iteration. Enabling this option will generally speed up fitting (by factors of many for large islands), but may result in somewhat higher residuals.

peak_maxsize

This parameter is a float (default is 30.0). If island size in beam area is more than this value, attempt to fit peaks iteratively (if peak_fit = True). The minimum value is 30.

rmsmean_map_filename

This parameter is a list (default is None) that sets the filenames of FITS files to use as the mean and rms maps, given as a list [<mean_map.fits>, <rms_map.fits>]. If supplied, the internally generated mean and rms maps are not used.

rmsmean_map_filename_det

This parameter is a list (default is None) that sets the filenames of FITS files to use as the mean and rms maps when a detection image is specified, given as a list [<mean_map.fits>, <rms_map.fits>]. If supplied, the internally generated mean and rms maps are not used.

rms_value

This parameter is a float (default is None) that sets the value of constant rms in Jy/beam to use if rms_map = False. If None, the value is calculated inside the program.

spline_rank

This parameter is an integer (default is 3) that sets the order of the interpolating spline function to interpolate the background rms and mean maps over the entire image.

Note

Bicubic interpolation (the default) can cause ringing artifacts to appear in the rms and mean maps in regions where sharp changes occur. These artifacts can result in regions with negative values. If you find such artifacts, try changing the spline_rank parameter.

split_isl

This parameter is a Boolean (default is True). If True, an island is split if it is too large, has a large convex deficiency and it opens well. If it doesn’t open well, then isl.mean = isl.clipped_mean, and is taken for fitting. Splitting, if needed, is always done for wavelet images

splitisl_maxsize

This parameter is a float (default is 50.0). If island size in beam area is more than this, consider splitting island. The minimum value is 50.

src_ra_dec

This parameter is a list of tuples (default is None) that defines the center positions at which fitting will be done. The size of the region used for the fit is given by the src_radius_pix parameter. Positions should be given as a list of RA and Dec, in degrees, one set per source. These positions will override the normal island finding module.

src_radius_pix

This parameter is a float (default is None) that determines the size of the region used to fit the source positions specified by the src_ra_dec parameter. If None, the radius is set to the FWHM of the beam major axis.

stop_at

This parameter is a string (default is None) that stops an analysis after: ‘isl’ = island finding step or ‘read’ = image reading step.

trim_box

This parameter is a tuple (default is None) that defines a subregion of the image on which to do source detection. It is specified as (xmin, xmax, ymin, ymax) in pixels. E.g., trim_box = (120, 840, 15, 895). If None, the entire image is used.

Flagging options

If flagging_opts = True, a number of options are listed for flagging unwanted Gaussians that occur durring a fit. Flagged Gaussians are not included in any further analysis or catalog. They may be viewed using the show_fit task (see show_fit: visualizing the fit results). A flag value is associated with each flagged Gaussian that allows the user to determine the reason or reasons that it was flagged. If multiple flagging conditions are met by a single Gaussian, the flag values are summed. For example, if a Gaussian is flagged because it is too large (its size exceeds that implied by flag_maxsize_bm, giving a flag value of 64) and because it is too bright (its peak flux density per beam exceeds that implied by flag_maxsnr, giving a flag value of 2) then the final flag value is 64 + 2 = 66.

Note

If a fit did not produce good results, it is often useful to check whether there are flagged Gaussians and adjust the flagging options as necessary. Flagged Gaussians can be viewed by setting ch0_flagged = True in the show_fit task.

The options for flagging of Gaussians are:

flagging_opts ......... True : Show options for Gaussian flagging
  flag_bordersize ........ 0 : Flag Gaussian if centre is outside border -
                               flag_bordersize pixels
  flag_maxsize_bm ..... 25.0 : Flag Gaussian if area greater than flag_maxsize_bm
                               times beam area
  flag_maxsize_isl ..... 1.0 : Flag Gaussian if x, y bounding box around
                               sigma-contour is factor times island bbox
  flag_maxsnr .......... 1.5 : Flag Gaussian if peak is greater than flag_maxsnr
                               times max value in island
  flag_minsize_bm ...... 0.7 : Flag Gaussian if flag_smallsrc = True and area
                               smaller than flag_minsize_bm times beam area
  flag_minsnr .......... 0.9 : Flag Gaussian if peak is less than flag_minsnr
                               times thresh_pix times local rms
  flag_smallsrc ...... False : Flag sources smaller than flag_minsize_bm times
                               beam area

flag_bordersize

This parameter is an integer (default is 0). Any fitted Gaussian whose centre is flag_bordersize pixels outside the island bounding box is flagged. The flag value is increased by 4 (for x) and 8 (for y).

flag_maxsize_bm

This parameter is a float (default is 25.0). Any fitted Gaussian whose size is greater than flag_maxsize_bm times the synthesized beam is flagged. The flag value is increased by 64.

flag_maxsize_fwhm

This parameter is a float (default is 0.3). Any fitted Gaussian whose contour of flag_maxsize_fwhm times the FWHM falls outside the island is flagged. The flag value is increased by 256.

flag_maxsize_isl

This parameter is a float (default is 1.0). Any fitted Gaussian whose maximum x-dimension is larger than flag_maxsize_isl times the x-dimension of the island (and likewise for the y-dimension) is flagged. The flag value is increased by 16 (for x) and 32 (for y).

flag_maxsnr

This parameter is a float (default is 1.5). Any fitted Gaussian whose peak is greater than flag_maxsnr times the value of the image at the peak of the Gaussian is flagged. The flag value is increased by 2.

flag_minsize_bm

This parameter is a float (default is 0.7). If flag_smallsrc is True, then any fitted Gaussian whose size is less than flag_maxsize_bm times the synthesized beam is flagged. The Gaussian flag is increased by 128.

flag_minsnr

This parameter is a float (default is 0.7). Any fitted Gaussian whose peak is less than flag_minsnr times thresh_pix times the local rms is flagged. The flag value is increased by 1.

flag_smallsrc

This parameter is a Boolean (default is False). If True, then fitted Gaussians whose size is less than flag_minsize_bm times the synthesized beam area are flagged. When combining Gaussians into sources, an error is raised if a 2x2 box with the peak of the Gaussian does not have all four pixels belonging to the source. Usually this means that the Gaussian is an artifact or has a very small size.

If False, then if either of the sizes of the fitted Gaussian is zero, then the Gaussian is flagged.

If the image is barely Nyquist sampled, this flag is best set to False. This flag is automatically set to False while decomposing wavelet images into Gaussians.

Output options

If output_opts = True, options to control the output generated by process_image are listed. By default, only a log file is generated and output is controlled with the export_image (see export_image: exporting internally derived images) and write_catalog (see write_catalog: writing source catalogs) tasks. However, the user can specify that a number of optional output files be made automatically whenever process_image is run. These options are most useful for debugging or when running PyBDSF non-interactively in a Python script (see Scripting PyBDSF).

The output options are:

output_opts ........... True : Show output options
  bbs_patches ......... None : For BBS format, type of patch to use: None => no
                               patches. 'single' => all Gaussians in one patch.
                               'gaussian' => each Gaussian gets its own patch.
                               'source' => all Gaussians belonging to a single
                               source are grouped into one patch. 'mask' => use mask
                               file specified by bbs_patches_mask
  bbs_patches_mask .... None : Name of the mask file (of same size as input image)
                               that defines the patches if bbs_patches = 'mask'
  indir ............... None : Directory of input FITS files. None => get from
                               filename
  opdir_overwrite .. 'overwrite': 'overwrite'/'append': If output_all=True,
                               delete existing files or append a new directory
  outdir .............. None : Directory to use for all output files
                               (including log files). None => parent directory of the
                               input filename
  output_all ......... False : Write out all files automatically to directory
                               'outdir/filename_pybdsf'
  plot_allgaus ....... False : Make a plot of all Gaussians at the end
  plot_islands ....... False : Make separate plots of each island during fitting
                               (for large images, this may take a long time and a
                               lot of memory)
  print_timing ....... False : Print basic timing information
  quiet .............. False : Suppress text output to screen. Output is still
                               sent to the log file as usual
  savefits_det_meanim  False : Save detection background mean image as fits file
  savefits_det_rmsim . False : Save detection background rms image as fits file
  savefits_meanim .... False : Save background mean image as fits file
  savefits_modelim ... False : Save Gaussian model image as fits file
  savefits_normim .... False : Save norm image as fits file
  savefits_rankim .... False : Save island rank image as fits file
  savefits_residim ... False : Save residual image as fits file
  savefits_rmsim ..... False : Save background rms image as fits file
  solnname ............ None : Name of the run, to be appended to the name of the
                               output directory
  verbose_fitting .... False : Print out extra information during fitting

bbs_patches

This parameter is a string (default is None) that sets the type of patch to use in BBS-formatted catalogs. When the Gaussian catalogue is written as a BBS-readable sky file, this option determines whether all Gaussians are in a single patch ('single'), there are no patches (None), all Gaussians for a given source are in a separate patch ('source'), each Gaussian gets its own patch ('gaussian'), or a mask image is used to define the patches ('mask').

If you wish to have patches defined by island, then set group_by_isl = True before fitting to force all Gaussians in an island to be in a single source. Then set bbs_patches = 'source' when writing the catalog.

bbs_patches_mask

This parameter is a string (default is None) that sets the file name of the mask file to use to define patches in BBS-formatted catalogs. The mask image should be 1 inside the patches and 0 elsewhere and should be the same size as the input image (before any trim_box is applied). Any Gaussians that fall outside of the patches will be ignored and will not appear in the output sky model.

indir

This parameter is a string (default is None) that sets the directory of input FITS files. If None, the directory is defined by the input filename.

opdir_overwrite

This parameter is a string (default is 'overwrite') that determines whether existing output files are overwritten or not.

outdir

This parameter is a string (default is None) that sets the directory to use for all output files (including log files). If None, the parent directory of the input image filename is used.

output_all

This parameter is a Boolean (default is False). If True, all output products are written automatically to the directory 'outdir/filename_pybdsf'.

plot_allgaus

This parameter is a Boolean (default is False). If True, make a plot of all Gaussians at the end.

plot_islands

This parameter is a Boolean (default is False). If True, make separate plots of each island during fitting (for large images, this may take a long time and a lot of memory).

print_timing

This parameter is a Boolean (default is False). If True, print basic timing information.

quiet

This parameter is a Boolean (default is False). If True, suppress text output to screen. Output is still sent to the log file as usual.

savefits_det_rmsim

This parameter is a Boolean (default is False). If True, save detection background rms image as a FITS file.

savefits_det_meanim

This parameter is a Boolean (default is False). If True, save detection background mean image as a FITS file.

savefits_meanim

This parameter is a Boolean (default is False). If True, save background mean image as a FITS file.

savefits_modelim

This parameter is a Boolean (default is False). If True, save Gaussian model image as a FITS file.

savefits_normim

This parameter is a Boolean (default is False). If True, save norm image as a FITS file.

savefits_rankim

This parameter is a Boolean (default is False). If True, save island rank image as a FITS file.

savefits_residim

This parameter is a Boolean (default is False). If True, save residual image as a FITS file.

savefits_rmsim

This parameter is a Boolean (default is False). If True, save background rms image as a FITS file.

solnname

This parameter is a string (default is None) that sets the name of the run, to be appended to the name of the output directory.

verbose_fitting

This parameter is a Boolean (default is False). If True, print out extra information during fitting.

Multichannel options

If multichan_opts = True, the options used to control the way PyBDSF handles images with more than one frequency channel are listed. In particular, these options control how the multichannel image is collapsed to form the ch0 image on which source detection is done.

The options concerning multichannel images are:

multichan_opts ........ True : Show options for multi-channel images
  beam_sp_derive ...... True : If True and beam_spectrum is None, then assume
                               header beam is for lowest frequency and scales
                               with frequency for channels
  beam_spectrum ....... None : FWHM of synthesized beam per channel. Specify as
                               [(bmaj_ch1, bmin_ch1, bpa_ch1), (bmaj_ch2,
                               bmin_ch2, bpa_ch2), etc.] in degrees. E.g.,
                               beam_spectrum = [(0.01, 0.01, 45.0), (0.02, 0.01,
                               34.0)] for two channels. None => all equal to beam
  collapse_av ........... [] : List of channels to average if collapse_mode =
                               'average'; None => all
  collapse_ch0 ........... 0 : Number of the channel for source extraction, if
                               collapse_mode = 'single'
  collapse_mode ... 'average': Collapse method: 'average' or 'single'. Average
                               channels or take single channel to perform source
                               detection on
  collapse_wt ....... 'unity': Weighting: 'unity' or 'rms'. Average channels with
                               weights = 1 or 1/rms_clip^2 if collapse_mode =
                               'average'
  frequency_sp ........ None : Frequency in Hz of channels in input image when
                               more than one channel is present. E.g., frequency
                               = [74e6, 153e6]. None => get from header

beam_sp_derive

This parameter is a Boolean (default is True). If True and the parameter beam_spectrum is None and no channel-dependent beam parameters are found in the header, then we assume that the beam in the header is for the lowest frequency of the image cube and scale accordingly to calculate the beam per channel. If False, then a constant value of the beam (that given in the header or specified with the beam parameter) is taken instead.

Note

This parameter is used only in the Spectral index module module.

beam_spectrum

This parameter is a list of tuples (default is None) that sets the FWHM of synthesized beam per channel. Specify as [(bmaj_ch1, bmin_ch1, bpa_ch1), (bmaj_ch2, bmin_ch2, bpa_ch2), etc.] in degrees. E.g., beam_spectrum = [(0.01, 0.01, 45.0), (0.02, 0.01, 34.0)] for two channels.

If None, then the channel-dependent restoring beam is searched for in the header. If not found, the beam is either assumed to be a constant or to scale with frequency, depending on whether the parameter beam_sp_derive is False or True.

Note

This parameter is used only in the Spectral index module module. For normal fitting (outside of the spectral index module), a constant value of the beam (that given in the header or specified with the beam parameter) is taken instead.

collapse_av

This parameter is a list of integers (default is []) that specifies the channels to be averaged to produce the continuum image for performing source extraction, if collapse_mode is 'average'. If the value is [], then all channels are used. Otherwise, the value is a Python list of channel numbers.

collapse_ch0

This parameter is an integer (default is 0) that specifies the number of the channel for source extraction, if collapse_mode = 'single'.

collapse_mode

This parameter is a string (default is 'average') that determines whether, when multiple channels are present, the source extraction is done on a single channel ('single') or an average of many channels ('average').

collapse_wt

This parameter is a string (default is 'unity'). When collapse_mode is 'average', then if this value is 'unity', the channels given by collapse_av are averaged with unit weights and if 'rms', then they are averaged with weights which are inverse square of the clipped rms of each channel image.

frequency_sp

This parameter is a list of floats (default is None) that sets the frequency in Hz of channels in input image when more than one channel is present. E.g., frequency_sp = [74e6, 153e6].

If the frequency is not given by the user, then it is looked for in the image header. If not found, then an error is raised. PyBDSF will not work without the knowledge of the frequency.

À trous wavelet decomposition module

If atrous_do = True, this module decomposes the residual image that results from the normal fitting of Gaussians into wavelet images of various scales. Such a decomposition is useful if there is extended emission that is not well fit during normal fitting. Such emission therefore remains in the Gaussian residual image and can be further fit by Gaussians whose size is tuned to the various wavelet scales. Therefore, wavelet decomposition should be used when there is significant residual emission that remains after normal Gaussian fitting.

The wavelet module performs the following steps:

The number of wavelet scales to be considered is set by the atrous_jmax parameter. By default, this number is determined automatically from the size of the largest island in the image. Wavelet images are then made for scales of order (j) ranging from 1 to jmax.
For each scale (j), the appropriate à trous wavelet transformation is made (see Holschneider et al. 1989 for details). Additionally, the “remainder” image (called the c_J image) is also made. This image includes all emission not included in the other wavelet images.
Depending on the value of the atrous_sum option, fitting is done to either an image that is a sum over all scales equal to or larger than the scale under consideration (atrous_sum = True) or to an image of a single scale (atrous_sum = False). Fitting to the sum over all larger scales will generally result in increased signal to noise.
If atrous_bdsm = True, an rms map and a mean map are made for each wavelet image and Gaussians are fit in the normal way. Gaussians can be optionally restricted to lie within islands found from the initial image. If a wavelet island overlaps spatially with an existing island, the two islands are merged together to form a single island. The wavelet Gaussians can then be included in source catalogs (see write_catalog: writing source catalogs).

The options for this module are as follows:

atrous_do ............. True : Decompose Gaussian residual image into multiple
                               scales
  atrous_bdsm_do ...... True : Perform source extraction on each wavelet scale
  atrous_jmax ............ 0 : Max allowed wavelength order, 0 => calculate
                               inside program
  atrous_lpf ........... 'b3': Low pass filter, either 'b3' or 'tr', for B3
                               spline or Triangle
  atrous_orig_isl .... False : Restrict wavelet Gaussians to islands found in
                               original image
  atrous_sum .......... True : Fit to the sum of images of the remaining wavelet
                               scales
  use_scipy_fft ....... True : Use fast SciPy FFT for convolution

atrous_bdsm_do

This parameter is a Boolean (default is False). If True, PyBDSF performs source extraction on each wavelet scale.

atrous_jmax

This parameter is an integer (default is 0) which is the maximum order of the à trous wavelet decomposition. If 0 (or <0 or >15), then the value is determined within the program. The value of this parameter is then estimated as the (lower) rounded off value of ln[(nm-l)/(l-1) + 1]/ln2 + 1 where nm is the minimum of the residual image size (n, m) in pixels and l is the length of the filter à trous lpf (see the atrous_lpf parameter for more info).

A sensible value is such that the size of the kernel is not more than 3-4 times smaller than the smallest image dimension.

atrous_lpf

This parameter is a string (default is 'b3') that sets the low pass filter, which can be either the B3 spline or the triangle function, which is used to generate the à trous wavelets. The B3 spline is [1, 4, 6, 4, 1] and the triangle is [1, 2, 1], normalised so that the sum is unity. The lengths of the filters are hence 5 and 3 respectively.

atrous_orig_isl

This parameter is a Boolean (default is False). If True, all wavelet Gaussians must lie within the boundaries of islands found in the original image. If False, new islands that are found only in the wavelet images are included in the final fit.

atrous_sum

This parameter is a Boolean (default is True). If True, fitting is done on an image that is the sum of the remaining wavelet scales. Using the sum will generally result in improved signal. If False, fitting is done on only the wavelet scale under consideration.

use_scipy_fft

This parameter is a Boolean (default is True). If True, the SciPy FFT function will be used instead of the custom version. The SciPy version is much faster but also uses much more memory.

PSF variation module

If psf_vary_do = True, then the spatial variations in the PSF are estimated and their effects corrected for. To this end, PyBDSF performs the following steps:

A list of sources that are likely to be unresolved is constructed. This is done by first selecting only type ‘S’ sources by default (see Definition of output columns for details of source types), but this restriction can be overridden using the psf_stype_only option) and sources with SNRs that exceed psf_snrcut. Next, a function is fit to determine how the size of sources (normalized by the median size) varies with the SNR. The function used is defined as \(\sigma / median = \sqrt(c_1^2 + c_2^2/SNR^2)\), where \(\sigma\) is the size of the Gaussian and \(c_1\) and \(c_2\) are free parameters. Clipping of outliers is done during this fitting, controlled by the psf_nsig parameter. Lastly, unresolved sources are selected by choosing sources that lie within psf_kappa2 times the rms of this best-fit sigma-SNR relation. As this last step can be unreliable for high-SNR sources, an additional selection can be made for the highest SNR sources using the psf_high_snr parameter. All sources with SNRs above psf_high_snr will be taken as unresolved.
Next the image is tessellated using Voronoi tessellation to produce tiles within which the PSF shape is calculated (and assumed to be constant). The list of probable unresolved sources is filtered to select “calibrator” sources to use to determine the tessellation tiles. These sources are the brightest sources (known as the primary generators), defined as those sources that have SNRs in the top fraction of sources defined by psf_snrtop and that also have SNRs greater than psf_snrcutstack. These sources are then grouped by their proximity, if they are within 50% of the distance to third closest source.
The unresolved sources within each tile that have SNRs greater than psf_snrcutstack are then stacked to form a high-SNR PSF. For each tile, this PSF is fit with a Gaussian to recover its size. The significance of the variation in the sizes across the image is quantified.
If the variation is significant, the major axis, minor axis, and position angle are then interpolated across the image. Smoothing can be applied to these images to smooth out artifacts due to noise and the interpolation. Additionally, images are made of the ratio of peak-to-total flux and peak-to-aperture flux (if an aperture is specified). These ratio images provide conversions from total flux to peak flux for point sources. In the absence of smearing effects, these ratios should be around unity. However, if ionospheric effects are present, significant smearing can be present. In this case, these ratio images can be useful, for example, in determining the sensitivity at a particular location in the image to a point source with a given total flux.
Lastly, the deconvolved source sizes are adjusted to include the PSF variation as a function of position.

The options for this module are as follows:

psf_vary_do ........... True : Calculate PSF variation across image
  psf_high_snr ........ None : SNR above which all sources are taken to be
                               unresolved. E.g., psf_high_snr = 20.0. None => no
                               such selection is made
  psf_itess_method ....... 0 : 0 = normal, 1 = 0 + round, 2 = LogSNR, 3 =
                               SqrtLogSNR
  psf_kappa2 ........... 2.0 : Kappa for clipping for analytic fit
  psf_nsig ............. 3.0 : Kappa for clipping within each bin
  psf_over ............... 2 : Factor of nyquist sample for binning bmaj, etc. vs
                               SNR
  psf_smooth .......... None : Size of Gaussian to use for smoothing of
                               interpolated images in arcsec. None => no smoothing
  psf_snrcut .......... 10.0 : Minimum SNR for statistics
  psf_snrcutstack ..... 15.0 : Unresolved sources with higher SNR taken for
                               stacked psfs
  psf_snrtop .......... 0.15 : Fraction of SNR > snrcut as primary generators
  psf_stype_only ...... True : Restrict sources used in PSF variation
                               estimating to be only of type 'S'

psf_high_snr: This parameter is a float (default is None). Gaussians with SNR greater than this are used to determine the PSF variation, even if they are deemed to be resolved. This corrects for the unreliability at high SNRs in the algorithm used to find unresolved sources. The minimum value is 20.0. If None, then no such selection is made.
psf_itess_method: This parameter is an integer (default is 0) which can be 0, 1, 2 or 3, which corresponds to a tessellation method. If 0, 2 or 3, then the weights used for Voronoi tessellation are unity, log(SNR) and sqrt[log(SNR)] where SNR is the signal to noise ratio of the generator in a tile. If 1, then the image is tessellated such that each tile has smooth boundaries instead of straight lines, using pixel-dependent weights.
psf_kappa2: This parameter is a float (default is 2.0). When iteratively arriving at a statistically probable set of ‘unresolved’ sources, the fitted major and minor axis sizes versus SNR are binned and fitted with analytical functions. Those Gaussians which are within psf_kappa2 times the fitted rms from the fitted median are then considered ‘unresolved’ and are used further to estimate the PSFs.
psf_nsig: This parameter is a float (default is 3.0). When constructing a set of ‘unresolved’ sources for psf estimation, the (clipped) median, rms and mean of major and minor axis sizes of Gaussians versus SNR within each bin is calculated using kappa = psf_nsig.
psf_over: This parameter is an integer (default is 2). When constructing a set of ‘unresolved’ sources for psf estimation, this parameter controls the factor of nyquist sample for binning bmaj, etc. vs SNR.
psf_smooth: This parameter is a float (default is None) that sets the smoothing scale (in arcsec) used to smooth the interpolated images. Generally, artifacts due to noise and the interpolation can be significantly reduced if the smoothing scale is similar to the typical source separation scale.
psf_snrcut: This parameter is a float (default is 10.0). Only Gaussians with SNR greater than this are considered for processing. The minimum value is 5.0
psf_snrcutstack: This parameter is a float (default is 15.0). Only Gaussians with SNR greater than this are used for estimating PSF images in each tile.
psf_snrtop: This parameter is a float (default is 0.15). If psf_generators is ‘calibrator’, then the peak pixels of Gaussians which are the psf_snrtop fraction of the SNR distribution are taken as Voronoi generators.
psf_stype_only: This parameter is a Boolean (default is False). If True, sources are restricted to be only of type ‘S’.

Spectral index module

If spectralindex_do = True (and the input image has more than one frequency), then spectral indices are calculated for the sources in the following way:

The rms maps for the remaining channels are determined.
Neighboring channels are averaged to attempt to obtain the target SNR per channel for a given source, set by the specind_snr parameter.

Note

No color corrections are applied during averaging. However, unless the source spectral index is very steep or the channels are very wide, the correction is minimal. See Effect of neglecting color corrections for details.
Flux densities are measured for both individual Gaussians and for total sources. Once source flux densities have been measured in each channel, the SEDs are fit with the following power-law model: \(S_\nu \propto \nu^\alpha\), where \(\alpha\) is the spectral index. The resulting best-fit spectral index is then included in any catalogs that are written out (see write_catalog: writing source catalogs). In addition, plots of the fits can be viewed with the show_fit task (see show_fit: visualizing the fit results).

The options for this module are as follows:

spectralindex_do ...... True : Calculate spectral indices (for multi-channel
                               image)
  flagchan_list ......... [] : List of channels to flag before (averaging and)
                               extracting spectral index
  flagchan_rms ........ True : Flag channels before (averaging and) extracting
                               spectral index, if their rms if more than 5
                               (clipped) sigma outside the median rms over all
                               channels, but only if <= 10% of channels
  flagchan_snr ........ True : Flag channels that do not meet SNR criterion set
                               by specind_snr
  specind_maxchan ........ 0 : Maximum number of channels to average for a
                               given source when when attempting to meet target
                               SNR. 1 => no averaging; 0 => no maximum
  specind_snr .......... 3.0 : Target SNR to use when fitting power law. If
                               there is insufficient SNR, neighboring channels
                               are averaged to obtain the target SNR

flagchan_list: This parameter is a list of integers (default is []) that specifies the channels to flag before (averaging and) extracting the spectral indices. Flagged channels will not be used during fitting. If the value is an empty list ([]), then all channels are used. Otherwise, the value is a Python list of channel numbers, starting from 0 (i.e., the first channel has number 0, the second has number 1, etc.).
flagchan_rms: This parameter is a Boolean (default is True). If True, then the clipped rms and median (r and m) of the clipped rms of each channel is calculated. Those channels whose clipped rms is greater than 4r away from m are flagged prior to averaging and calculating spectral indices from the image cube. However, these channels are flagged only if the total number of these bad channels does not exceed 10% of the total number of channels themselves.
flagchan_snr: This parameter is a Boolean (default is True). If True, then flux densities in channels that do not meet the target SNR are not used in fitting.
specind_maxchan: This parameter is an integer (default is 0) that sets the maximum number of channels that can be averaged together to attempt to reach the target SNR set by the specind_snr parameter. If 0, there is no limit to the number of channels that can be averaged. If 1, no averaging will be done.
specind_snr: This parameter is a float (default is 3.0) that sets the target SNR to use when fitting for the spectral index. If there is insufficient SNR, neighboring channels are averaged to obtain the target SNR. The maximum allowable number of channels to average is determined by the specind_maxchan parameter. Channels (after averaging) that fail to meet the target SNR are not used in fitting.

Polarization module

If polarisation_do = True, then the polarization properties of the sources are calculated. First, if pi_fit = True, source detection is performed on the polarized intensity (PI) image [4] to detect sources without Stokes I counterparts. The polarization module then calculates the I, Q, U, and V flux densities, the total, linear, and circular polarisation fractions and the linear polarisation angle of each Gaussian and source. The linear polarisation angle is defined from North, with positive angles towards East. Flux densities are calculated by fitting the normalization of the Gaussians found using the Stokes I or PI images.

For linearly polarised emission, the signal and noise add vectorially, giving a Rice distribution instead of a Gaussian one. To correct for this, a bias is estimated and removed from the polarisation fraction using the same method used for the NVSS catalog (see ftp://ftp.cv.nrao.edu/pub/nvss/catalog.ps). Errors on the linear and total polarisation fractions and polarisation angle are estimated using the debiased polarised flux density and standard error propagation. See Sparks & Axon (1999) [5] for a more detailed treatment.

The options for this module are as follows:

polarisation_do ....... True : Find polarisation properties
  pi_fit .............. True : Check the polarized intesity (PI) image for
                               sources not found in Stokes I
  pi_thresh_isl ....... None : Threshold for PI island boundary in number
                               of sigma above the mean. None => use thresh_isl
  pi_thresh_pix ....... None : Source detection threshold for PI image:
                               threshold for the island peak in number of sigma
                               above the mean. None => use thresh_pix

pi_fit: This parameter is a Boolean (default is True). If True, the polarized intensity image is searched for sources not present in the Stokes I image. If any such sources are found, they are added to the the Stokes I source lists. Use the pi_thresh_pix and pi_thresh_isl parameters to control island detection in the PI image.
pi_thresh_isl: This parameter is a float (default is None) that determines the region to which fitting is done in the polarized intensity (PI) image. If None, the value is set to that of the thresh_isl parameter. A higher value will produce smaller islands, and hence smaller regions that are considered in the fits. A lower value will produce larger islands. Use the pi_thresh_pix parameter to set the detection threshold for sources. Generally, pi_thresh_isl should be lower than pi_thresh_pix.
pi_thresh_pix: This parameter is a float (default is None) that sets the overall detection threshold for islands in the polarized intensity (PI) image (i.e. pi_thresh_pix = 5 will find all sources with peak flux densities per beam of 5-sigma or greater). If None, the value is set to that of the thresh_pix parameter. Use the pi_thresh_isl parameter to control how much of each island is used in fitting. Generally, pi_thresh_pix should be larger than pi_thresh_isl.

Shapelet decomposition module

If shapelet_do = True, then islands are decomposed into shapelets. Shapelets are a set of 2-D basis functions (for details, see Refregier 2003 [6]) that can be used to completely model any source, typically with far fewer parameters than pixels in the source. Shapelets are useful in particular for modeling complex islands that are not well modeled by Gaussians alone. PyBDSF can currently fit cartesian shapelets to an image. The shapelet parameters can be written to a catalog using write_catalog (see write_catalog: writing source catalogs).

For each island of emission, a shapelet decomposition is done after estimating the best values of the center, the scale \(\beta\), and nmax in the following way. First, an initial guess of \(\beta\) is taken as \(2\sqrt{[m2(x)m2(y)]}\), where \(m2\) is the second moment over the island, based on shapeelt analysis of simulated images of resolved sources. Similarly, a guess for nmax is taken as the minimum of 14, and maximum of 10 and \(2n + 2\) where \(n=\sqrt{(n^2 + m^2)}/n_p^n - 1\), where (n, m) is the size of the island and \(n^m_p\) is the synthesized beam minor axis FWHM in pixels. This guess for nmax is based partly on simulations and partly on the requirememts of computing time, number of constraints, etc, for shapelet decomposition.

These initial values are then used to calculate the optimal central position around which to decompose the island. First, for every pixel in the island, the coefficients c12 and c21 are computed assuming that pixel as the centre of expansion. Next, the zero crossings for every vertical line of the c12 image and horizontal line of the c21 image are computed. The intersection point of these two zero-crossing vectors is then taken as the proper centre of the expansion for the image. If this procedure does not work, then the first moment is taken as the center.

This updated center position is used to compute the optimal \(\beta\), which is taken as the value of \(\beta\) that minimises the residual rms in the island area. Using this \(\beta\), the center is computed once more and the final shapelet deocmposition is then made.

The options for this module are as follows:

shapelet_do ........... True : Decompose islands into shapelets
  shapelet_basis .. 'cartesian': Basis set for shapelet decomposition:
                               'cartesian' or 'polar'
  shapelet_fitmode .... 'fit': Calculate shapelet coeff's by fitting ('fit') or
                               integrating (None)
  shapelet_gresid ..... 'False': Use Gaussian residual image for shapelet
                               decomposition?

shapelet_basis: This parameter is a string (default is 'cartesian') that determines the type of shapelet basis used. Currently however, only cartesian is supported.
shapelet_fitmode: This parameter is a string (default is 'fit') that determines the method of calculating shapelet coefficients. If None, then these are calculated by integrating (actually, by summing over pixels, which introduces errors due to discretisation). If ‘fit’, then the coefficients are found by least-squares fitting of the shapelet basis functions to the image.
shapelet_gresid: This parameter is a Boolean (default is True). If True, then the shapelet decomposition is done on the Gaussian residual image rather that the ch0 image.

Footnotes

process_image: processing an image

General reduction parameters

Adaptive box options

Advanced options

Flagging options

Output options

Multichannel options

À trous wavelet decomposition module

PSF variation module

Spectral index module

Polarization module

Shapelet decomposition module

`process_image`: processing an image