Efficient Source Finding for Radio Interferometric Images

TL;DR

This paper introduces new methods to accelerate and automate source finding in radio interferometric images, addressing the challenges posed by upcoming large-scale observatories like the SKA.

Contribution

It presents novel techniques for efficient and automated source detection tailored for fast imaging systems in radio astronomy.

Findings

Improved speed of source detection algorithms.

Enhanced automation in the source finding process.

Potential applicability to large-scale radio telescopes like SKA.

Abstract

Object detection in astronomical images, generically referred to as source finding, is often performed before the object characterisation stage in astrophysical processing work flows. In radio astronomy, source finding has historically been performed by bespoke off-line systems; however, modern data acquisition systems as well as those proposed for upcoming observatories such as the Square Kilometre Array (SKA), will make this approach unfeasible. One area where a change of approach is particularly necessary is in the design of fast imaging systems for transient studies. This paper presents a number of advances in accelerating and automating the source finding in such systems.

Figures

Tables

Table 1. Table 1 : Main functions of the labeling procedure and the required number of image passes for two implementation approaches: straightforward and optimised. An image pass implies writing operations on the image unless otherwise stated, e.g “only read”.

	No of image passes
Labelling functions	Straightforward approach	Optimised approach
Positive thresholding	1	1
Scanning phase - positive sources	1
Negative thresholding	1
Scanning phase - negative sources	1
Labelling phase - positive sources	1	1
Detection of valid positive sources	1 (only read)
Labelling phase - negative sources	1
Detection of valid negative sources	1 (only read)
Removal of invalid positive sources	1	1
Removal of invalid negative sources	1
Combination of label maps	1
Computation of source barycentre	1 (only read)
Total number of image passes	9 + 3 (only read)	3

Equations

$$\sigma^{2}=\Bigg{(}\frac{1}{n}\underbrace{\sum_{i=1}^{n}x_{i}^{2}}_{\rm Acc2}\Bigg{)}-\Bigg{(}\frac{1}{n}\underbrace{\sum_{i=1}^{n}x_{i}}_{\rm Acc1}\Bigg{)}^{2}$$

$$S_n=S_{n-1}+(x_n-{\mu }_{n-1})(x_n-{\mu }_n),$$

$$S_n=n{\sigma }_n^2.$$

$${\mu }_n={\mu }_{n-1}+\frac{1}{n}(x_n-{\mu }_{n-1}),$$

Full text

††journal: Astronomy and Computing

1 Introduction

Object detection, known as source finding in radio astronomy, is a component of the standard analysis pipeline for most astronomical surveys. Source finding is used to identify individual astrophysical objects in survey data, while object characterisation estimates their properties such as location, intensity and shape. Once extracted, these properties are then summarised to form radio survey catalogues. These catalogues are constructed with the intention of providing a compressed version of the information contained in the survey and hence of the astrophysical systems represented in those data.

The software packages used for this process are known as source finders and typically provide object detection and characterisation capabilities to produce the catalogues of detected objects as their output. Historically, due to the faintness of the radio sky and the comparative sparsity of extragalactic objects in radio surveys, the development of these source finders has focused largely on statistical completeness (contain all sources present in the image) and reliability (all sources found and extracted are real; [1]), rather than optimization of the algorithms for run time and memory usage. Consequently, the source finding algorithms developed for radio astronomy over the last few decades have been designed to be highly reliable, returning lists of sources that are expected to be complete and include only a small or zero fraction of false positives.

For the new generation of radio telescopes such as MeerKAT [2], ASKAP [3] and ultimately the Square Kilometre Array (SKA; [4]) the improvements in the sensitivity of radio telescopes, coupled with increased survey speeds due to larger instantaneous fields-of-view (FOV), radio astronomy is starting to produce large area, sensitive maps of the sky at high cadence. As this trend increases, the data processing for these facilities will need to be fully automated, with minimal manual input, and be able to operate reliably and quickly on large volumes of data. Consequently, source finding algorithms will not only need to be accurate but will also be required to work efficiently on much larger image datasets as a matter of routine.

One specific case where this will be required is for the slow transients pipeline (STP) of the SKA telescope. This imaging pipeline is designed to identify transient astronomical sources, the brightness of which varies on timescales long enough for them to be detected in individual residual snap-shot images from the telescope, where a residual snapshot image is defined as an image that can be created from $\sim 1$ second of data subtracted by the sky model data. The STP is expected to produce images at this cadence and identify any transient sources in those residual images. However, the FOV (in terms of number of pixels) for the SKA will make it impossible to store each individual residual snapshot image and consequently the processing for such a pipeline must be performed in soft real time.

When searching for such short-lived radio transients it is not possible to re-observe in order to recover objects missed by initial processing, hence accurate source finding is a key component of this processing. Consequently, to perform the processing for the SKA STP requires not only the imaging of the data to be performed efficiently, but also the source finding in each residual snap-shot image to be optimized for speed.

In this paper we describe a new algorithm for computationally efficient source finding in interferometric radio images: ceres. This algorithm has been designed as part of the development for the SKA slow transients pipeline with its codebase being available on GitHub [5]. The source finding is a component of the STP workflow being applied to the residual dirty image generated previously in the imager stage specifically by computing the residual visibility data (i.e. by subtracting the known sky model visibility data), gridding the residual visibility data and computing the FFT of the gridded data. In the SKA STP, this procedure should be performed at $\sim 1$ second cadence.

An important feature of the algorithm developed is the ability to find both positive and negative sources in the residual dirty image. The potential to search for negative sources is important for imaging applications in radio astronomy that monitor fading transient sources, as well as those that use polarization data, since polarized sources can appear as negative features in Stokes images. While the algorithm discussed performs both object detection and characterisation, it is important to recall that this paper is focused on the detection stage, for which significant speedup techniques are proposed.

The paper proceeds as follows: in § 2 we give an overview of current source finding algorithms and the computational steps required in the source finding process, in § 3 we describe the changes and optimizations in each step that have been implemented in this new algorithm, and in § 4 we present benchmarking of the ceres algorithm in comparison with other widely used source finders. In § 6 we draw our conclusions.

2 Source Finding in Radio Interferometry

2.1 Image-plane vs uv-plane source analysis

Source finding is the process of analysing data to separate signal from noise and infer the flux intensity, position, and (in the case of extended sources) structure of radiation-emitting sources on-sky. For optical data this analysis is performed on pixelated images, and essentially consists of locating clusters of pixels which are outside the expected range of values recorded due to noise alone, then processing each cluster to estimate the required source parameters. Knowledge of the system characteristics may be used to enhance the analysis (e.g. noise-suppression with matched-filters, deconvolution using knowledge of the point-spread function to better parameterise adjacent sources, etc).

