3.1.1. Homados

The main purpose of Homados is to add noise to images. We shall often refer to the original image, as the $\mathcal{D}$ -image, and the noise-added image, as the $\mathcal{S}$ -image. These ‘deep-shallow’ ( $\mathcal{DS}$ ) image pairs can be used to create statistics such as $\mathcal{DS}$ -completeness ( $\mathcal{C_{DS}}$ ) and $\mathcal{DS}$ -reliability ( $\mathcal{R_{DS}}$ ), based on the assumption that the sources detected in the $\mathcal{D}$ -image are real. These statistics are used for real images, where the source inputs are unknown.

An $\mathcal{S}$ -image is created by adding to the $\mathcal{D}$ -image a Gaussian noise map that has been convolved with the corresponding synthesised beam (i.e., BMIN, BMAJ, and BPA). The noise map is created with mean noise, $\mu_{image}$ , and RMS noise, $n\sigma_{image}$ ( $\equiv\sigma_{noise\;\;map}$ ), where n is the desired noise level (i.e., factor), and $\mu_{image}$ and $\sigma_{image}$ are obtained from the $\mathcal{D}$ -image using $\sigma$ -clipping (Akhlaghi & Ichikawa Reference Akhlaghi and Ichikawa2015). This is then convolved with the synthesised beam, from which its RMS noise, $\sigma_{convolved}$ , is computed. For convergence, this process is repeated using the convolved image as input, but with n replaced by $\sigma_{noise\;\;map}\,n/\sigma_{convolved}$ . The final convolved image is then added to the $\mathcal{D}$ -image, obtaining the $\mathcal{S}$ -image.

Fig. 2 shows an example of an $\mathcal{S}$ -image generated by Homados from an Australia Telescope Large Area Survey (ATLAS) Chandra Deep Field South (CDFS) Data Release 1 (DR1) tile (Norris et al. Reference Norris2006). The noise level was scaled by a factor of $n=5$ . This factor is assumed for the $\mathcal{S}$ -image generation in the rest of this paper.

Figure 2. Homados $\mathcal{S}$ -image generation example, using an image cutout sample from an ATLAS CDFS DR1 $2.2^{\circ}\times2.7^{\circ}$ tile (Norris et al. Reference Norris2006). The figures show $\mathcal{D}$ (left) and $\mathcal{S}$ (right) images, zoomed in. The noise level scale factor, n, was set to 5 to generate the shallow image.

In addition, Homados uses $\sigma$ -clipping to compute image statistics, such as, m (median), $\mu$ (mean), $\sigma$ (RMS), $I_{\min}$ (minimum pixel value), and $I_{\max}$ (maximum pixel value). It also does image inversion for FDR calculations.

3.1.2. Cerberus

Cerberus is an extensible interface for running SF modules within the Hydra software suite. It currently supports Aegean, Caesar, ProFound, PyBDSF, and Selavy, as indicated by its command-line interface.Footnote ^e

Table 2. Cerberus RMS and Island parameter definitions in units of $\sigma$ with respect to the background, with soft constraint $\sigma_{Island}<\sigma_{RMS}$ .

$^{\rm a}$ https://github.com/PaulHancock/Aegean/wiki/Simple-usage.

$^{\rm b}$ https://caesar-doc.readthedocs.io/en/latest/usage/app_options.html#input-options.

$^{\rm c}$ https://cran.r-project.org/web/packages/ProFound/ProFound.pdf.

$^{\rm d}$ https://pybdsf.readthedocs.io/en/latest/process_image.html.

$^{\rm e}$ https://www.atnf.csiro.au/computing/software/askapsoft/sdp/docs/current/analysis/selavy.html.

New modules are added through code generation, using Jinja template-codeFootnote ^f in conjunction with DockerFootnote ^g and YAMLFootnote ^h configuration files. The workflow is as follows:

• • Create a containerised SF wrapper:

1. – Create a SF wrapper script

2. – Create a Docker build file wrapper

3. – Update the master docker-compose build file

4. – Build the container image

• • Update the cerberus.py script:

1. – Create a YAML metadata file

2. – Update the master YAML metadata file

3. – Run the Jinja script generator tool

• • Test the Hydra software suite

• • Update the GitFootnote ⁱ repository

  

  Figure 3. Cerberus code generation workflow.

Fig. 3 summarises this high-level workflow: containers for each SF are shown under Container Images and the Docker and YAML configuration files are shown under Configuration Management. The developer must follow a fixed set of rules when adding a new SF, in order for the Jinja template-driven script generator to update Cerberus. (For the purpose of reproducibility, Appendix A provides architectural design notes, using Aegean as an example.) All of this is transparent to the user, who has access to a simple interface, so one does not have to be an expert at using SFs in order to use Hydra.

Hydra’s modular design requires that the user has access to the key elements of a SF’s interface; in particular, access to its ‘RMS-like’ and ‘Island-like’ parameters. In the case of Aegean, for example, this would be seedclip and floodclip (Hancock et al. Reference Hancock, Murphy, Gaensler, Hopkins and Curran2012), respectively. It is important to note that the parameters are not necessarily equivalent between SFs;Footnote ^j regardless, they do affect thresholding and island formation. Consequently, they have the strongest influence on FDR calculations. Table 2 summarises the parameters for the currently supported SFs. These parameters are used by Typhon to baseline the SFs, by minimising their FDRs.