Source finding in radio-synthesis maps has evolved in parallel to optical-data processing techniques, but has a distinct set of challenges and requirements. Firstly, the transform of radio-data from recorded antenna-voltages through correlation and Fourier-transform to the image-map produces images with a strong inter-pixel noise correlation on the scale of the synthesized beam. Second, incomplete coverage of the uv-plane may produce a synthesized beam with a complex, large-scale, multimodal structure, where the secondary peaks are known as ‘side-lobes’. As a result, a bright source may obscure or distort the brightness-profile of a a fainter one even if their sky-positions are separated by many multiples of the synthesized-beam width.

The latter problem is often addressed through deconvolution, typically using a variant of the clean algorithm [6]. This ameliorates the issue of side-lobes, at the cost of adding an additional, iterative, non-linear processing step to the data, which may introduce additional degeneracies and can result in significantly different images depending on the clean-parameters. Source finding algorithms are then typically run against the cleaned image.

In some circumstances - particularly if the uv-plane coverage is very sparse, or if the source has excellent signal-to-noise ratio and structure of interest on small scales, it may be more effective to model the sources at the level of the visibility data directly. In [7], the authors discuss a procedure to ‘degrid’ a sky-model to the visibility plane in order fit model parameters through least-squares residuals in the visibilities. However, this requires a prior ‘guess’ at a suitable model of the source structure, and only provides best fit parameters with limited information about the suitability of the model. [8] takes the concept further and use a Monte-Carlo process to perform a full Bayesian inference on both the the source parameters and the systematic noise effects, but are still limited to simple models (single point-source, two point-sources, or an extended Gaussian). Consequently traditional sourcefinding is still the only option for most blind-field surveys, and we limit the rest of this discussion to standard image-plane techniques.

2.2 Image-plane source finding algorithms

There are several widely-used source finding software packages in astronomy which take a largely similar approach. Here we give an overview of the steps involved, and briefly cover differences between the following selection: Bertin’s source-extractor (‘S-Extractor’) [9], the LOFAR-TKP Python Source-Extractor (‘PySE’) [10, 11], and Hancock’s ‘Aegean’ package [12].

2.2.1 Background and noise estimation

The first step in the source finding algorithms is to characterise the background and noise levels of an image, in order to choose a sensible threshold for determining pixels of interest. This step can be remarkably compute intensive due to the need to process all the image-pixels and ideally discard pixels containing significant source-flux, as discussed in §3. Generally, we expect the background pixel distribution to be well approximated by a normal distribution, so this is equivalent to estimating the mode and root mean square (RMS) of the background-pixel sample. Usually some form of sigma-clipping is applied to reject outlier pixels containing source flux. The RMS and mode can then be estimated from the remaining pixels, with either the mean, median, or some combination of both being used to estimate the mode [see e.g. 9]. For radio-synthesis maps the situation is complicated by the strong inter-pixel noise correlation, and so the routine may be modified to only use a sparse sampling of the pixels for RMS-estimation, as in Aegean, or to apply a correction factor based on knowledge of the synthesized-beam parameters, as in PySE [see 11, §4.3.1 therein].

2.2.2 Connected-component labeling

Once the background and noise-variance are known we can pick a suitable threshold-level to reject most background-noise pixels. Typical values are a factor 5–10 times the RMS above the background level (referred to as a 5–10$\sigma$ threshold in the following text). Alternatively a threshold may be chosen by specifying a desired false-detection rate [13, 14]. Using the threshold of choice, the pixels can be converted to a binary-map (above/below threshold), which is then subdivided into labelled regions representing sources or clusters of sources. The above-threshold regions are often referred to as ‘islands’, by way of analogy to islands of land sticking out of the ocean. This binary-map labelling procedure is more generally known as ‘connected components labelling’, and has a fairly extensive treatment in the computer-science literature. S-Extractor uses a custom implementation of Lutz’s algorithm [15], while PySE and Aegean make use of the SciPy function scipy.ndimage.measurements.label [16]. One issue with this approach as described — using a single threshold — is that an island may quite frequently consist of a single pixel above the threshold level, for a source close to the detection limit. This makes it difficult to estimate any parameters other than the value and location of the island’s peak pixel. PySE addresses this edge-case by making use of two thresholds referred to as the ‘detection and analysis’. The ‘detection’ threshold is chosen at the same 5–10$\sigma$ level as before, while the ‘analysis’ threshold is chosen at a lower level that would result in many more false detections if used alone, say 3–4$\sigma$. The thresholding and connected-component labelling procedure is performed using the analysis threshold, but any islands with a peak value less than the detection threshold are rejected as marginal detections likely due to noise fluctuations.

2.2.3 Deblending

One failing of the ‘island-finding’ approach to source-extraction is that multiple bright sources located in close proximity may result in an above-threshold connected region that encompasses several sources. To circumvent this issue, all three of the sourcefinders from the literature implement some form of user-configurable ‘deblending’ step, wherein every island is (optionally) further analyzed to determine if it can be plausibly broken up into multiple ‘sub-islands’. S-Extractor and PySE do this by allocating a number of thresholds on a logarithmic spacing between the original threshold and the island peak, and then re-running a connectivity analysis on the island using each of the higher thresholds in turn. If multiple sub-islands are detected at this new, higher threshold then further criteria are applied (contained flux as a portion of the original island-flux and peak-value in the sub-island) before decided whether to accept the proposed island-split. The process then repeats at the next (higher) threshold, sub-dividing further as appropriate, for a user-defined number of iterations [9],[10, pp 40–41]. Aegean takes a novel alternative approach, computing a map of local curvature across the image and the designating connected-regions of negative curvature (convex areas) as ‘summits’. These ‘summits’ are then used in conjunction with the standard threshold-map to designate sub-island regions for separate Gaussian component fitting [1, §7 therein].

Since the work described in this paper is primarily focused on high-speed transient-detection in residual images we do not expect multiple close sources, and so have not implemented a deblending step, but the code could be extended in a straightforward manner to apply either deblending method if required for regular blind-field surveys.

2.2.4 Island characterization and fitting

Once island regions have been identified we can attempt to extract summary metrics about each island (hopefully corresponding to a single source). At the simplest level, we can use the position and value of the brightest pixel within an island, but this is a rather crude characterisation. The approach adopted by by S-Extractor, and refined in PySE, is to use the first and second statistical moments of the island pixels, which in the case of a Gaussian profile correspond to the mean and variance of the Gaussian (we can also estimate the rotation angle through some reparameterisation formulae). It is worth noting that sources with flux-level close to the cut-off threshold (and therefore with small islands of designated pixels) will result in underestimates of the variance (similar to a truncated Gaussian); [10] derives correction factors for this systematic bias and these correction factors are applied in PySE. S-Extractor then returns these moment-estimate parameters and halts. However, it is often possible to improve accuracy through least-squares fitting of a Gaussian profile, particularly for compact sources at high signal-to-noise ratio [see e.g. 10, figures 2.4–2.5]. PySE implements Gaussian fitting as a secondary step, using the moments-parameter estimates as initial parameters for the fit and applying scipy.optimize.leastsq, a wrapper around the MINPACK implementation of lmdif (a Levenberg-Marquardt optimisation algorithm). Aegean also implements Gaussian fitting, employing some sensible constraints on the Gaussian parameter values, while simply setting the initial parameters using the peak pixel position and flux, with initial beam parameters set to the synthesized beam profile. Aegean makes use of another Python library, ‘lmfit’, which extends the scipy.optimize.leastsq routine to provide fitting under parameter constraints.

3 Efficient Source Finding

Source finding is the final stage of an imaging workflow. For sparse residual images, such as those proposed for the SKA slow transients pipeline (STP), or data with high filling factors in the uv-plane, this detection step may be performed directly over the dirty image; the only direct data dependency between the previous work flow elements and the source finding stage.

The relevant steps of the source finding algorithm addressed in this work are represented in the block diagram of Figure 1. The algorithm as presented in this paper can be used to find both positive and negative sources in the residual dirty image, although the search for negative sources can be disabled. Searching for negative sources is important in STP since the polarized sources can appear as negative features in Stokes images.

The main bottleneck of a source finding algorithm is generally considered to be the RMS estimation, mainly due to the sigma clipping step performed internally, which requires several passes over the image data. However, other components of source finding, such as the background level estimation and the connected component labelling also present a noticeable computational complexity which can be ignored.

To minimise the run time, we here propose a number of optimisations and improvements to the standard source finding process. Specifically we focus on reducing the number of passes over image data, avoiding image copies and redundant computations, and taking advantage of parallel processing.

In the following sections, we describe each step of the source finding algorithm described pictorially in Fig. 1 and we present the proposed solutions for an efficient implementation using the C++ language.

3.1 Background level estimation

The first step in source finding is to estimate the background level of the image by computing the median of the image. The fastest existing algorithm to compute the median is quickselect [17]. It consists of a selection algorithm that finds the $k^{th}$ smallest element from an array of $n$ elements. When $n$ is odd the median is the $k^{th}$ smallest element with $k=(n+1)/2$, and when $n$ is even it is the mean of the elements $k=n/2$ and $k=n/2+1$. The quickselect algorithm is derived from the well known sorting algorithm quicksort. In brief, quickselect chooses a pivot element $p$ from the array and rearranges the array so that all elements smaller than $p$ are to its left and all elements larger than $p$ are to its right. It then acts recursively on one of the sub-arrays (left or right), depending on where the median element lies.

The quickselect algorithm can perform efficiently in practice and when the pivot element is selected randomly from the array, it presents $\mathcal{O}(n)$ average computational complexity. The main disadvantage of quickselect is the use of in-place partitioning, which rearranges the input array. In applications like imaging, where the input array is required to keep its original order, it is necessary to make a scratch copy of the array.

Typically, in C++ programs, the median function is implemented using the Nth_element function provided in C++ standard library for $k^{th}$ smallest element selection based on quickselect algorithm. The parallel implementation of the C++ standard library [18] includes a parallel version of the Nth_element function that can be used for faster median computation.

An alternative approach for median computation is based on the work of [19], where two methods for fast median computation are described based on the concept of successive binning. The first method, referred to as binmedian, presents $\mathcal{O}(n)$ average computational complexity, similarly to quickselect. Its main advantage relative to quickselect is the significantly faster update of the median as more data are added to the array. However, in the context of source finding application, this is feature not relevant. The second method, referred to as binapprox, computes an approximate median value. It tends to perform faster than binmedian, presenting a worst-case complexity of $\mathcal{O}(n)$.

In this work, both the binapprox method and a simplified version of binmedian were implemented and evaluated. The principle behind these methods is based on the following lemma (proof in [19]):

Lemma 1

If X is a random variable having mean $\mu$, variance $\sigma^{2}$, and median $m$, then $m\in[\mu-\sigma,\mu+\sigma]$

In this context, the basic algorithm of binapprox is given by:

Compute the mean $\mu$ and standard deviation $\sigma$; 2. 2.

Form $B$ bins across $[\mu-\sigma,\mu+\sigma]$ and map each data point to a bin; 3. 3.

Find the bin $b$ that contains the median. This consist in summing the size of each bin until the result goes by the middle value, i.e. the total number of samples divided by 2; 4. 4.

Return the midpoint of bin $b$.

The estimated median can differ from the exact median value by at most half the width of the interval, or $\sigma/B$. Typically, the number of bins is $B=1000$, which results in a maximum error of $1/1000^{th}$ of a standard deviation. The main advantage of binapprox is not only its reduced running times, but also the fact that the input data are not modified, contrary to quickselect. Conversely, its main disadvantage is that the computed median may not correspond to the true median. However, given that the maximum error is $1/1000^{th}$ of the image noise it should be a quite small and irrelevant error.

The first three steps of the binmedian method are similar to those of binapprox; however, instead of returning the midpoint of bin $b$ in the fourth step, binmedian recurses on the set of samples mapped to $b$. This recursion is applied until the number of remaining samples in median bin $b$ is smaller than a fixed constant. The quickselect method is then used to find the median based only on the few remaining samples. More details of this approach are given in [19].

Although it presents $\mathcal{O}(n)$ average computational complexity, the binmedian method requires additional memory access and copy operations that reduce its performance. At each iteration, the full image array needs to be accessed to perform binning of the selected bin samples spread along the array. To avoid multiple accesses to the full array, the samples of the selected bin can be copied into a smaller array. However, this solution implies data copies at each iteration of the algorithm.