Hydra also requires that SF modules provide optional RMS box and step size parameters, even if they are dummies. Some SF software manuals recommend these parameters be externally optimised, under certain conditions. PyBDSF is an example of such a case.Footnote ^k Regardless, this is also a good way of baselining (i.e., calibrating, Huynh et al. Reference Huynh, Hopkins, Norris, Hancock, Murphy, Jurek and Whitting2012; Riggi et al. Reference Riggi2016) SFs for comparison purposes.

3.1.3. Typhon

Typhon is a tool for optimising the SFs to a standard baseline that can be used for comparison purposes. We have adopted the Percent Real Detections (PRD) metric, as used by Hale et al. (Reference Hale, Robotham, Davies, Jarvis, Driver and Heywood2019) in a comparative study of Aegean, ProFound, and PyBDSF: that is,

(1) \begin{equation} \mbox{PRD} = \frac{N_{image}-N_{inv.\;image}}{N_{image}}\,100 \end{equation}

where $N_{image}$ is the number of detections in the original image, and $N_{inv.\;image}$ is the number of detections in the inverted image.Footnote ^l Basically, if one assumes the image noise is predominately Gaussian, then the peaks detected in the inverted image should statistically match the noise-peaks detected in the non-inverted image. Thus the FDR can be reduced by optimising the PRD. This approach is not suitable for non-Gaussian (or non-symmetric) noise properties, such as the Ricean noise distribution in polarisation images.

Figure 4. Example Typhon PRD of a $2^\circ\times2^\circ$ simulated $\mathcal{D}$ -image (left) and its corresponding $\mathcal{S}$ -image (right). The variation in the PRD with the RMS parameters (in $\sigma$ units) is represented along the horizontal axis, and the variation in the PRD with the island parameters is represented by the error bars. The data points represent average values: that is, Aegean and ProFound indicate the true shape of the curves, due their insensitivity to their island parameters. The SF parameters are listed in Table 2. The dotted horizontal lines indicate the 90% PRD levels.

Typhon uses the RMS and island parameters to optimise the PRD for each SF. Fig. 4 shows Typhon generated PRD curves for a $2^\circ\times2^\circ$ simulated $\mathcal{D}$ -image along with its corresponding $\mathcal{S}$ -image. Typhon identifies the optimal parameters to be those that correspond to the 90% PRD threshold. This threshold is motivated by the desire to use the knee of the PRD curve, whose position appears to be scale-invariant above a certain image size. Although the shape of the curve is not always guaranteed to be smooth, this crude method appears to be quite effective at framing the region of interest around the desired 90% PRD. The 90% to 98% PRD range has been investigated and the former threshold seems to provide reasonable results. Hale et al. (Reference Hale, Robotham, Davies, Jarvis, Driver and Heywood2019) use a 98% PRD to baseline their SFs, beyond which the detection rate degrades. At that cut-off, however, we tend to find a non-scale-invariant increase in the RMS threshold with image size.

Typhon uses the image statistics output from Homados to determine the RMS parameter range over which to optimise the PRD:

(2) \begin{equation} 1.5\sigma\le\mbox{RMS}\le\mbox{RMS}_{\max}\,,\end{equation}

where

\begin{equation*} \mbox{RMS}_{\max} = \left\lceil{\frac{I_{\max}-\mu}{\sigma}}\right\rceil\sigma\,. \nonumber \end{equation*}

with $\mu$ , $\sigma$ , and $I_{\max}$ determined through $\sigma$ -clipping (re. Section 3.1.1). The $1.5\sigma$ lower limit is where the FDR starts to degrade. In general, Typhon samples the PRD from high values to low values in the RMS parameter, while varying the island parameter at each step, until the 90% threshold is reached.

Figure 5. Clustering Algorithm Infographic: The left panel shows the results of 3 hypothetical SFs (red, green, and gold). The middle panel shows the results after clustering, resulting in two clumps, assigned clump_id 1 (upper panel) and clump_id 2 (lower panel). A clump is defined through the spatial overlap between SF detections (i.e., components), in the $\mathcal{D}$ and $\mathcal{S}$ -images together. The components are numbered independently and can be associated with the clump they end up in. For instance, components 1, 2, and 3 are linked together in the $\mathcal{D}$ -image, and, in addition, 1 overlaps with 6, and 3 overlaps with 8, 9, and 10 in the $\mathcal{S}$ -image. Together this set of components populate clump_id 1. Similarly, clump_id 2 is composed of components 4, 5, and 7. Clumps are centred in the Hydra Viewer (Fig. 7), with unassociated components greyed out. The right panel shows the results after the clumps are decomposed into closest (i.e., overlapping centre-to-centre) matches between SFs, in the $\mathcal{D}$ and $\mathcal{S}$ -images, such that there is only one SF with a match in the $\mathcal{D}$ and $\mathcal{S}$ -images. These matched sets are assigned match_ids, with boxes enclosing the extremities of the components. The Hydra Viewer displays these numbers at the centre of the boxes (which are coloured differently here, for emphasis). So clump_id 1 contains match_id 1 = {1, 2, 6} and match_id 2 = {3, 8, 9, 10}, while clump_id 2 contains match_id 3 = {4, 5, 7}.