In this work it was observed that stopping successive binning after the first iteration and using quickselect over the remaining samples of the selected bin performed more efficiently most of the times. This solution can be viewed as an extension to the binapprox method. Instead of returning the midpoint value of the selected bin, all the samples of the selected bin are copied to a new array and inputted to the quickselect function that returns the median value. Hereafter, we will use the term binmedian to refer to this modified solution for computing exact median. The main advantage of this method, relative to the use of quickselect on the original array is that the binning step significantly reduces the number of samples to be handled by quickselect.

A common characteristic of the binapprox and binmedian methods is the calculation of the mean and standard deviation quantities in the first step of their algorithms. Since the standard deviation is also required for the sigma clipping step performed in the context of RMS estimation, this quantity can be reused thus avoiding redundant computations in the source finding algorithm. Consequently, in addition to the reduced running time for median computation when compared with quickselect, the binapprox and binmedian functions proposed here allow for a reduction of the computational complexity in the RMS estimation step, by reusing the precomputed standard deviation quantity.

A further advantage of the binning methods is that they can be easily implemented using multi-core processing. We used the Intel Threading Building Blocks (TBB) [20] library to perform the binning step in parallel, specifically by assigning distinct threads to separated regions of the image and then combining the per-thread binning results. A parallel approach was also used to estimate the mean and standard deviation of the image specifically by computing the auxiliary quantities (the accumulation values) for distinct regions of the image in parallel and then combining the results, as discussed in § 3.2. To minimise the number of passes over the image, the mean and standard deviation are calculated jointly using a single reading of the image.

For source finding execution, the median function can be selected from the input configuration file. While the binapprox method does not compute exact median, having a negligible error, the alternative binmedian function provides the exact median at the cost of a higher computational complexity.

3.2 RMS estimation

Unless an RMS value is provided a priori as input to the source finding function, estimation of the RMS is performed globally as part of the source finding process. For RMS estimation, sigma clipping is first applied to remove the outlier samples. The sigma clipping function requires computing the standard deviation statistical quantity and subtracting the background level (median) from the image. Since the background level was previously estimated, its value is directly used by the sigma clipping function. Furthermore, as described in § 3.1, the binning-based approach used for background level estimation also computed the standard deviation quantity in its first step. This value is thus kept in memory and passed to the RMS estimation function for the sigma clipping function.

Sigma clipping is the main bottleneck of the RMS estimation procedure. This function iterates over all data, rejecting samples that are discrepant by more than a specified number of standard deviations from a central value, the median. By default, we use 5 iterations for sigma clipping, each resulting in one pass through the data. At each iteration the standard deviation is recomputed, after rejecting the discrepant samples from the calculations.

Here we propose an efficient solution for sigma clipping and RMS estimation that reduces both the number of computations and memory writes. Since the input data image should not be modified, the clipped samples at each iteration of sigma clipping cannot be removed or marked in the image, for instance by using the NaN (Not-a-Number) symbol. A copy of the image could be used for this purpose but, due to its large size and the required number of writes, this solution is not recommended. Instead, the proposed solution compares all image samples, including discrepant samples detected in previous iterations, against the current and previous threshold values for sigma clipping. Although this approach requires additional comparisons in order to know which samples were clipped at previous iterations, it is less costly than creating a scratch copy of the image and removing or marking the clipped samples at each iteration.

Given that sigma clipping recomputes the standard deviation at each iteration using data that are quite similar, with only slight differences due to the clipped samples, we here consider a computationally efficient solution that minimises the number of computations. This proposed method uses auxiliary accumulation variables that are updated each time a sample is clipped.

For an array of $n$ samples, the standard deviation ($\sigma$) can be expressed as:

[TABLE]

As can be seen from this equation, the standard deviation is computed using two auxiliary accumulation values, specifically Acc1 and Acc2. Considering this fact, we may subtract the samples clipped at each iteration from the accumulation variables to recompute the standard deviation quantity. This solution minimises the number of computations required to update the standard deviation quantity, since it does not require one to sum all the point values and their squares repeatedly at each iteration. Furthermore, the computation of these auxiliary accumulation values can be easily performed in parallel, specifically by having multiple threads summing separate regions of the image into local auxiliary accumulation values that are combined in the end. For the first iteration of sigma clipping, these accumulation variables are derived from the mean and standard deviation quantities previously computed during background level estimation.

As an additional optimisation, the sigma clipping and RMS estimation functions are merged. Since the RMS corresponds to the standard deviation of the clipped data, we use the standard deviation value computed in the last iteration of sigma clipping function as the estimated RMS value.

A disadvantage of this method is that it relies on a standard deviation equation that is not considered numerically stable. This equation can suffer from loss of precision due to the difference in magnitude between one sample and the sum of all samples in the accumulation variables. A numerically stable equation for the standard deviation is:

[TABLE]

assuming

[TABLE]

and with $\mu$ being the mean quantity that can be computed using the following numerically stable equation:

[TABLE]

However the drawback of using these equations is their high impact on the computational performance. It can be observed that these expressions involve a larger number of operations, in particular multiplication and division operations, that would result in a longer running time. In the interest of maintaining the performance we do not currently implement this more stable method.

In summary, we use an efficient implementation of the RMS estimation function, which only requires 5 passes through the image, one per iteration of the sigma clipping function. For a faster performance, we use multi-threaded processing in which different threads compute the auxiliary accumulation values for separated regions of the image and then combine the results. The number of performed computations is minimised by making use of the previously computed standard deviation and median quantities determined for background level estimation and adopting an efficient solution to update their values when some samples are clipped. Furthermore, the algorithm is designed so that the input image matrix does not need to be copied.

3.3 Source Labeling

Based on the estimated RMS value, the analysis and detection threshold values are determined. The image is then thresholded using the analysis threshold and the connected components (the sources) are labeled. Finally, the detection threshold is used to identify valid sources, following the same method as described for PySE in § 2.2.2.

A Connected Component Labeling [21] algorithm is used to search all the regions (connected components) in the binary image generated by the thresholding step, assigning a unique label to each identified region. A connected component in a binary image is an area of pixels with value 1, which are connected per a pre-defined connectivity type. In astronomical source finding, the 8-connectivity type is used, i.e. it is assumed that two samples are connected if they are adjacent in the horizontal, vertical or diagonal directions. There are several CCL algorithms proposed in literature, some of which exploit parallel processing [22, 23]. The existing algorithms can be classified into three major groups:

Multi-pass

algorithms usually require several passes over the image before reaching the final labels. The number of passes depends on the content. These algorithms may use specific techniques to reduce the number of passes, such as a label connection table.

Two-pass

algorithms require only two passes over the image. These algorithms usually operate in three distinct phases:

• •

Scanning phase: the image is scanned in the first pass to assign provisional labels to object pixels and to record the equivalence information among provisional labels.

• •

Analysis phase: this phase analyses the label equivalence information to determine the final labels.

• •

Labelling phase: this phase assigns final labels to object pixels using a second pass through the image.

One-pass

algorithms scan the image to find an unlabelled object pixel and then assign the same label to all connected object pixels. As they perform a single pass over the image, these methods require irregular access to pixels, which usually leads to poor performance.

The computational performance of CCL algorithms depends not only on their approach and the number of passes over the image, but also on the input data characteristics (e.g. the number of connected components present in the image). The YACCLAB (Yet Another Connected Components Labeling Benchmark; [24]) is one example of a C++ open source framework which evaluates various CCL algorithms for different data types.

3.3.1 Efficient Component Labeling

Here we implement a multithreaded version of the 8-connectivity CCL algorithm provided in OpenCV [25], a reference open source library for computer vision. This algorithm is based on the approach of [26], which is the two-pass CCL algorithm briefly described in § 3.3. Although this solution requires two passes over the image, it tends to be significantly faster than one-pass algorithms that use irregular and slow memory access patterns. In addition to the use of multithreaded processing for CCL, the overall source find algorithm performance was improved by merging some related source finding steps into the processing loop of the two-pass CCL function. These optimisations are summarised here and described in more detail below:

Label non-shifted FFT output image as if its quadrants were properly shifted; 2. 2.

Use multi-threaded implementation for faster performance; 3. 3.

Merge the thresholding step into the labelling function to avoid allocating memory for the binary image; 4. 4.

Merge the valid source detection step with the last stage of the labelling function to save one full pass over the label map; 5. 5.

Merge labelling of positive and negative sources in a single CCL procedure using only two image passes;

Optimisation 1

Interferometric images are made by Fourier transforming the native measurement data from radio interferometers. Typically this is done using the Fast Fourier Transform (FFT), implementations of which require a quadrant shift to be made on the output data grid in order to correctly align data values. Since this shift requires copying the image data with non-sequential accesses, it tends to be computationally expensive. For example, in the STP, it was observed that FFT quadrant shift may increase the imager running time up to 20%. In this context, the first optimisation implemented here intends to avoid the need to perform quadrant shifting prior to the source finding process.

An example of a binary image that results from thresholding a non-shifted image is shown in the left image of Figure 2. The corresponding shifted version is presented in the right image of Figure 2. If the non-shifted image was used as input to the original CCL function, without the proposed changes, the labelling result would be as illustrated in Figure 3. Six components would be detected in this example (see the left image of Figure 3). When the labelled image is shifted the results are not as expected, because multiple labels appear as being assigned to the same object (see right image of Figure 3).

To fix this issue, the CCL algorithm has been modified so that it produces a label map as shown in Figure 4. As can be observed, the shifted version of the label map (right image of Figure 4) presents only 3 labels as expected. To explain the proposed change, we first review the phases of the CCL algorithm: the scanning, analysis and labelling. In the scanning phase, provisional labels are assigned to object pixels. In this phase, multiple labels can be assigned to one connected component, resulting in a provisional label map as illustrated in the left side of Figure 5, that contains 4 provisional labels. After determining the final labels in the analysis phase and performing the final labelling phase, the resulting label map is the expected one, with 3 labels, as illustrated in the right side of Figure 5.

The proposed solution modifies the scanning phase of CCL to assume discontinuities at the horizontal and vertical dashed lines in the middle of the image (see the left image of Figure 6). Thus, when assigning the provisional labels, different numbers are assigned to samples separated by these lines. An example illustrating the result of the modified scanning phase for the non-shifted image is given in the left side of Figure 6. In this example, labels 4 and 6 are created due to the horizontal dashed line in the image centre representing a discontinuity.

A second change is related to the label equivalence information created during the scanning phase. An equivalence is defined between provisional labels that are assigned to the same connected component (or object). Since the objects in the non-shifted image can be partitioned at the image margins (see for instance, the object with label 2 in the right image of Figure 6), an equivalence between provisional labels assigned to objects at the image margins shall be created. In practice, it is added an new step that creates equivalences between provisional labels touching opposite image margins. This process is denominated as border merging and it is also used on parallel implementations of CCL. After border merging, the analysis and labelling phases are applied, producing the expected result as shown in the right image of Figure 6.

It is important to refer that the changes introduced by this optimisation may slightly increase the computational complexity of the labelling stage, typically up to a maximum of 3%. This increase is mostly due to the additional conditions required in the scanning phase. The image data continues being scanned sequentially in the memory thus avoiding a high performance penalty.

Optimisation 2

The second proposed improvement is the use of parallel processing to reduce the running time of CCL. The proposed algorithm uses a multi-threaded implementation of CCL, that parallelises the scanning and labeling phases. Parallelisation of the analysis phase is not possible due to the way the equivalence information is processed to determine the final labels. This is not a significant issue as the analysis phase presents a low complexity, especially when the number of provisional labels is not significant.

For an efficient parallel processing, the image should be partitioned in horizontal or vertical strips, depending whether the image is stored in row-major or column-major order, and each thread assigned to a distinct strip. Other partitioning approaches, like the use of quadrants are not recommended and would create a high performance penalty, since the threads would be assigned to non-contiguous blocks of data in the system’s memory.

While the labeling phase can be easily implemented in parallel by assigning different strips of the binary image to each thread, the scanning phase requires some extra work. For the scanning phase, the image also can be divided into strips that are assigned to different threads. However, this procedure changes the labeling results of the algorithm, specifically when the one connected component is defined across two strips. In that case, different provisional labels are assigned by each thread to each partition of the object. An example is illustrated in the right side of Figure 7 for a image represented in row-major order. This issue can be resolved by performing border merging over the horizontal strip borders, with the exception of the dashed lines in the middle of the image. Figure 7 shows the border merging lines in the single-threaded and multi-threaded implementations represented by the thick solid lines.

Optimisation 3

The third proposed optimisation is to perform the image thresholding during the scanning phase of CCL. Instead of generating a new array to store the binary image, the CCL function receives the original dirty image as input. When it scans the image, each sample is compared with the analysis threshold, the binary value is determined and the provisional label is generated.

Optimisation 4