The island parameters are SF-specific, and are typically defined over a finite range. Table 3 shows the parameter ranges used by Hydra, which are stored in its Configuration Management (Fig. 3). Typhon uses this information along with the constraint $\sigma_{island}<\sigma_{RMS}$ (otherwise, $\sigma_{island}=0.999\,\sigma_{RMS}$ ), as it searches the parameter space.

Table 3. Configured island parameters.

Typhon will also perform an initial RMS box optimisation before optimising the PRD, if it is configured to do so. This is of particular importance for extended objects or around bright sources (Mohan & Rafferty Reference Mohan and Rafferty2015), especially for Gaussian-based extraction SFs such as Aegean, PyBDSF, and Selavy. Typhon uses Aegean’s background/noise image generation tool, bane (Hancock et al. Reference Hancock, Cathryn and Hurley-Walker2018), to search the RMS box size (box_size) and step size (step_size) parameter space,

(3) \begin{equation} \left.\begin{array}{c} \mbox{$\displaystyle3\le\frac{\mbox{box_size}}{[4(\texttt{BMAJ}+\texttt{BMIN})/2]}\le6$}\\[10pt] \mbox{$\displaystyle\frac{1}{4}\le\frac{\mbox{step_size}}{\mbox{box_size}}\le\frac{1}{2}$}\\ \end{array}\right\}\,,\end{equation}

for the lowest background level, $\mu$ (c.f., Riggi et al. Reference Riggi2016). The $4(\texttt{BMAJ}+\texttt{BMIN})/2$ factor represents the bane default box size, where we assume a square box, for simplicity. The limits 3 and 6 are consistent with the rule of thumb that the box size should be 10–20 times larger than the beam size (Riggi et al. Reference Riggi2016). The $1/4$ and $1/2$ bounds are used for providing a smoothly sliding box (Mohan & Rafferty Reference Mohan and Rafferty2015).

The Typhon optimisation algorithm can be summarised as follows.

• • If the RMS box $\mu$ -optimisation is desired:

1. – Minimise $\mu$ over a $6\times3$ box_size by step_size search grid, constrained by Equation (3)

• • Select a centralised $n\times n$ image sample-cutout:

1. – Use an $n^2$ -area rectangle, if non-square image

2. – Use the full area, if image is too small

• • Determine the RMS parameter bounds (Equation (2))

• • For each SF:

1. – If applicable, set the RMS box parameters to the optimised values

2. – Extract the island parameter bounds from Configuration Management (re. Table 3)

3. – Optimise the PRD of the sample-cutout:

• * Iterate the RMS parameter backwards from RMS $_{\max}$

• * At each RMS step, iterate the island parameter, such that, $\sigma_{island}<\sigma_{RMS}$ , otherwise set $\sigma_{island}=0.999\,\sigma_{RMS}$

• * For each RMS-island parameter pair compute the PRD

• * Terminate iterations just before the PRD passes below 90%

• If the PRD is always below 90% choose the highest value.

• • Run the SFs on the full image using their optimised parameters

• • Archive the results in a tarball

For our initial studies, we have chosen to set $n=2.5^\circ$ to provide a sufficiently large region of sky to ensure the finder parameters are not biased by small-scale structure in a given image. Also Aegean, PyBDSF, and Selavy are configured to use the RMS box optimisation inputs from bane, with Selavy only accepting the RMS box size. Appendix B provides details of the SFs and their settings used herein.

For the purposes of placing the SFs on equal footing, we have chosen to restrict Aegean, Caesar and Selavy to single threaded mode, so as to keep the background/noise statistics consistent, at the cost of computational speed. In addition, we keep all of the internal parameters of all of the SFs fixed, instead of tweaking them by hand for each use case. Every effort has been made to keep each SF module as generic as possible.

3.1.4. Hydra

Hydra is the main tool that glues everything together, by running Typhon for $\mathcal{D}$ and $\mathcal{S}$ images, and producing data products such as diagnostics plots and catalogues. One of main underlying features of Hydra is that it uses a clustering algorithm (Boyce Reference Boyce2018) to relate information between SF components in both $\mathcal{D}$ and $\mathcal{S}$ images. In addition, Hydra will also accept simulated catalogue input, under a source-finder pseudonym.

Fig. 5 shows an example of how the clustering algorithm works. All components between all $\mathcal{D}$ and $\mathcal{S}$ SF detections (i.e., catalogue rows) are spatially grouped together, with their overlaps forming clumps with unique clump_ids. The clumps are also decomposed into the closest $\mathcal{DS}$ matches, and assigned unique match_ids. The matches are further broken down by SF into subclumps (not shown), and assigned unique subclump_ids. All of this information is compiled into a cluster catalogue (or table), containing the following key reference elements (columns): cluster catalogue ID, clump ID, match ID, subclump ID, SF $\mathcal{D/S}$ catalogue cross-reference ID, and image depth ( $=\mathcal{D},\,\mathcal{S}$ ). In addition, the catalogue contains common SF output parameters, such as RA, Dec, flux density, etc. There is also a clump catalogue, consisting of rows by unique clump_id, of cluster centroid positions, cutout sizes, total number of components, number of components per SF, SFs with the best residual statistics, etc.