The fourth proposed optimisation merges the last step of valid label detection with the labeling phase of CCL. For detecting valid labels, the maximum sample value of each object or source is compared with the detection threshold. A positive source that contains at least one sample above the positive detection threshold is considered valid. A negative source is valid only if it contains at least one sample below the negative detection threshold. Merging this procedure with the last phase of the CCL allows one to avoid an extra pass over the label map and thus to reduce the number of memory accesses. Its implementation uses parallel processing, performed in the same loop as the labeling phase of CCL.

Optimisation 5

As a final optimisation, the proposed CCL function was designed to perform the search for positive and negative sources at the same time. This optimisation will avoid calling the CCL function two times for each of both positive and negative sources, which would result in four passes over image data. The proposed solution only requires two passes to label both positive and negative sources. To make it possible, we use two structures to independently store the equivalence information for positive and negative sources. An advantage of this solution is that the same label map can be used to represent the positive and negative sources, as they do not overlap. This allows one to significantly reduce the amount of memory usage as well as the number of memory reads and writes. We note that the alternative is to perform two runs of CCL to detect the positive and negative sources and then combine the two resulting label maps.

3.3.2 Data Passes

The optimisations described above result in a significant improvement of the performance of the detection algorithm for positive and negative sources, reducing the number of image passes to two, which is the minimum required for the CCL algorithm. After detection of the valid sources, the algorithm computes the first and second statistical moments of each valid source. To compute these moments, a third pass over the label map is required. In this pass, some auxiliary data required for the Gaussian fitting step is also computed, namely the bounding box and the number of samples for each region. Optionally, the third pass may be also used to remove invalid sources (those below the detection threshold) from the label map. The final label map may be desired for instance for debugging purposes, e.g. to be observed and analysed by the user.

Table 1 summarises the major steps of the labeling algorithm and the number of required passes over the image (which may refer to the dirty image or label map), specifically for two implementation approaches: the straightforward and the optimised. An image pass implies write operations on the image unless otherwise stated, e.g “only read”. The straightforward approach refers to a trivial implementation of the labeling function, while the optimised approach corresponds to the proposed solution that minimises the number of data writings and readings.

It can be observed that the proposed optimised approach requires only 3 image passes while the straightforward implementation uses 12 image passes, although not all of them involve write operations. Such a large reduction in the number of passes, relative to the straightforward implementation, results in a significant performance improvement for the source finding. We do not believe that it is possible to reduce the number of passes further due to the data dependencies in the processing.

In the first image pass of the optimised approach the thresholding step and the scanning phase of CCL are performed for both positive and negative sources. The second image pass of the optimised solution is used for the final labeling phase of CCL and the detection of valid sources. Finally, the third pass is used to compute the statistical moments of each source, compute auxiliary data about each source and optionally to update the label map, removing any invalid sources. The optimised solution does not need to combine the label maps of negative and positive sources, since the labeling function directly generates a combined source map.

3.4 Gaussian fitting

As a final step in the source finding algorithm, Gaussian functions are fitted to each detected source or island. This procedure consists of solving non-linear least squares problems to estimate the parameters of the Gaussian functions that fit each island. Gaussian fitting has been implemented using the Ceres Solver library [27].

The Ceres Solver first models the problem and then solves it. For modeling, the Ceres Solver provides three main methods to compute the derivatives:

• •

Analytic: The derivatives shall be manually implemented by the user.

• •

Numeric: Ceres Solver numerically computes the derivatives using finite differences.

• •

Automatic: Ceres Solver automatically computes the analytic derivatives using C++ templates and operator overloading.

In the algorithm developed, analytic or automatic derivatives can be used, depending on the chosen option in the input configuration file to the program. We did not implement numeric derivatives because they present a larger computational complexity. In this work it was observed that the analytic derivatives tended to perform more efficiently than automatic derivatives.

For the solver step, the Ceres Solver provides several methods. Through the configuration file, it is possible to choose between the trust region and line search minimiser methods to solve the non-linear least squares problem. The trust region method approximates the objective function using a model function (often a quadratic) over a subset of the search space known as the trust region. By default, we use the Levenberg-Marquardt algorithm and the dense QR factorization routines of the Eigen library to solve the trust region problem. For the line search approach, the BFGS or LBFGS method may be used to choose a search direction. These are a generalisation of the conjugate gradient method for non-linear functions.

The Gaussian fitting step computes the following list of parameters for each detected island:

amplitude

(double) the amplitude parameter of the Gaussian function;

x_center, y_center

(double) the x, y coordinates of the Gaussian function centre in the image;

semimajor

(double) the length of the major Gaussian axis in pixels;

semiminor

(double) the length of the minor Gaussian axis in pixels;

theta

(double) the position of the semimajor axis measured counter-clockwise from the x axis (in radians).

In addition to the Gaussian parameters derived by least-squares optimisation, an initial Gaussian estimation based on the first and second statistical moments is outputted. Thus, each detected island is accompanied by two sets of Gaussian parameters. Jointly, the following data about each detected island is provided:

sign

(integer) represents whether the transient source is positive or negative. Possible values are $+1$ and $-1$;

val

(double) represents the extremum (maximum or minimum) pixel value for the (positive or negative) transient source;

x_idx, y_idx

(integer) represent the x, y index values (i.e. pixel coordinate in the image) of the extremum value;

num_samples

(integer) indicates the number of pixels in the island;

In terms of computational complexity, the Gaussian fitting step usually does not present a significant impact on the overall running time of the pipeline, especially when the number of detected islands is not significant. However, when the number of detected islands increases the computational cost of this step may be noticeable.

4 Tests and benchmarks

Here we present and analyse the benchmark results of the main functions that constitute the source finding algorithm developed. Benchmarks were performed on an High Performance Computer (HPC) with two processor sockets, specifically the Intel Xeon CPU E5-2650 v4 working at 2.20GHz and using 504 GB of RAM. Each CPU has 12 physical cores that provide 24 threads due to Hyper-Threading (HT) technology, summing a total of 48 threads.

The HT technology doubles the number of available threads in a multi-core CPU by duplicating certain parts of the processor which allows each processor core to appear as the usual physical core and an extra logical core to the operating system. While hyper-threading tends to provide a superior computing performance for certain cases, the parallel efficiency of each thread tends to be inferior due to hardware resources of one physical core being shared between two threads. For this reason, this technology has been disabled for the performed experiments. The amount of threads used in the HPC was thus 24.

The benchmarking tests discussed here are based on the measured running times, which provide a rough idea of the relative computational complexity of the tested functions. It is important to note that these timings may vary between distinct benchmarking instances, even on the same hardware, as they may be affected by the operating system tasks running concurrently and as well as its status.

The following experiments are presented and discussed in this section:

• •

Median - evaluate the running time of the proposed median functions discussed in § 3.1 for different sized matrices of random data, comparing with other median implementations, as well as the single- and multi-threaded execution.

• •

Standard deviation - evaluate the running time of the standard deviation function proposed in § 3.2 for different sized matrices of random data, comparing with other standard deviation implementations, as well as the single- and multi-threaded execution.

• •

Steps of source find - evaluate the running times of the main steps of the source find algorithm proposed here using simulated astronomical images.

• •

Source find - evaluate the total running time of the source find algorithm proposed here using simulated astronomical images and comparing single- and multi-threaded executions.

The simulated astronomical data used for source find testing contains 3 gaussian-shaped sources. Such reduced number of sources was used because we are using residual images generated in the context of the STP development. There is not any limitation in the algorithm relative to the number of sources present in the image. In fact, we used a larger number of sources in the experiments of § 5 where the algorithm is compared with other source finders.

Figure 8 presents the running time of different implementations of the median function used for background level estimation, specifically:

• •

Nth_element method - provided in the C++ standard library for $k^{th}$ smallest element selection based on quickselect algorithm;

• •

arma::median method of the Armadillo library [28] - a high quality linear algebra library for the C++ language used in this work for matrix representation. It aims towards a good balance between speed and ease of use;

• •

binmedian method - proposed in § 3.1;

• •

binapprox method - proposed in § 3.1;

The respective speed up gains, using the Nth_element function as reference method are given in Figure 9. From these results, it may be concluded that the binapprox method provides the best performance, presenting speed up gains between 12 and 18 relative to the reference Nth_element function. The exceptional running performance for the $2^{11}\times 2^{11}$ image size can be related with the cache memory optimisation. As such small image sizes mostly fit into the cache memory, the associated execution performance is easily affected by the cache optimisation, which depends on the hierarchical cache sizes, number of threads and allocated buffer sizes.

Despite its superior performance, it is important to note that the binapprox method does not compute the exact median value. If exact median is required, the binmedian method provides the best results, presenting a speed up gain between 6 and 8 relative to the Nth_element function. Regarding the method provided by Armadillo, its performance is close to the one that uses the Nth_element function. This is not a surprising result since Armadillo’s implementation of median is based also on the Nth_element function.

It is important to remember that although all the tested algorithms shown in Figures 8 and 9 use multi-core processing, they do not present the same parallel efficiency. To estimate the speedup gains provided by parallel implementation using TBB, the previously presented median functions were run in single-threaded mode and the speedup gains were computed. The achieved multi-threading speedup gains are shown in Figure 10. It can be observed that the methods based on the Nth_element function (which includes Armadillo) present a multi-threading speedup between 2 and 3, corresponding to a small parallel efficiency of between 8% and 12.5% in the 24-core HPC. These methods use a parallel version of the Nth_element function provided by the Standard C++ Library, which does not scale efficiently due to the characteristics of its algorithm based on quickselect. Differently, the binmedian and binapprox methods benefit significantly from parallel processing, presenting speedup gains between 12 and 17, which correspond to a parallel efficiency between 50% and 70% in the used computer.

In fact, the good performance of the binmedian and binapprox methods is possible due to the high parallel efficiency of these functions. When using a single-core system, the advantage of the binmedian method over the Armadillo median function is irrelevant. In some cases, the single-threaded version of binmedian can be even worse than the median methods based on the Nth_element function. The use of multi-core processing is thus the key to take advantage of the binmedian and binapprox algorithms.

Another important function of STP, that is used also as a step of the binmedian and binapprox methods is the standard deviation function. Figure 11 presents the running times of two implementations of standard deviation function, namely the one provided by Armadillo library and the proposed solution based on TBB (see § 3.2) executed in both single-threaded and multi-threaded modes. The speed up results of the proposed function relative to the Armadillo method are presented in Figure 12. It can be observed that the proposed TBB-based solution presents significant speed up gains, being approximately 4 and 12 times faster than Armadillo for the single- and multi-threaded modes, respectively. For the smaller image sizes, the speed up results present a larger variation, due to the fact that these sizes mostly fit into the cache memory, being easily affected by the cache optimisation techniques.

Usually the mean quantity is computed as the first step of the standard deviation function. However, the proposed function uses an alternative approach that avoids the mean computation. The developed function computes two accumulation values as explained in § 3.2 using a single pass over the array. These accumulation values are then used to derive both the mean and standard deviation of the data. The fact that the developed function provides both these quantities at once is another important advantage, since it avoids the necessity to call the mean function when both the mean and standard deviation quantities are required.

Figure 13 presents the benchmarking results of the full source finding algorithm, performing the background level estimation (based on binapprox function), RMS estimation (with sigma clipping using 5 iterations), source labelling and gaussian fitting. We used simulated images with sizes of $2^{12}\times 2^{12}$, $2^{14}\times 2^{14}$ and $2^{16}\times 2^{16}$. As can be observed in the referred figures, the gaussian fitting step tends to use an insignificant running time. This is mostly because the test-data has few sources. Its timing shall increase in the presence of a larger number of sources.

The most complex step of the source finder is the RMS estimation that needs to perform 5 passes over the image to perform sigma clipping. Due to the optimisations described in this paper, the labelling step performs quite efficiently, presenting execution times inferior to both the background level and RMS estimation functions, for most image sizes. When comparing the results for different image sizes, one observes that the relative computational cost of each step does not vary significantly, although the respective absolute timings increase for larger images sizes as expected.

To analyse the parallel performance of the source finding algorithm, Figure 14 presents the running times of both the single-threaded and multi-threaded execution of the algorithm. The running times of the multi-threaded algorithm vary between 4 milliseconds and 4 seconds, while the single-threaded algorithm varies between 10 milliseconds and 44 seconds, for the image sizes $2^{10}\times 2^{10}$ and $2^{16}\times 2^{16}$, respectively. As can be observed in Figure 14, these running times tend to increase exponentially with the image width, increasing linearly with the number of image samples.

The performance gains of the multi-threaded implementation relative to the single-threaded one vary between 3 and 11, as shown in Figure 15, corresponding to a parallel efficiency of 12.5% and 46%, respectively. The multi-threading gains are mostly achieved by the parallel implementation of the median (binapprox method), RMS estimation and labelling functions, being more noticeable for larger image sizes. The fact that the smaller images better fit into CPU cache memory may reduce the parallel efficiency of the algorithm, for instance, due to the false sharing issue caused when multiple threads write data into the same cache line.