Figure 6. Derivation of the distance metric used for clustering. Here we assume that the space is locally flat, so that $\Delta(C_i,C_j)\approx(\delta^2_{\texttt{RA}_{ij}}+\delta^2_{\texttt{Dec}_{ij}})^{1/2}$ , where $\delta_{\texttt{RA}_{ij}}=(\texttt{RA}_j-\texttt{RA}_i)\cos(\texttt{Dec}_i)$ and $\delta_{\texttt{Dec}_{ij}}=\texttt{Dec}_j-\texttt{Dec}_i$ . The distances from the centres of components $C_i$ and $C_j$ to their edges, along a ray between them, is given by $r_i$ and $r_j$ , respectively: that is, $r_\mu$ , is a standard geometrical expression in terms of angle $\beta_\mu=\pi/2-(\theta_\mu-\eta)$ with respect to the ray and the semi-major axis $a_\mu$ , where $\theta_\mu$ is the position angle, $\eta=-\tan^{-1}(\delta_{\texttt{Dec}_{ij}}/\delta_{\texttt{RA}_{ij}})$ , and $\mu=i,\,j$ . The grey area outside the ellipses is the skirt, whose extent is determined by f.

Figure 7. Hydra Viewer Infographic: The cutout viewing section of Hydra’s local web-viewer interface. At the top is the navigation bar, which allows the user to navigate by clump ID, go to a specific clump, turn on/off cutout annotations, or examine S/N bins of diagnostic plots such as $\mathcal{C}$ and $\mathcal{R}$ , by using the Mode button. The main panel contains $\mathcal{D}$ (top) and $\mathcal{S}$ (bottom) square image cutouts (first column) and SF residual image cutouts (remaining columns), centred about a given clump’s centroid. Here the annotation is turned on, with the neighbouring clumps greyed out. The table at the bottom (not to scale) is the cluster table rows for the clump, with the following columns: cluster catalogue ID, SF catalogue cross-reference ID, clump ID, subclump ID, match ID, SF or $\mathcal{J}$ catalogue name, image depth, RA ( $^\circ$ ), Dec ( $^\circ$ ), semi-major axis ( $^{\prime\prime}$ ), semi-minor axis ( $^{\prime\prime}$ ), position angle ( $^\circ$ ), total flux density (mJy), bane RMS noise (mJy), S/N (total flux over bane RMS noise), peak flux (mJy beam $^{-1}$ ), normalised-residual RMS (mJy (arcmin $^2$ beam) $^{-1}$ ), normalised-residual MADFM (mJy (arcmin $^2$ beam) $^{-1}$ ), and normalised-residual $\Sigma I^2$ ((mJy (arcmin $^2$ beam)) $^{-2}$ ). The normalised-residual statistics are normalised by the cutout area (arcmin $^2$ ). This statistical information is also shown below each cutout, along with the number of components (N), and cutout size (Size, in arcmin). This figure is to illustrate the layout of the Hydra viewer, not the details. It shows screen shots from the Hydra Viewer pasted together, hence the fonts appear small. The data at the bottom is raw output from the cluster table, which is not rounded in this version of the software.

Fig. 6 shows the derivation of the distance metric used in the clustering algorithm. The algorithm uses the following distance metric constraint to determine the overlap between two elliptical components, $C_i$ and $C_j$ , with centre-to-edge distances, $r_i$ and $r_j$ , along an adjoining ray.

(4) \begin{equation} \Delta(C_i,C_j)\le r_i+r_j\end{equation}

where

\begin{equation*} \Delta(C_i,C_j)=\sqrt{(\texttt{RA}_j\!-\!\texttt{RA}_i)^2\cos^2(\texttt{Dec}_i)\!+\!(\texttt{Dec}_j\!-\!\texttt{Dec}_i)^2}\,,\nonumber\end{equation*}

is the distance metric, and

\begin{equation*} r_\mu=\frac{a_\mu b_\mu}{\sqrt{a_\mu^2\cos^2(\theta_\mu-\eta)+b_\mu^2\sin^2(\theta_\mu-\eta)}}\,,\nonumber\end{equation*}

are the centre-to-edge distances, for $\mu=i,j$ , with

\begin{equation*} \eta = -\tan^{-1}\left[\frac{\texttt{Dec}_j-\texttt{Dec}_i}{(\texttt{RA}_j-\texttt{RA}_i)\cos(\texttt{Dec}_i)}\right]\,,\nonumber\end{equation*}

$a_\mu$ is the semi-major axis, $b_\mu$ is the semi-minor axis, and $\theta_\mu$ is the position angle (defined in the same manner as the beam PA, Greisen Reference Greisen2017). So components satisfying this constraint are clumped together.

Hydra also provides a web-viewer (known as the Hydra Viewer) for exploring image and residual image cutouts by clump_id, along with corresponding cluster table information. Fig. 7 provides a detailed description of the Hydra Viewer’s cutout interface. As indicated in the figure, the viewer has radio component annotations that can be toggled on/off. Fig. 8 provides a more detailed example. Table 4 describes the annotation colours, which are stored as metadata in Hydra’s Configuration Management (Fig. 3).

The following is a list of the data products produced by Hydra:

• • Typhon Metrics

1. – $\mathcal{D/S}$ Diagnostic Plots of

1. * PRD

2. * PRD CPU Times

3. * Residual RMS

4. * Residual MADFM

5. * Residual $\Sigma I^2$

1. – Table of $\mathcal{D}$ and $\mathcal{S}$ optimised RMS and island parameters

• • $\mathcal{D/S}$ Catalogues

1. – SF Catalogues

2. – Cluster Catalogue

3. – Clump Catalogue

• • Optional Simulated Input Catalogue

• • $\mathcal{D/S}$ Cutouts

1. – Un/annotated Images

2. – Un/annotated Residual Images

• • $\mathcal{D/S}$ Diagnostic Plots

1. – Clump Size Distributions

2. – Detections vs S/N

3. – $\mathcal{C}$ vs S/N

4. – $\mathcal{R}$ vs S/N

5. – Flux-Ratios: $S_{out}/S_{in}$ vs S/N

6. – False-Positives vs S/N

• • Hydra Viewer: Local Web-browser Tool

All of this information is stored in a tarball. The Hydra Viewer allows the user to view all of these data products. The cutout viewer portion is linked only to the cluster catalogue. It is accessible through an index.html file in the main tar directory.

5.1.1. Simulated compact sources, CMP

The simulated image, shown in Fig. 11, is produced in two steps; generation of a noise image, followed by the addition of artificial sources. We use miriad (Sault, Teuben, & Wright Reference Sault, Teuben, Wright, Shaw, Payne and Hayes1995) to generate the artificial noise image, using the following steps. The imgen task was used to produce a 1800 $\times$ 1800 pixel image, with 4” pixels, (i.e., a $2^{\circ} \times 2^{\circ}$ field) populated by random Gaussian noise of RMS $20\,\unicode{x03BC}$ Jy beam $^{-1}$ . This image was convolved using convol to mimic a 15” FWHM beam, which has the effect of increasing the noise level by a factor of 2.8, so the resulting image is then scaled using maths to divide by this factor, restoring the original noise level of $20\,\unicode{x03BC}$ Jy beam $^{-1}$ .

Figure 11. Simulated map with point-like (compact) sources. The coordinates are arbitrarily set, and the FWHM is set to 15”.

To generate the properties of the artificial sources, we use the 6th order polynomial fit to the 1.4 GHz source counts from Hopkins et al. (Reference Hopkins, Afonso, Chan, Cram, Georgakakis and Mobasher2003), which is consistent with more recent source count determinations for the flux density range considered here (e.g., Gürkan et al. Reference Gürkan2022, and references therein). A sequence of 34 exponentially spaced bins in flux density was defined, ranging from $50\,\unicode{x03BC}$ Jy to $1\,$ Jy, and within each bin the number of sources was calculated from the source count model. The flux density for each artificial source was assigned randomly between the extrema of the bin in which it lies. Source positions were also assigned randomly, with no attempt to mimic the clustering properties of real sources. For the $2^{\circ} \times 2^{\circ}$ field with a flux density limit of $50\,\unicode{x03BC}$ Jy, we end up with a list of 9075 artificial sources.

These sources were added to the noise image using the Python module Astropy (Astropy Collaboration et al. Reference Astropy Collaboration2013, Reference Astropy Collaboration2018) by constructing 2D Gaussian models with the FWHM of the restoring beam (15”) and scaling the amplitude to represent the randomly assigned peak flux density of the source. Given the sources modelled here are assumed to be point-like (compact), the peak flux density for a source has the same amplitude as the integrated flux density. Using this Gaussian model for each source, we generated an image array for each source to be added into the simulated image. We used Astropy again to convert the RA/Dec location of the source to pixel locations and each source was added to the simulated image at the desired location.

5.1.2. Simulated extended sources, EXT

Following a similar procedure as in Section 5.1.1, we generated a sky model of size 1800 $\times$ 1800 pixels ( $4^{\prime\prime}$ pixel size, 2 $^{\circ}\times$ 2 $^{\circ}$ field of view) with both point-like and extended sources injected. The noise level is again set to 20 $\unicode{x03BC}$ Jy beam $^{-1}$ . Extended sources are 2D elliptical Gaussians with randomised position angle and axis ratio, with axis ratio varying between 0.4 to 1. A maximum major axis size was set at three times the synthesised beam size ( $45^{\prime\prime}$ , with a $15^{\prime\prime}$ FWHM beam, as in Section 5.1.1).

A total of 9974 artificial sources were injected, corresponding to a source density of 2500 deg $^{-2}$ , with 10% being extended sources. The peak flux densities S of both point-like and extended sources were set to follow an exponential distribution $10^{-\lambda S}$ with $\lambda$ =1.6, consistent with that seen in early ASKAP observations of the Stellar Continuum Originating from Radio Physics In Our Galaxy (Scorpio, Umana et al. Reference Umana2015) field (Riggi et al. Reference Riggi2021). The minimum peak flux density for all sources was set at $50\,\unicode{x03BC}$ Jy and the maximum fixed at 1 Jy for compact sources and 1 mJy for extended sources. The final simulated image, shown in Fig. 12, was produced by convolving this input sky model using the CASAFootnote ^p imsmooth task and a target resolution of $15^{\prime\prime}$ .