5 Comparative benchmarks

In addition to the detailed benchmarking of our C++ implementation described above, we also performed some comparative benchmarks with the other widely-used source-finding packages, some of which were outlined in § 2.2, namely:

• •

PySE - LOFAR-TKP Python Source-Extractor [10, 11];

• •

SEP - fork of the S-Extractor package [9];

• •

Python ref. code - reference Python code we implemented as a guide for the C++ implementation;

• •

Ceres - the C++ source find algorithm proposed in this paper;

Test-data consisted of randomly generated images of different sizes containing 64 gaussian-shaped sources. The generated test images were correctly shifted, so the labelling algorithm was compiled considering this fact, i.e. these tests are not accounting for the additional gains that can be achieved by avoiding the quadrant shift operation after FFT usually performed in the radio imaging pipeline.

For these benchmarks we focused on two key stages, the background and noise estimation, and source-extraction including gaussian fitting. We configured the packages to perform these tasks on in-memory arrays of image-data — this reflects the likely ‘real-world’ usage of a high-throughput source-extraction pipeline, and removes the contaminating (and potentially benchmark-dominating) effects of different input-output data-access strategies. However, the original S-Extractor package [9] provides only a command-line interface, with no library API for direct calls, and processes images in a piecewise manner to reduce memory overheads. We instead benchmark the SEP library [29], a recent fork of S-Extractor which modifies the original S-Extractor codebase to run directly on in-memory arrays and provides Python bindings. We do not provide benchmarks for the Aegean package as the interface provides no easy method for supplying image-arrays directly. In the interests of simplicity the benchmarks are all run via Python scripts, with C/C++ routines accessed via Python bindings. Consequently, the performance of these packages will be slightly reduced compared to an end-to-end C++ pipeline due to the overhead involved in data-passing. 111It is worth noting that the array-ordering can have a significant influence on performance when using Python bindings — SEP will only accept arrays in the C-style ordering (the numpy default), but the Python bindings developed for the C++ library described herein are significantly more performant (performance improved by around 1.5–2x) when passed numpy arrays in Fortran ordering, as used internally by the Armadillo routines.

Figures 16 and 17 depict the results of timing benchmarks for the source-extraction and fitting process, as performed on images of varying size. The experiments not applying source fitting are labelled with the abbreviation “ex.”, while the full process of source-extraction and fitting is labelled with “ex.+fit”. Figure 16 shows the results for all packages tested. The first striking result is how much slower the Python implementations are compared to the C/C++ packages — this is unexpected given that the Python implementations make heavy use of the Numpy / Scipy libraries which in turn make heavy use of C-extensions. We can observe that the improvement from a dedicated C++ routine is fairly dramatic, with a speedup factor of several hundreds. We also note that performance of our reference Python implementation is significantly slower than PySE — this seems to be due to performing additional array copies and highlights the challenges of working efficiently with large arrays in Python.

Figure 17 depicts the results for the C/C++ packages only, with tests performed at varying image sizes of widths from $2^{10}$ – $2^{14}$ (1024 – 16384) pixels. The novel C++ implementation is markedly faster than SEP for large image sizes, completing the source-extraction approximately 13 times faster than SEP for an image width of $2^{14}$ pixels. For the proposed Ceres package we plot two lines, one for extraction and moments-estimation alone — which is directly comparable to the results given by SEP — and a second line which also includes the time required for Gaussian fitting of each source-island. Since the number of sources is constant, we expect an efficient implementation of the island-fitting to perform in constant time for a constant number of sources regardless of the image size, and this is true to a good approximation for the Ceres implementation, with the fitting procedure adding a near-constant and quite small additional execution time across image-sizes (some milliseconds). We note that we were unable to perform benchmarks of SEP at image-widths above $2^{14}$ as attempting to do so resulted in a segmentation fault.

It is important to note that the timing results of the proposed Ceres package in these experiments include the parallel processing techniques. However, further experiments using a single thread revealed that the developed package continues presenting a superior performance than SEP.

6 Conclusions

Upcoming observatories, such as the Square Kilometre Array, demand high efficient algorithms for object detection in astronomical images. In this paper, we discuss a new high performant solution for source finding, which incorporates a number of advanced techniques to accelerate its algorithm and adapt it for multi-core processing systems.

Main contributions include: fast median computation for background noise estimation based on data binning which presents better parallel scalability than traditional quick-select approach; fast sigma clipping algorithm for RMS estimation which minimises the amount of computations and memory accesses; and an efficient and highly parallelisable CCL algorithm that supports detection of positive and negative sources with reduced computational complexity and processes non-shifted images as outputted by FFT from the previous imaging pipeline stage. Timing benchmarks demonstrate the superior computational performance of the techniques proposed for source finding, in particular when using multi-core processing. Comparative benchmarks against existing source finding algorithms show that the developed C++ solution presents a state-of-the-art performance when processing in-memory arrays of image-data.

Acknowledgements

The authors gratefully acknowledge support from the UK Science and Technology Facilities Council.

The authors gratefully acknowledge the support of the ENGAGE SKA HPC cluster at IT-Aveiro. ENGAGE SKA (POCI-01-0145-FEDER-022217) is funded by COMPETE 2020 and FCT, Portugal.

The authors gratefully acknowledge the developers of the following packages that were used in this work: Armadillo222http://arma.sourceforge.net [28], Google Test333https://github.com/google/googletest, Google Benchmark444https://github.com/google/benchmark, cnpy555https://github.com/rogersce/cnpy, FFTW3666http://www.fftw.org [30], pybind11777https://github.com/pybind/pybind11, TBB888https://www.threadingbuildingblocks.org [20], RapidJSON999https://github.com/miloyip/rapidjson, TCLAP101010http://tclap.sourceforge.net, spdlog111111https://github.com/gabime/spdlog, Eigen121212http://eigen.tuxfamily.org [31], Ceres Solver131313http://ceres-solver.org, Astropy141414http://www.astropy.org [32, 33], Numpy151515http://www.numpy.org [34], Scipy161616http://www.scipy.org [16], Matplotlib171717https://matplotlib.org/ [35], Jupyter Notebook181818https://jupyter.org/, Spack191919https://spack.io [36], GNU/Linux operating system and tools202020https://www.gnu.org.