Figure 12. Simulated image with both point-like (compact) and extended sources. The sky coordinates are arbitrarily chosen.

It is important to note here that, unlike the compact source simulation above where the faintest injected source lies at $\mbox{(S/N)}_{\min}\sim2.5$ , in the extended source simulated image 20.6% of the injected sources fall below $\mbox{S/N}=1$ . This is by design, to provide a more realistic image, and to test the impact on the SFs of having real sources lying below the noise level (re. M. M. Boyce et al., in preparation).

Table 6. Hydra $\mu$ -optimised box_size and step_size inputs for Aegean, PyBDSF, and Selavy, $^{\rm a}$ using CMP and EXT $\mathcal{D/S}$ -image data.

$\!\!^{\rm a}$ Selavy only accepts box_size.

5.2.1. Optimisation run results

Hydra uses Typhon to set the RMS box and island parameters for each SF (Aegean, Caesar, ProFound, PyBDSF, and Selavy) to ensure a 90% PRD cutoff. The RMS box parameters, obtained from Typhon’s $\mu$ -optimisation routine, were used by Aegean, PyBDSF, and Selavy. An $\mathcal{S}$ -image was generated by Homados, and RMS box and island parameters similarly estimated. Tables 6 and 7 summarise these results for our CMP and EXT images.

For the CMP source $\mathcal{D}$ -image, $\mu\sim22\,\unicode{x03BC}$ Jy beam $^{-1}$ (Table 6) which is consistent with the design RMS noise of $20\,\unicode{x03BC}$ Jy beam $^{-1}$ (re. Section 5.1). For the EXT source $\mathcal{D}$ -image, $\mu\sim68\,\unicode{x03BC}$ Jy beam $^{-1}$ which is higher due to the inclusion of extended structures with a slightly higher source density, and the fact that the box size has doubled. This is consistent with an average source size increase from $15^{\prime\prime}$ to $30^{\prime\prime}$ . We also note that $\mu_{\mathcal{S}}/\mu_{\mathcal{D}}\sim 5$ for all images, which is consistent with the factor of 5 noise increase for the $\mathcal{S}$ -images (Section 3.1.1). So the RMS box statistics are consistent with what might be expected from the simulated images.

In Table 7 we notice, in the broadest sense, that the RMS and island parameters are similar for the $\mathcal{D/S}$ -images, row-wise. Fig. 13 shows stacked plots of $\mathcal{D/S}$ source counts for the images. In the CMP and EXT $\mathcal{D}$ -images, there is some variability in the source counts, that is, $N\sim4\,650\pm890$ , except for Selavy being consistently lower by $N\sim1\,770\pm950$ . In the $\mathcal{S}$ -image case, the variability is fairly tight, with $N_{\small {\rm CMP}}\sim661\pm54$ and $N_{\small {\rm EXT}}\sim1\,258\pm70$ , except for Selavy which is consistently low with $N_{\small {\rm CMP}}=427$ and $N_{\small {\rm EXT}}=787$ .

Table 7. Typhon run statistics for CMP and EXT images, with SF, image depth, SF RMS parameter ( $n_{rms}$ [ $\sigma$ ]), SF island parameter ( $n_{island}$ [ $\sigma$ ]), source counts (N), residual RMS ( $\unicode{x03BC}$ Jy beam $^{-1}$ ), and residual MADFM ( $\unicode{x03BC}$ Jy beam $^{-1}$ ) columns. $^{\rm a}$

$\!\!^{\rm a}$ The MADFM estimators are normalised by 0.6744888 (Whiting Reference Whiting2012).

Figure 13. SF CMP and EXT $\mathcal{D/S}$ -image detection stacked plots (re. the N columns of Table 7).

For the CMP $\mathcal{D}$ -image, the RMS and MADFM residual statistics are all $\sim$ $19\,\unicode{x03BC}$ Jy beam $^{-1}$ , with the exception of Selavy having a significantly higher RMS ( $\sim$ $45\,\unicode{x03BC}$ Jy beam $^{-1}$ ), likely the primary cause for the reduction in numbers of detected sources it reports. For the EXT $\mathcal{D}$ -image, the RMS values increase in order from ProFound, Caesar, Aegean, PyBDSF, to Selavy, with the latter three being comparable to the former being at extreme end. This is also reflected in their MAFDMs, although here the values are somewhat comparable. For the $\mathcal{S}$ -image case, everything is similar within each image data set, except for Selavy having a higher RMS in the EXT image case. These observed discrepancies for Selavy are likely due to its use of robust statistics in the background estimation (Section B.5). In contrast, the MADFM is similar in both cases for all SFs.

Table 8 shows ratios of deep-to-injected ( $\mathcal{D}\,:\,\mathcal{J}$ ) and shallow-to-deep ( $\mathcal{S}\,:\,\mathcal{D}$ ) source counts: that is, $\mathcal{D}\,:\,\mathcal{J}$ indicates the fraction of sources recovered in the simulated images, whereas $\mathcal{S}\,:\,\mathcal{D}$ indicates the recovery rate assuming the deep detections are ‘real’. Also included are $\mathcal{S}\,:\,\mathcal{J}$ recovery rates, for comparison.Footnote ^q The $\mathcal{D}\,:\,\mathcal{J}$ recovery rates are not expected to reach 100%, as the simulations includes low S/N sources, and for EXT sources, some lying below S/N $=1$ . The $\mathcal{S}\,:\,\mathcal{D}$ recovery rate is lower for CMP sources than EXT, which both track consistently with $\mathcal{D}\,:\,\mathcal{J}$ .Footnote ^r This may be due to some confusion in the EXT image (Fig. 12), given the $\mathcal{S}\,:\,\mathcal{J}$ recovery rates are similar.

Table 8. Ratios of deep-to-injected ( $\mathcal{D}\,:\,\mathcal{J}$ ) shallow-to-injected ( $\mathcal{S}\,:\,\mathcal{J}$ ), and shallow-to-deep ( $\mathcal{S}\,:\,\mathcal{D}$ ) sources. The $\mathcal{D}/\mathcal{S}$ source counts (N) are provided in Table 7, and the injected source counts are 9075 and 9974 for CMP and EXT sources, respectively (re. Section 5.1).

The consistency of the RMS box, MADFM, $\mathcal{D}\,:\,\mathcal{J}$ , $\mathcal{S}\,:\,\mathcal{J}$ , and $\mathcal{S}\,:\,\mathcal{D}$ statistics provides a good indication that Hydra’s optimisation routines are performing robustly.

5.2.2. Source size distributions

Fig. 14 shows the major-axis distribution for our simulated image data. Both the $\mathcal{D}$ and $\mathcal{S}$ source detections are combined, as the $\mathcal{S}$ provides no contrasting information and its statistics are quite low (Fig. 13). Note that the size estimates for different SFs use different methods and are not necessarily directly comparable. Those SFs that fit Gaussians to source components report size as a FWHM, while others (such as ProFound) use different measures, such as a flux-weighted fit (see Appendix B for details).

Figure 14. Major-axis distributions for (a) CMP and (b) EXT sources (with $\mathcal{D}$ and $\mathcal{S}$ both included). The vertical dashed-line represents the beam size. The distributions of the injected sources are shown in grey. Recall that size estimates between SFs are not necessarily directly comparable as they are estimated using different methods.

For the CMP source case, ProFound and PyBDSF tend to overestimate the source sizes. This is likely due to deblending issues or the incorporation of noise spikes. As for EXT sources, PyBDSF tends to most accurately recover the extended source sizes, but, along with Caesar also has the most outliers. These could be attributed to size overestimates due to inclusion of noise spikes or adjacent sources in the fitted sizes. All other SFs tend to marginally underestimate the EXT source sizes. Components smaller than $15^{\prime\prime}$ are attributed primarily to noise spikes, but also occasionally to underestimating the source sizes.

5.3.1. Simulated sources

Fig. 15 shows $\mathcal{C_{D}}$ vs S/N (top) and $\mathcal{R_{D}}$ vs S/N (bottom) for CMP (left) and EXT (right) $\mathcal{D}$ -images, where the signal-to-noise (S/N) is the ratio of the $\mathcal{D}$ -signal to $\mathcal{D}$ -noise. Fig. 16 shows the corresponding results for $\mathcal{S}$ -images, where S/N is the ratio of the $\mathcal{S}$ -signal to $\mathcal{S}$ -noise.

Figure 15. Simulated $\mathcal{D}$ -image compact (left) and extended (right) source $\mathcal{C_{D}}$ (top) and $\mathcal{R_{D}}$ (bottom) vs S/N ( $\mathcal{D}$ -signal/Deep-noise) plots. The $\mathcal{D}$ -noise is computed using bane.

Figure 16. Simulated $\mathcal{S}$ -image compact (left) and extended (right) source $\mathcal{C_{S}}$ (top) and $\mathcal{R_{S}}$ (bottom) vs S/N ( $\mathcal{D}$ -signal/ $\mathcal{S}$ -noise) plots. The $\mathcal{S}$ -noise is computed using bane.

For CMP sources all SFs show a similar behaviour in general. The completeness metric starts to decline slowly from 100% at S/N $\lesssim 30$ , and drops rapidly toward zero below S/N $\sim 5$ . False detections are typically limited to about 10% of the sample down to S/N $\sim 5$ , but can appear in some SFs as high as S/N $\sim 30$ . In general terms, Aegean provides the highest level of completeness at any given S/N, with a well-behaved decline in reliability below S/N $\sim 5-6$ . At the other end of the scale, Selavy has the lowest level of completeness at any S/N, but the best reliability (fewest false detections).

All SFs seem to miss some bright sources, with Selavy standing out as the poorest in this regard. This is likely due in part to the handling of overlapping sources. Both PyBDSF and ProFound report the largest numbers of false sources (seen in $\mathcal{R_D}$ and $\mathcal{R_S}$ ) at high S/N. For PyBDSF this is a consequence of overestimating source sizes (Fig. 14), especially in the presence of closely neighbouring sources, or nearby noise spikes, by quite significant amounts in some cases (see Paper II). For ProFound this arises due to the blending of neighbouring sources (see Paper II).

For EXT sources there is generally poorer performance overall compared to those for the CMP sources, most clearly seen in the $\mathcal{D}$ -image results (Fig. 15). There are also several artefacts appearing at unphysically low S/N ( $\mbox{S/N}<1$ ) arising from spurious faint detections. Even at reasonable S/N ( $10<\mbox{S/N}<100$ ) there is measurable incompleteness, and reliability that dips as low as 80% for some SFs. Here Aegean appears to perform the best, with Selavy showing much poorer performance.

5.3.2. Metrics for real sources

Fig. 17(a) and (b) show $\mathcal{C_{DS}}$ vs S/N and $\mathcal{R_{DS}}$ vs S/N, respectively, for CMP sources. As we know what the true sources are, we can explore the validity of $\mathcal{C_{DS}}$ and $\mathcal{R_{DS}}$ by removing any false detections in the $\mathcal{D}$ and $\mathcal{S}$ -images, that is, we can compute $\mathcal{\tilde{C}_{DS}}$ (Equation (11)) and $\mathcal{\tilde{R}_{DS}}$ (Equation (12)). The results are shown in Fig. 17(c) and (d), respectively.

Figure 17. $\mathcal{C_{DS}}$ (a and g), $\mathcal{R_{DS}}$ (b and h), $\mathcal{\tilde{C}_{DS}}$ (c and i), $\mathcal{\tilde{R}_{DS}}$ (d and j), $\delta\mathcal{C_{DS}}$ (e and k), and $\delta\mathcal{R_{DS}}$ (f and l) vs S/N for CMP (top set) and simulated-extend (bottom set) sources. The S/N are expressed as $\mathcal{D}$ -signal/ $\mathcal{S}$ -noise and $\mathcal{S}$ -signal/ $\mathcal{S}$ -noise for completeness and reliability, respectively, where the $\mathcal{S}$ -noise is computed using bane.

Comparing $\mathcal{\tilde{C}_{DS}}$ and $\mathcal{\tilde{R}_{DS}}$ with $\mathcal{C}_{\mathcal{DS}}$ and $\mathcal{R}_{\mathcal{DS}}$ , respectively, $\tilde{\mathcal{C}}_{\mathcal{DS}}$ is largely unchanged, while $\tilde{\mathcal{R}}_{\mathcal{DS}}$ does not show the dip in $\mathcal{R}_{\mathcal{DS}}$ around S/N $\sim$ 10 seen for most SFs. It also excludes the decline seen by Aegean at low S/N. This suggests that the apparently poorer estimated reliability in $\mathcal{R}_{\mathcal{DS}}$ arises from the existence of spurious sources detected in the $\mathcal{D}$ -image that are (not surprisingly) missed in the $\mathcal{S}$ -image.

To quantify these results we define residuals between our unconstrained ( $\mathcal{C}_{\mathcal{DS}}$ and $\mathcal{R}_{\mathcal{DS}}$ ) and constrained ( $\mathcal{\tilde{C}_{DS}}$ and $\mathcal{\tilde{R}_{DS}}$ ) metrics: that is, the residual completeness,

(13) \begin{equation} \delta\mathcal{C_{DS}}=\tilde{\mathcal{C}}_{\mathcal{DS}}-\mathcal{C}_{\mathcal{DS}}\,, \end{equation}

and residual reliability,

(14) \begin{equation} \delta\mathcal{R_{DS}}=\tilde{\mathcal{R}}_{\mathcal{DS}}-\mathcal{R}_{\mathcal{DS}}\,. \end{equation}

These quantities are shown in Fig. 17(e) and (f), respectively. In general, they are expected to be positive, as there should be an excess of false detections in $\mathcal{C_{DS}}$ and $\mathcal{R_{DS}}$ compared to $\mathcal{\tilde{C}_{DS}}$ and $\mathcal{\tilde{R}_{DS}}$ , by construction. Negative values may appear when real input sources are detected but poorly fit, leading to inconsistent flux densities between the $\mathcal{D}$ and $\mathcal{S}$ -images. The figure shows $\delta\mathcal{C_{DS}}$ predominantly highlighting a small number of spurious PyBDSF detections, and $\delta\mathcal{R_{DS}}$ emphasising the distribution in S/N of the false detections.

Recognising these limitations, while also noting that for several SFs $\delta\mathcal{C_{DS}}$ and $\delta\mathcal{R_{DS}}$ are small, the approach of estimating completeness and reliability for a given finder based on real images is not unreasonable. Clearly it is not as robust as doing so using known injected sources as a reference, but it may be a useful addition to analyses comparing finder performance on real data containing imaging artefacts and other hard to simulate systematics. Similar conclusions can be drawn for the EXT source case (Fig. 17(g) through (l)).

Finally, we note that $\mathcal{C_{DS}}$ and $\mathcal{R_{DS}}$ in Fig. 17(a) and (b), respectively, for CMP sources fare much better than their EXT counterparts, (g) and (h), respectively. The metrics in the latter case perform even more poorly than their $\mathcal{C_{D}}$ , $\mathcal{R_{D}}$ , $\mathcal{C_{S}}$ , and $\mathcal{R_{S}}$ counterparts (Figs. 15 and 16). This could perhaps be due to confusion, given the source density (Fig. 12) and high MADFM statistics (Table 7).