© A. Men’shchikov 2021

1. Introduction

Multiwavelength far-infrared and submillimeter dust continuum observations with the large space telescopes Spitzer, Herschel, and Planck in the past decades greatly increased the amount and improved the quality of the available data in various areas of astrophysical research. Observed images with diffraction-limited angular resolutions and high sensitivity reveal an impressive diversity of the enormously complex structures, covering orders of magnitude in intensities and spatial scales. The images feature foremost the bright fluctuating backgrounds, omnipresent filaments, and huge numbers of sources of different physical nature, all blended with each other, whose appearance and resolution are often markedly different at short and long wavelengths. The massive amount of information that is coded in the fine structure of the observed images must contain clues to the complex physical processes taking place in space, but these clues are extremely difficult to decipher. It is quite clear that the era of these high-quality data from space telescopes and large ground-based interferometers, such as the Atacama Large Millimeter/submillimeter Array (ALMA), requires more sophisticated tools for their accurate analysis and correct interpretation than those developed for the lower-quality images of the past. Adequate extraction methods must be explicitly designed for the multiwavelength imaging observations with highly dissimilar angular resolutions across wavebands. They must also be able to handle the bright filamentary backgrounds that vary on all spatial scales, whose fluctuation levels differ by several orders of magnitude across the observed images.

The source- and filament-extraction methods are growing in numbers. In the area of star formation, a new method was published every year or two within the seven-year period after the launch of Herschel. Rosolowsky et al. (2008) devised dendrograms to describe the hierarchical structure of clumps observed in the data cubes from molecular line observations. The method carries out topological analysis of image structures by isophotal contours at varying intensity levels and represents them graphically as a tree. Molinari et al. (2011) created cutex to extract sources in star-forming regions observed with Herschel. The method analyzes multidirectional second derivatives of the observed image to detect sources, and it measures them by fitting elliptical Gaussians on a planar background to their peaks. Men’shchikov et al. (2012) developed getsources, the multiwavelength source extraction method for the Herschel observations of star-forming regions. The method spatially decomposes images, combines them into wavelength-independent detection images, subtracts the backgrounds of detected sources, and measures the sources, deblending them when they overlap. Kirk et al. (2013) presented csar for the Herschel images. The method analyzes areas of connected pixels that are bound by closed isophotal contours, descending to a predefined background level and partitioning peaks at their lowest isolated contours into sources. Berry (2015) created fellwalker to identify clumps in submillimeter data cubes. The method finds image peaks by tracing the line of the steepest ascent and identifies sources as the hill with the highest value found for all pixels in its neighborhood. Sousbie (2011) produced disperse to identify structures in the large-scale distribution of matter in the Universe. The method applies the computational topology to trace filaments and other structures. Men’shchikov (2013) developed getfilaments to improve the source extraction with getsources on the filamentary backgrounds observed with Herschel. The method separates filaments from sources in spatially decomposed images and subtracts them from the detection images, thereby reducing the rate of spurious sources. Schisano et al. (2014) devised a Hessian matrix-based approach to extract filaments in Herschel observations of the Galactic plane. The method analyzes multidimensional second derivatives to identify filaments and determine their properties. Clark et al. (2014) presented rht to characterize fibers in the interstellar H I medium. The method has been applied to various observations of diffuse H I, revealing alignment of the fibers along magnetic fields. Koch & Rosolowsky (2015) published filfinder to identify filaments in the Herschel images of star-forming regions. the method applies a mathematical morphology approach to isolating filaments in observed images. Juvela (2016) presented tm to trace filaments in observed images. The method matches a predefined template (stencil) of an elongated structure at each pixel of an image by shifting and rotating the template and analyzing the parameters of the matches.

These extraction methods all have several important issues. Sources and filaments are handled completely independently by these methods, although numerous Herschel observations have demonstrated that there is a tight physical relation between them. Most sources are found in filamentary structures, and the corresponding starless, prestellar, and protostellar cores are thought to form inside the structures that are created by dynamical processes, magnetic fields, and gravity within a molecular cloud. All major structural components of the observed images, that is, the background cloud, filamentary structures, and sources, are heavily blended with each other; curved filaments are even blended with themselves. The degree of their blending increases at longer wavelengths with lower angular resolutions, which increases the inaccuracies in their detections, measurements, and interpretations.

Most of the extraction methods focus on detecting structures, whereas the most important and difficult problem is measuring them accurately. Numerous algorithms only partition the image between sources and do not allow them to overlap, although deblending of the mixed emission of the structural components is an indispensable property of an accurate extraction method. For best detection and measurement results, source-extraction methods must be able to separate underlying filamentary structures and filament-extraction methods must be able to separate sources. The existing source- and filament-extraction methods use completely different approaches, and the quality of their results is expected to be very dissimilar.

It seems unlikely that methods that are based on very different approaches would give consistent results in terms of detection completeness, number of false-positive detections, and measurement accuracy. In practice, various methods do perform very differently, as can be shown quantitatively on simulated benchmarks for which the properties of all components are known. This highlights the need of systematic comparisons of different methods in order to understand their qualities, inaccuracies, and biases. The danger is real that numerous uncalibrated methods are applied for the same type of star formation studies, which would give inconsistent, contradictory results and incorrect conclusions. This would create serious, lasting problems for the science.

Source- and filament-extraction methods are the critically important tools that must be calibrated and validated using a standard set of benchmark images with fully known properties of all components before they are applied in astrophysical contexts. It would be desirable to use the same extraction tool to exclude any biases or dissimilarities that are caused by different methods. If a new extraction method is to be used, it must be tested on standard benchmarks to ensure that its detection and measurement qualities are consistent or better. This approach is usually practiced within research consortia, but this does not solve the global problem that the results obtained from the same data by different consortia or research groups using different tools may still be affected by the uncalibrated (or suboptimal) tools that were used.

This paper presents getsf, a new multiwavelength method for extracting sources and filaments. It also describes a realistic simulated benchmark, resembling the Herschel images of star-forming regions, which is used below to illustrate the method and in a separate paper (Men’shchikov 2021) to quantitatively evaluate its performance. The multiwavelength benchmark simulates the images of a dense cloud with strong nonuniform fluctuations, a wide dense filament with a power-law intensity profile, and hundreds of radiative transfer models of starless and protostellar cores with wide ranges of sizes, masses, and profiles. The simulated benchmark with fully known parameters allows quantitative analyses of extraction results and conclusive comparisons of different methods by evaluating their extraction completeness, reliability, and goodness, along with the detection and measurement accuracies. The multiwavelength images can serve as the standard benchmark problem for other source- and filament-extraction methods, allowing researchers to perform their own tests and choose the most reliable and accurate extraction method for their studies. Instead of publishing benchmarking results for some of the existing methods, it seems a better idea to provide researchers with the benchmark¹ and a quality evaluation system (Men’shchikov 2021) to enable comparisons of the methods of their choice. In practice, this approach of having own experience is much more convincing and it allows a consistent evaluation of newly developed methods.

The new source and filament extraction method getsf represents a major improvement over the previous algorithms getsources, getfilaments, and getimages (Men’shchikov et al. 2012; Men’shchikov 2013, 2017, hereafter referred to as Papers I, II, and III); throughout this paper, the three predecessors are collectively referred to as getold. The new method (Fig. 1) consistently handles two types of structures, sources and filaments, that are important for studies of star formation, separating the structural components from each other and from their backgrounds. All major processing steps of getsf employ spatial decomposition of images into a number of finely spaced single-scale images to better isolate the contributions of structures with various widths. The method produces accurately flattened detection images with uniform levels of the residual background and noise fluctuations. To detect sources and filaments, getsf combines independent information contained in the multiwaveband single-scale images of the structural components, preserving the higher angular resolutions. Then getsf measures and catalogs the detected sources and filaments in their background-subtracted images. The fully automated method needs only one user-defined parameter, the maximum size of the structures of interest to extract, constrained by users from the input images on the basis of their research interests.

This work follows Papers I–III in advocating a clear distinction between the words source and object, unlike many publications in which it is implicitly assumed that the two are completely equivalent. “Source” is used in the context of the source extractions and statistical analysis of their results, and “object” is only used in the context of the physical interpretation of the extracted sources. In this paper, the sources are defined as the emission peaks (mostly unresolved) that are significantly stronger than the local surrounding fluctuations, indicating the presence of the physical objects in space that produced the observed emission. The implicit assumption that an unresolved far-infrared source on a complex fluctuating background contains emission of just one single object is invalid in general. Too often, an emission peak is actually a blend of many components, produced by different physical entities. This is illustrated by the recent images of the massive star-forming cloud W43-MM1 (Motte et al. 2018), obtained with the ALMA interferometer. This object appears as a single source in the Herschel images, even with the 5.6 and 11.3″ resolutions at 70 and 160 μm. However, the ALMA image (Sect. 4.8) displays a rich cluster of much smaller sources that are unresolved or just slightly resolved even at the 0.44″ resolution.

Section 2 describes the new multiwavelength benchmark for source- and filament-extraction methods, resembling the Herschel observations of star-forming regions. Section 3 presents getsf, the new source- and filament-extraction method, employing separation of the structural components. Section 4 illustrates the performance of getsf on a large variety of images that were obtained with different telescopes in a wide spectral range, from X-rays to millimeter wavelengths. Section 5 describes all strengths and limitations of getsf. Section 6 presents a summary of this work. Appendix A discusses inaccuracies of the surface densities and temperatures, derived by spectral fitting of the images. Appendix B describes the single-scale spatial decomposition that is used by getsf in its processing steps. Appendix C gives details on the software.

In this paper, images are represented by capital calligraphic characters (e.g., 𝒜,  ℬ,  𝒞) and software names and numerical methods are typeset slanted (e.g., getsf) to distinguish them from other emphasized words. The curly brackets are used to collectively refer to either of the characters, separated by vertical lines. For example, {a|b} refers to a or b and {A|B}_{a|b}c expands to A_{a|b}c or B_{a|b}c, as well as to A_ac, A_bc, B_ac, or B_bc.

2. Benchmark for source and filament extractions

Realistic multiwavelength, multicomponent images of a simulated star-forming region were computed to present getsf in this paper and to compare its performance with the previous benchmark that was used in Papers I and III. The benchmark images were created for all Herschel wavebands (at λ of 70, 100, 160, 250, 350, and 500 μm). They consist of independent structural components: a background cloud ℬ_λ, a long filament ℱ_λ, round sources 𝒮_λ, and small-scale instrumental noise 𝒩_λ:

$$\begin{array}{c}\hfill {\mathcal{H}}_{\lambda }={\mathcal{B}}_{\lambda }+{\mathcal{F}}_{\lambda }+{\mathcal{S}}_{\lambda }+{\mathcal{N}}_{\lambda }={\mathcal{C}}_{\lambda }+{\mathcal{S}}_{\lambda }+{\mathcal{N}}_{\lambda },\end{array}$$(1)

where 𝒞_λ = ℬ_λ + ℱ_λ is the emission intensity of the filamentary background. All simulated images were computed on a 2″ pixel grid with 2690 × 2690 pixels, covering 1.5° ×1.5° or 3.7 pc at a distance D = 140 pc of the nearest star-forming regions (e.g., those in Taurus or Ophiuchus).

2.1. Simulated filamentary background

An image of the background surface density was computed from a purely synthetic scale-free background 𝒟_A (cf. Paper I), with N_H2 ∼ 2.7 × 10²⁰ to 5 × 10²² cm⁻² that had uniform fluctuations across the entire image. To simulate complex astrophysical backgrounds with strongly nonuniform fluctuations (e.g., Könyves et al. 2015), 𝒟_A was multiplied by a circular shape 𝒫 with a radial profile defined by Eq. (2) below (with Θ = 1500″ and ζ = 2), normalized to unity and centered on the image; finally, a constant value of 1.5 × 10²¹ cm⁻² was added to increase the minimum value. The surface densities of the resulting background cloud image 𝒟_B (Fig. 2) are 1.5 × 10²¹ to 4.8 × 10²² cm⁻² and the fluctuations differ by approximately two orders of magnitude. The total mass of the cloud is M_B = 1.78 × 10³ M_⊙.

To simulate filamentary backgrounds, a long spiral filament was added to the background cloud 𝒟_B. The spiral shape was chosen so that the filament occupied various areas of the cloud with very different surface densities and to cause the filament to be blended (with itself) to some extent. The spiral filament image 𝒟_F has a crest value of N₀ = 10²³ cm⁻², a full width at half-maximum (FWHM) W = 0.1 pc, and a radial profile similar to those observed with Herschel in star-forming regions (e.g., Arzoumanian et al. 2011, 2019),

$$\begin{array}{c}\hfill N_{{\mathrm{H}}_2}(\theta )=N_0{(1+(2^{1/\zeta }-1){(\theta /\mathrm{\Theta })}^2)}^{-\zeta },\end{array}$$(2)

where θ is the angular distance, Θ is the structure half-width at half-maximum, and ζ is a power-law exponent. With Θ = 75″ (or 0.05 pc at D = 140 pc) and ζ = 1.5, this Moffat (Plummer) function approximates a Gaussian of 0.1 pc (FWHM) in its core and it transforms into a power-law profile N_H2(θ) ∝ θ⁻³ for θ ≫ Θ. The filament mass M_F = 3.04 × 10³ M_⊙ and length L_F = 10.5 pc correspond to the linear density Λ_F = 290 M_⊙ pc⁻¹. The resulting surface densities 𝒟_C = 𝒟_B + 𝒟_F of the filamentary cloud are in the range of 1.7 × 10²¹ to 1.4 × 10²³ cm⁻² (Fig. 2), and its total mass is M_C = 4.82 × 10³ M_⊙.

To approximate the nonuniform line-of-sight dust temperatures of the star-forming clouds observed with Herschel (e.g., Men’shchikov et al. 2010; Arzoumanian et al. 2019), an image of average line-of-sight temperatures was improvised as

$$\begin{array}{c}\hfill {\mathcal{T}}_{\mathrm{C}}=200{({10}^{-20}{\mathcal{D}}_{\mathrm{C}}+20)}^{-1}+15\mathrm{K}.\end{array}$$(3)

The pixel values of the resulting temperature image 𝒯_C range between 15 K in the innermost areas of the filamentary cloud and 20 K in its outermost parts (Fig. 2). The temperatures from Eq. (3) were used to simulate the cloud images 𝒞_λ in all Herschel wavebands, assuming optically thin dust emission:

$$\begin{array}{c}\hfill {\mathcal{C}}_{\nu }=B_{\nu }({\mathcal{T}}_{\mathrm{C}}){\mathcal{D}}_{\mathrm{C}}{\kappa }_{\nu }\eta \mathrm{\mu }{m}_{\mathrm{H}},\end{array}$$(4)

where B_ν is the blackbody intensity, κ_ν is the dust opacity, η = 0.01 is the dust-to-gas mass ratio, μ = 2.8 is the mean molecular weight per H₂ molecule, and m_H is the hydrogen mass. The dust opacity was parameterized as a power law κ_ν = κ₀(ν/ν₀)^β with κ₀ = 9.31 cm² g⁻¹ (per gram of dust), λ₀ = 300 μm, and β = 2.

2.2. Simulated starless and protostellar cores

To populate the filamentary cloud with realistic sources, 156 radiative transfer models were computed by a numerical solution of the dust continuum radiative transfer problem in spherical geometry (using modust, Bouwman 2001). The models adopted tabulated absorption opacities κ_abs for dust grains with thin ice mantles (Ossenkopf & Henning 1994), corresponding to a density n_H = 10⁶ cm⁻³ and coagulation time t = 10⁵ yr. The opacity values at λ >  160 μm were replaced with a power law κ_λ ∝ λ⁻², consistent with the parameterization used in Eq. (4).

The models of three populations of starless cores and one population of protostellar cores cover wide ranges of masses (from 0.05 to 2 M_⊙) and half-maximum sizes (from ∼0.001 to 0.1 pc). Density profiles of the critical Bonnor–Ebert spheres were adopted for starless cores, whereas the protostellar cores have power-law densities ρ(r) ∝ r⁻². Starless cores consist of low-, medium-, and high-density subpopulations, following the M ∝ R relation for the isothermal Bonnor–Ebert spheres (with T_BE = 7,  14,  28 K) in the area of the mass–radius diagram occupied by prestellar cores observed in the Ophiuchus and Orion star-forming regions (Motte et al. 1998, 2001).

Both types of cores were embedded in background spherical clouds with a uniform surface density of 3 × 10²¹ cm⁻² and outer radius of 1.4 × 10⁵ AU (1000″ or 0.68 pc). In an isotropic interstellar radiation field (Black 1994) with the strength parameter G₀ = 10 (e.g., Parravano et al. 2003), the embedding clouds acquired temperatures of T ≈ 22 K at their edges, consistent with the highest values of 𝒯_C from Eq. (3). The embedding clouds lowered T(r) toward the interiors of both starless and protostellar cores. Accreting protostars in the centers of the protostellar cores, however, produced luminosity L_A ∝ M and thus sharply peaked temperature distributions deeper in their central parts.

2.3. Complete simulated images

Individual surface density images of the models of 828 starless and 91 protostellar cores were distributed in the dense areas (N_H2 ≥ 5 × 10²¹ cm⁻²) of the filamentary cloud 𝒟_C. They were added quasi-randomly, without overlapping, to the 𝒟_C image at positions, where their peak surface density exceeded that of the cloud N_H2 value. An initial mass function (IMF)-like power-law mass function with a slope dN/dM of −1.7 was used to determine the numbers of models per mass bin δ log₁₀​M  ≈  0.1 in each of the four populations. This resulted in the surface densities 𝒟_S, the intensities 𝒮_λ of sources (Fig. 3), and in the complete simulated images 𝒞_λ + 𝒮_λ.

The final simulated Herschel images ℋ_λ from Eq. (1) of the modeled star-forming region were obtained by adding different realizations of the random Gaussian noise 𝒩_λ at 70, 100, 160, 250, 350, and 500 μm and convolving the resulting images to the angular resolutions O_λ of 8.4, 9.4, 13.5, 18.2, 24.9, and 36.3″, respectively (Fig. 4). The resulting images ℋ_λ have σ noise levels of 6, 6, 5.5, 2.5, 1.2, and 0.5 MJy sr⁻¹, resembling the actual noise measured in the Herschel images of the Rosette molecular complex (Motte et al. 2010).

3. Source- and filament-extraction method

The main processing steps of getsf are outlined in Fig. 1, where several major blocks of the algorithm are highlighted. The method may be summarized as follows: (1) preparation of a complete set of images for an extraction, (2) separation of the structural components of sources and filaments from their backgrounds, (3) flattening of the residual noise and background fluctuations in the images of sources and filaments, (4) combination of the flattened components of sources and filaments over selected wavebands, (5) detection of sources and filaments in the combined images of the components, and (6) measurements of the properties of the detected sources and filaments.

Like its predecessors, getsf has just a single, user-definable parameter: the maximum size (width) of the structures of interest to extract. Internal parameters of getsf have been carefully calibrated and verified in numerous tests using large numbers of diverse images (both simulated and real-life observed images) to ensure that getsf works in all cases. This approach rests on the conviction that high-quality extraction methods for scientific applications must not depend on the human factor. It is the responsibility of the creator of a numerical method to make it as general as possible and to minimize the number of free parameters as much as possible. An internal multidimensional parameter space of complex numerical tools must never be delegated to the end user to explore if the aim is to obtain consistent and reliable scientific results.

3.1. Preparation of images for extraction

The multiwavelength extraction methods must be able to use all available information contained in the observed images across various wavebands with different angular resolutions. It is usually beneficial to collect all available images for a specific region of the sky under study.

3.1.1. Original observed set of images

To prepare multiwavelength ℋ_λ for processing with getsf, it is necessary to convert them into the images ℐ_λ, all on the same grid of pixels. To this end, getsf resamples all images (using swarp, Bertin et al. 2002) on a pixel size, chosen to be optimal for the highest-resolution images available. It is very important to carefully verify alignment of the resampled ℐ_λ and correct it (if necessary) to ensure that all unresolved intensity peaks remain on the same pixel across all wavebands. To reveal possible misalignments, it is sufficient to open each pair of prepared images in ds9 (Joye & Mandel 2003) and blink the two frames, going from the highest-resolution to the lowest-resolution images.

Most astronomical images have irregularly shaped coverage and limited usable areas that differ between wavebands. To include only the “good” parts of the ℐ_λ coverage in the image processing, it is necessary to create masks ℳ_λ (with pixel values 1 or 0). With these masks, getsf can process only the good areas of ℐ_λ that have a mask value of 1. To facilitate the image preparation, getsf always creates default masks ℳ_λ = 1. However, for most real observations, the masks must be prepared very carefully and independently for each image. To manually create the masks, one can use imagej (Abràmoff et al. 2004) or gimp² that allows users to create a polygon over an image, convert the polygon into a mask, and save it in the FITS format.

3.1.2. Derived high-resolution images

The multiwavelength far-infrared Herschel images open the possibility of computing maps of surface density and dust temperature by fitting the spectral shapes Π_λ of the image pixels. The standard procedure assumes that (1) the original images represent optically thin thermal emission of dust grains with a power-law opacity κ_λ ∝ λ^−β and a constant β value, (2) the dust temperature is constant along the lines of sight passing through each pixel of the images, and (3) the lines of sight are not contaminated by unrelated radiation at either end, in front of the observed structures and behind them. Unfortunately, one or more of the assumptions are likely to be invalid, especially the stipulation of the opacity law and the constant line-of-sight temperatures (e.g., Men’shchikov 2016). The values of the derived surface densities and temperatures therefore must be considered as fairly unreliable and implying large error bars.

When we assume that the observations include the Herschel images, the spectral shapes Π_λ of each pixel can be fit at several wavelengths (160−500 μm) and resolutions (18.2−36.3″), which results in three sets of surface densities and dust temperatures. The highest-resolution derived images are the least reliable because they are obtained from fitting only two images (at 160 and 250 μm), whereas the lowest-resolution maps are the most accurate because they come from fitting four independent images (at 160, 250, 350, and 500 μm).

In an attempt to combine the higher accuracy of the lower-resolution images with the higher angular resolutions of the less accurate images, Palmeirim et al. (2013) published a simple algorithm that uses complementary spatial information contained in the observed images to create a surface density image with the resolution O_P = 18.2″ of the 250 μm image. When this approach is extended to temperatures, the sharper images can be computed by adding the higher-resolution information to the low-resolution images as differential terms,

$$\begin{array}{c}\hfill {\{D|T\}}_{\mathrm{P}}={\{D|T\}}_4+\delta {\{D|T\}}_3+\delta {\{D|T\}}_2,\end{array}$$(5)

where the base surface density and temperature {𝒟|𝒯}₄ are derived by fitting the 160, 250, 350, and 500 μm images at the lowest resolution O₅₀₀ = 36.3″. The additional terms, containing the higher-resolution contributions, are produced by unsharp masking,

$$\begin{array}{c}\hfill \delta {\{D|T\}}_{\{2|3\}}={\{\mathcal{D}|\mathcal{T}\}}_{\{\mathcal{2}|\mathcal{3}\}}-{\mathcal{G}}_{\{\mathcal{3}|\mathcal{4}\}}\ast {\{\mathcal{D}|\mathcal{T}\}}_{\{\mathcal{2}|\mathcal{3}\}},\end{array}$$(6)

where {𝒟|𝒯}₃ are computed by fitting the three images at 160, 250, and 350 μm at the resolution O₃₅₀ = 24.9″, and {𝒟|𝒯}₂ are obtained by fitting the two images at 160 and 250 μm at the resolution O₂₅₀ = 18.2″; the Gaussian kernels 𝒢_{3|4} convolve the images to the next lower resolutions O_{350|500}.

The following generalization of the above algorithm allows deriving surface densities and temperatures with any (arbitrarily high) angular resolution existing among the observed ℐ_λ. The three independently derived maps of temperatures 𝒯_{2|3|4} with the resolutions of 18.2−36.3″ and six observed Herschel images with their native resolutions O_λ of 8.4−36.3″ define 18 surface densities,

$$\begin{array}{c}\hfill {\mathcal{D}}_{O_{\lambda }\{2|3|4\}}=\frac{{\mathcal{I}}_{\nu }}{B_{\nu }({\mathcal{T}}_{\{2|3|4\}}){\kappa }_{\nu }\eta \mathrm{\mu }{m}_{\mathrm{H}}},\end{array}$$(7)

with the assumptions and parameterizations of Eq. (4). It is required that the resolution of temperatures must not be higher than O_λ, which excludes 𝒟_O3502 and 𝒟_O500{2|3} from the algorithm and provides 15 independently computed variants of the surface densities of the observed region, with different resolutions. The high-resolution surface density image is computed as

$$\begin{array}{c}\hfill {\mathcal{D}}_{O_{\mathrm{H}}}={\mathcal{D}}_{O_{500}}+{\displaystyle \sum _{\lambda ={\lambda }_{\mathrm{H}}}^{500}}max(\delta {\mathcal{D}}_{O_{\lambda }2},\delta {\mathcal{D}}_{O_{\lambda }3},\delta {\mathcal{D}}_{O_{\lambda }4}),\end{array}$$(8)

where λ_H denotes the wavelength of the image ℐ_​λH with the desired angular resolution O_H ≡ O_λH and the differential terms with higher-resolution information are obtained by the same unsharp masking,

$$\begin{array}{c}\hfill \delta {\mathcal{D}}_{O_{\lambda }\{2|3|4\}}={\mathcal{D}}_{O_{\lambda }\{2|3|4\}}-{\mathcal{G}}_{O_{\lambda +}}\ast {\mathcal{D}}_{O_{\lambda }\{2|3|4\}},\end{array}$$(9)

where 𝒢_Oλ+ is the Gaussian kernel (regarded as the delta function at 500 μm), convolving 𝒟_Oλ{2|3|4} to a lower resolution of the next longer wavelength. For images at λ <  250 μm, only the positive values of δ𝒟_Oλ{2|3|4} are used in Eq. (8) to circumvent the problem of creating artificial depressions and negative pixels around strong peaks due to the resolution mismatch (O_λ <  O₂₅₀) between ℐ_​λ and the lower-resolution 𝒯_{2|3|4} in Eq. (7).

The problem is caused by the sharp radial temperature gradients toward the unresolved centers of protostellar cores (cf. Men’shchikov 2016). They are smeared out by the low resolutions of 18.2−36.3″, hence the fitting of Π_λ leads to underestimated temperatures 𝒯_{2|3|4} and overestimated values (within an order of magnitude) of peak surface densities at higher resolutions O_λ <  O₂₅₀. This means that unsharp masking of the overestimated peaks could create negative annuli in δ𝒟_Oλ{2|3|4} and negative pixels in 𝒟_OH. Fortunately, the surface densities are quite accurate outside the unresolved peaks (Appendix A).

A slight modification of Eq. (8) allows deriving the high-resolution surface densities with an enhanced contrast of all unresolved or slightly resolved structures,

$$\begin{array}{c}\hfill {\mathcal{D}}_{O_{\mathrm{H}}}^{+}={\mathcal{D}}_{O_{\mathrm{H}}}+{\displaystyle \sum _{\lambda ={\lambda }_{\mathrm{H}}}^{500}}{\displaystyle \sum _{n=2}^4}max(\delta {\mathcal{D}}_{O_{\lambda }n},0),\end{array}$$(10)

where the positive parts of the differential high-resolution terms from Eq. (9) are added to 𝒟_OH. These high-contrast images may be useful for detection of unresolved structures because the latter are usually diluted by the observations with insufficient angular resolution; a higher contrast improves their visibility.

A high-resolution temperature 𝒯_​OH, consistent with the high-resolution surface density 𝒟_OH, is computed by numerically inverting the Planck function,

$$\begin{array}{c}\hfill {\mathcal{T}}_{O_{\mathrm{H}}}=B_{{\nu }_{\mathrm{H}}}^{-1}\left(\frac{{\mathcal{I}}_{{\nu }_{\mathrm{H}}}}{{\mathcal{D}}_{O_{\mathrm{H}}}{\kappa }_{{\nu }_{\mathrm{H}}}\eta \mathrm{\mu }{m}_{\mathrm{H}}}\right),\end{array}$$(11)

with ${\nu }_{\mathrm{H}}=c{\lambda }_{\mathrm{H}}^{-1}$, where c is the speed of light. The high-resolution images {𝒟|𝒯}_13″ are shown in Fig. 5 along with the true simulated 𝒟_C + 𝒟_S (Sect. 2.3). A comparison demonstrates that the pixel-fitting procedure reduces visibility of many unresolved or slightly resolved starless cores, which is the manifestation of the invalid assumption of the uniform line-of-sight temperatures. Starless (prestellar) cores have lower temperatures in their centers, and their smearing by an insufficient resolution leads to overestimated temperatures that suppress the surface density peaks (cf. Fig. A.1).

The new hires algorithm, outlined by Eqs. (7)–(11), brings the benefits of a resolution O_H ≈ 8″, twice better than O_P and four times better than O₅₀₀, if the image quality at the shortest wavelengths permits this. Moreover, the angular resolutions of the Herschel images at 70, 100, and 160 μm, obtained with a slow scanning speed of 20″ s⁻¹, are even higher: 6, 7, and 11″, respectively. These observations, illustrated in Fig. 6, allow deriving the surface densities and temperatures with O_H ≈ 6″, a three times better resolution than when using Eq. (5). If the 70 μm image is too noisy or there is evidence of its strong contamination by emission unrelated to that of the adopted dust grains (e.g., polycyclic aromatic hydrocarbons or transiently heated very small dust grains), then the derived images may still have the 7−13″ resolution of the 100 or 160 μm wavebands. In addition to the high resolution, the images from Eq. (8) also have a better quality than those from Eq. (5) because they accumulate all available high-resolution information from the (up to 15) independently computed images 𝒟_Oλ{2|3|4} that use all three temperatures 𝒯_{2|3|4} with each original ℐ_​λ.

The hires algorithm works with any number 2 ≤ N ≤ 6 of Herschel wavebands. If the 160 μm image is unavailable or disabled, then the temperature 𝒯₂ at the resolution O₂₅₀ is removed from Eq. (7) and 𝒯_{3|4} at the resolutions O_{350|500} are obtained from fitting of only the 250, 350, and 500 μm images. If the 250 μm image is also unavailable or disabled, then only the single temperature 𝒯₄ at the lowest resolution O₅₀₀ remains in Eq. (7), obtained from the 350 and 500 μm images. Although the algorithm is unaffected by the changes, the reduction in the number of independently derived temperatures would lower the angular resolution and accuracy of the resulting surface density image. The improved algorithm can use the realistic, non-Gaussian point-spread functions (PSF) published by Aniano et al. (2011). However, the surface densities are largely determined by the SPIRE bands with nearly Gaussian PSFs, whereas only the PACS 160 μm band is used in the pixel fitting. The benchmark tests have shown that effects of the realistic PSFs on surface densities are very small, at percent levels, much smaller than the general uncertainties of the pixel-fitting methods (Appendix A). It is therefore sufficient to use the Gaussian PSFs when surface densities are derived.

For some studies, it might be useful to have all images at the same wavelength-independent angular resolution. With the high-resolution surface densities and temperatures from Eqs. (8) and (11), it is straightforward to obtain such images:

$$\begin{array}{c}\hfill {\mathcal{J}}_{\nu O_{\mathrm{H}}}=B_{\nu }({\mathcal{T}}_{O_{\mathrm{H}}}){\mathcal{D}}_{O_{\mathrm{H}}}{\kappa }_{\nu }\eta \mathrm{\mu }{m}_{\mathrm{H}},\end{array}$$(12)

with the assumptions and parameterizations of Eq. (4). For example, the intensities 𝒥_​λ13″ at 250, 350, and 500 μm would be sharper than ℐ_​λ by the factors 1.3, 1.8, and 2.7, respectively.

When the available original set of images ℐ_​λ allows creation of 𝒟_OH, it is advantageous to have it complement the original data set, handling it as an image ℐ_ƛ “observed” in a fictitious waveband ƛ. In the multiwavelength extractions with getsf, it may be recommended to use 𝒟_OH for better detections and deblending of dense structures. The surface densities are not accurate enough for source measurements, as demonstrated in Appendix A and Men’shchikov (2016).

The following presentation and discussion of getsf implicitly assumes that the additional detection images are contained in the set of images ℐ_​λ. In other words, all supplementary wavebands are included in the set of λ prepared for extraction. The latter was done for a multiwavelength data set that included all images in the Herschel wavebands (Fig. 4) and the high-resolution surface density ℐ_ƛ ≡ 𝒟_13″ (Fig. 5), a total of N_W = 7 wavelengths.

3.1.3. Practical definition of maximum size

Before starting any extraction with getsf, it is necessary to formulate the aim of the study and determine what structures of interest are to be extracted. The method knows and is able to separate three types of structures: sources, filaments, and backgrounds. To separate the structural components with getsf, the maximum size {X|Y}_λ of the sources (X_λ) and filaments (Y_λ) of interest needs to be manually (visually) estimated from the prepared ℐ_​λ independently for each waveband, which can be accomplished by opening an image in ds9 and placing a circular region fully covering the width of the largest structure. The maximum size of structures is the single physical parameter that the method needs to know for each observed image. Being a function of the type of structures (sources, filaments) and the waveband λ, it is split into X_λ and Y_λ in this paper for convenience.

The maximum size {X|Y}_λ is defined as the footprint radius (in arcsec) of the largest source and the widest filament to be extracted. A footprint size has the meaning of a full width at zero (background) level: the largest two-sided extent from a source peak or filament crest at which this structure is still visible in ℐ_​λ against its background. For a Gaussian intensity distribution, the footprint radius is slightly larger than the half-maximum width H_λ of a structure. For a power-law intensity profile, the footprint radius may become much larger than H_λ. If the widest filaments of interest are blended (overlapping each other with their footprints), Y_λ must be increased accordingly to approximate the full extent of the blend. In contrast, it is not necessary to adjust X_λ for blended sources because their final background will be determined from their footprints at the measurement step (Sect. 3.4.6).

It is not necessary (also not possible) to evaluate the maximum size parameter {X|Y}_λ very precisely, a 50% accuracy is quite sufficient. Its purpose is to set a reasonable limit to the spatial scales when separating the structural components, and to the size of the structures to be measured and cataloged. The method works with spatially decomposed images, and it needs to know the maximum scale. It makes no sense to perform the decomposition up to a very large scale if the extraction is aimed at much smaller sources or narrower filaments. The method has no limitations with respect to the sizes (widths) of the structures to extract. However, it is important to avoid detecting, measuring, and cataloging the peaks that are unnecessarily too wide because they would likely overlap with other sources of interest, which potentially would make their measurements less accurate.

To extract all structures in the benchmark images presented in Sect. 2, the estimated X_λ values for sources are 16, 25, 30, 150, 150, and 150″, whereas the estimated Y_λ values for the filament are 350″ in all six Herschel wavebands (Fig. 4); in the additional surface density image ℐ_ƛ ≡ 𝒟_13″ (Fig. 5), the {X|Y}_ƛ values are the same as those for the 250−500 μm images.

3.2. Backgrounds of the structural components

Complex observed images may be radically simplified by subtracting backgrounds on spatial scales much larger than the maximum size {X|Y}_λ. The independent largest sizes for sources and filaments effectively define two different backgrounds for the two scales. The X_λ-scale background ℬ_λX is derived to separate the component of sources 𝒮_λ, whereas the Y_λ-scale background ℬ_λY is obtained to separate the component of filaments ℱ_λ. The backgrounds are collectively referred to as ℬ_λ{X|Y}.

The true background under the observed structures is fundamentally unknown, and it is a major source of large uncertainties and measurement inaccuracies, especially for the faintest structures. In practice, the backgrounds ℬ_λ{X|Y} are defined in getsf as the smooth intensity distributions on spatial scales S_​j larger than 4{X|Y}_λ that remain in ℐ_​λ after a complete removal of all sources or filaments with the maximum size of {X|Y}_λ. In contrast to the background derivation by median filtering (getimages, Paper III), which may become extremely slow for very large images and wide structures, getsf employs a more direct, precise, and effective clipping algorithm to separate the structures.

3.2.1. Decomposition of the original images

In general, observed images are very complex blends of various structural components on different spatial scales, and great advantages are obtained when a spatial decomposition is used to simplify the images (cf. Papers I and II). Following the getold approach, getsf employs successive unsharp masking (Appendix B) to decompose the original images ℐ_​λ into a set of single-scale images ℐ_​λj (Fig. 7). It also uses an iterative algorithm (Appendix B) to determine a single-scale standard deviation σ_​λj, as well as its total value σ_​λ, which are used to separate the structural components present in ℐ_​λ.

3.2.2. Separation of the structural components

The backgrounds ℬ_λ{X|Y} are computed by cutting small round peaks and elongated structures off the decomposed images ℐ_​λj and recovering the full images using Eq. (B.4). It is important to note that the appearance of the structures in the decomposed images depends on both the spatial scale and intensity level.

To remove the structural components, getsf slices ℐ_​λj by a number N_L of intensity levels I_λjl, spaced by δlnI_λj = 0.05 from the image maximum down to σ_​λj for sources and to 0.3σ_​λj for filaments. Each slice l cuts through all the structures present in ℐ_​λj on that intensity level, producing various shapes of connected pixels,

$$\begin{array}{c}\hfill {\mathcal{I}}_{\lambda jl}=min(max({\mathcal{I}}_{\lambda j},I_{\lambda jl}),I_{\lambda jl}),l=1,2,\cdots ,N_{\mathrm{L}}.\end{array}$$(13)

Relatively round source-like peaks in ℐ_​λj may be effectively distinguished from elongated structures by the number of connected pixels N_λjl that their shapes occupy in the slice ℐ_​λjl (cf. Papers I and II). The single-scale images indeed most clearly show the structures with matching sizes (H_λ ≈ S_​j), whereas the signals from much narrower and much wider structures are suppressed. As a consequence, the source-like shapes occupy relatively small areas of connected pixels in ℐ_​λjl that are comparable to the area $\pi S_j^2$ of the convolution kernel 𝒢_​j. In contrast to the round peaks, elongated shapes in ℐ_​λjl have greater lengths L_λ than widths W_λ ≈ S_​j, which means that the filamentary shapes in slices ℐ_​λjl extend over much larger areas than $\pi S_j^2$.

In addition to N_λjl, getsf uses two more quantities to discriminate between sources and filaments: elongation E_λjl and sparsity S_​λjl. They are defined by the major and minor sizes (a_λjl and b_λjl) of each cluster of connected pixels, obtained from intensity moments (cf. Appendix F in Paper I),

$$\begin{array}{c}\hfill E_{\lambda jl}\equiv \frac{a_{\lambda jl}}{b_{\lambda jl}},S_{\lambda jl}\equiv \frac{\pi a_{\lambda jl}{b}_{\lambda jl}}{N_{\lambda jl}\mathrm{\Delta }},\end{array}$$(14)

where Δ is the pixel size. Only simple and relatively straight filamentary shapes can be identified in ℐ_​λjl by their elongation. Most of the actually observed filaments in space are shaped quite irregularly on different scales and intensity levels. The elongation E_λjl alone cannot be used to quantify strongly curved, not very dense clusters of connected pixels that meander around (e.g., a spiral structure). Although E_λjl may well be close to unity for sparse shapes, high values of S_​λjl for these structures would indicate that they do not belong to sources.

The structural components are separated in single scales ℐ_​λj using the three quantities described above. The shapes produced by sources in a slice ℐ_​λjl are not very elongated, not very sparse, and not very large. In contrast, the shapes produced by filaments in a slice ℐ_​λjl are elongated or sparse. Hence, these definitions for the source-like and filament-like shapes are written as

$$\begin{array}{c}\hfill \begin{array}{c}{E}_{\lambda jl}\le 1.47,S_{\lambda jl}\le 1.39,N_{\lambda jl}\le \pi {\left({\xi }_{\lambda j}{S}_j\right)}^2{\mathrm{\Delta }}^{-2},\hfill \\ E_{\lambda jl}>3.00\vee S_{\lambda jl}>1.39,\hfill \end{array}\end{array}$$(15)

where the limiting values of elongation and sparsity were determined empirically from numerous benchmark extractions. The ξ_λj factor accounts for the fact that the area of a decomposed unresolved peak increases nonlinearly toward the smallest spatial scales S_​j ≲ O_λ. The factor may be determined empirically by decomposing an unresolved peak 𝒫 in single scales 𝒫_​j (Fig. B.1) and finding the distances θ, where the one-dimensional profile P_​j(θ) through the peak has dP_​j/dθ = 0 for P_​j <  0,

$$\begin{array}{c}\hfill {\xi }_{\lambda j}=0.47{\left(O_{\lambda }{S}_j^{-1}\right)}^{1.34}+0.83.\end{array}$$(16)

The ξ_λj factor ensures that N_λjl has appropriate values and that single-scale peaks are clipped cleanly on all spatial scales.

Various shapes formed by connected pixels are identified and analyzed in each single-scale slice using the tintfill algorithm (Smith 1979)³, previously employed by getold to detect sources and filaments (Papers I and II). Deriving the background ℬ_λX of sources, getsf decomposes ℐ_​λ and removes all source-like shapes from ℐ_​λjl, according to their definition in Eq. (15), in an iterative procedure (Sect. 3.2.3). Deriving the background ℬ_λY of filaments, getsf decomposes ℬ_λX and removes all filament-like shapes from ℬ_λXjl, according to their definitions in Eq. (15), in the same iterative procedure. The shapes are erased from each slice l by setting all their pixels to zero.

3.2.3. Reconstruction of the backgrounds

When we denote with ℬ_λ{X|Y}jlC either of the single-scale background slices ℬ_λXjlC or ℬ_λYjlC after the shape removal, the backgrounds on scale j are reassembled from the clipped slices as

$$\begin{array}{c}\hfill {\mathcal{B}}_{\lambda \{X|Y\}j\mathrm{C}}={\displaystyle \sum _{l=1}^{N_{\mathrm{L}}}}{\mathcal{B}}_{\lambda \{X|Y\}jl\mathrm{C}}.\end{array}$$(17)

To properly reconstruct the complete backgrounds ℬ_λ{X|Y} from ℬ_λ{X|Y}jC, it is not sufficient to just sum them over scales. The single-scale processing scheme requires that it must be done indirectly in several steps by reconstructing the complete images of sources and filaments.

In the first step, getsf recomputes the single-scale sources and filaments that have been clipped, removing all negative values from the reassembled single-scale backgrounds,

$$\begin{array}{c}\hfill {\{\mathcal{S}|\mathcal{F}\}}_{\lambda j}=\{{\mathcal{I}}_{\lambda }|{\mathcal{B}}_{\lambda X}\}-max({\mathcal{B}}_{\lambda \{X|Y\}j\mathrm{C}},0).\end{array}$$(18)

In the second step, getsf computes the full images of the sources and filaments over all scales, recursively summing the clipped structures from the largest to the smallest scales and removing all negative values from each partial sum,

$$\begin{array}{c}\hfill {\{\mathcal{S}|\mathcal{F}\}}_{\lambda }=max({\{\mathcal{S}|\mathcal{F}\}}_{\lambda }+{\{\mathcal{S}|\mathcal{F}\}}_{\lambda j},0),j=J_{\lambda \{X|Y\}},\cdots ,2,1,\end{array}$$(19)

where J_λ{X|Y} is the number of the largest spatial scales 4{X|Y}_λ for the backgrounds ℬ_λ{X|Y} and the initial value of the recursive sum is set to zero. The complete backgrounds are obtained by subtracting the structures from the original images,

$$\begin{array}{c}\hfill {\mathcal{B}}_{\lambda \{X|Y\}}^0=\{{\mathcal{I}}_{\lambda }|{\mathcal{B}}_{\lambda X}\}-{\{\mathcal{S}|\mathcal{F}\}}_{\lambda }.\end{array}$$(20)

The initial backgrounds in Eq. (20) are only the first approximations because they contain substantial residual contributions from the original structures. It is straightforward to define iterations to improve the backgrounds by decomposing them and clipping the residual shapes from each single scale. The algorithm described by Eqs. (17)–(20) remains the same, with two substitutions,

$$\begin{array}{c}\hfill {\{{\mathcal{I}}_{\lambda }|{\mathcal{B}}_{\lambda X}\}}^i\leftarrow {\mathcal{B}}_{\lambda \{X|Y\}}^{i-1},i=1,2,\cdots ,N_{\mathrm{I}},\end{array}$$(21)

where N_I is the number of iterations. Each successive iteration reduces contributions of the residual structures and improves the backgrounds until corrections in all pixels become small compared to the originals,

$$\begin{array}{c}\hfill \delta {\mathcal{B}}_{\lambda \{X|Y\}}^i<0.003(\{{\mathcal{I}}_{\lambda }|{\mathcal{B}}_{\lambda X}\}+10{\sigma }_{\lambda }),\end{array}$$(22)

where the additional term helps avoid unnecessary iterations in rare cases when the images contain extremely faint pixels.

The final background-subtracted structural components are computed as

$$\begin{array}{c}\hfill {\{\mathcal{S}|\mathcal{F}\}}_{\lambda }=\{{\mathcal{I}}_{\lambda }|{\mathcal{B}}_{\lambda X}\}-{\mathcal{B}}_{\lambda \{X|Y\}}.\end{array}$$(23)

The original images can be recovered by summing the three separated components: ℐ_​λ = 𝒮_λ + ℱ_λ + ℬ_λY. The positive parts of the small-scale background fluctuations and instrumental noise are contained in the component 𝒮_λ, hence the component ℱ_λ appears fairly smooth (Fig. 8).

3.3. Flattening of the structural components

Observations demonstrate that the levels of the large-scale backgrounds and their smaller-scale fluctuations often differ by orders of magnitude in various parts of large images. Although the subtraction of the smooth backgrounds ℬ_λ{X|Y} greatly simplifies the original images, it does not reduce the strong variations of the smaller-scale fluctuation levels across {𝒮|ℱ}_λ. As a consequence, many structures detected with global thresholds in the areas of stronger fluctuations may actually be spurious and unrelated to any real physical objects. On the other hand, faint real structures in the image areas with the lower levels of fluctuations may escape detection because the global threshold value is likely to be overestimated for those areas. To produce complete and reliable extractions using constant thresholds, it is necessary to make small-scale fluctuations uniform over the entire background-subtracted images.

The fluctuation levels are equalized using flattening images 𝒬_λ and ℛ_λ that are derived by getsf from the images 𝒰_λ and 𝒱_​λ of the standard deviations computed in the structural components with a circular sliding window of a radius O_λ,

$$\begin{array}{c}\hfill {\{\mathcal{U}|\mathcal{V}\}}_{\lambda }={\mathrm{sd}}_{O_{\lambda }}\left({\{\mathcal{S}|\mathcal{F}\}}_{\lambda \mathrm{R}}\right),\end{array}$$(24)

where 𝒮_λR and ℱ_λR are the regularized images 𝒮_λ and ℱ_λ, obtained using a smoother version of their backgrounds that is median-filtered using a sliding window of a radius 2O_λ and convolved with a Gaussian kernel of a half-maximum size O_λ,

$$\begin{array}{c}\hfill {\{\mathcal{S}|\mathcal{F}\}}_{\lambda R}=\{{\mathcal{I}}_{\lambda }|{\mathcal{B}}_{\lambda X}\}-{\mathcal{G}}_{O_{\lambda }}\ast {\mathrm{mf}}_{2O_{\lambda }}\left({\mathcal{B}}_{\lambda \{X|Y\}}\right).\end{array}$$(25)

This is done to improve the quality of {𝒰|𝒱}_λ for further processing because the structural components from by Eq. (23) are positively defined and have large areas of zero pixels. The regularized components in Eq. (25) acquire small-scale fluctuations resembling the background and noise fluctuations of the original images.

3.3.1. Decomposition of the standard deviations

The advantages of the spatial decomposition (Appendix B) apply also to the standard deviations {𝒰|𝒱}_λ. The getsf method produces the single-scales {𝒰|𝒱}_λj and employs the same iterative algorithm (Appendix B) to determine the single-scale standard deviation σ_​λj and its total value σ_​λ. This is done using the same procedure as was applied to ℐ_​λ in Sect. 3.2.1.

3.3.2. Removal of the structural features

The {𝒰|𝒱}_λ images sample local fluctuations and intensity gradients, revealing all sources and filaments present in ℐ_​λ (Figs. 9 and 10). To produce the corresponding flattening images, it is necessary to remove all such features from {𝒰|𝒱}_λ, hence to determine their {X|Y}_λ-scale backgrounds. Deriving the latter, getsf creates single-scale slices {𝒰|𝒱}_λjl, in a complete analogy with I_λjl in Sect. 3.2.2, and clips from them all source- and filament-like shapes according to their definitions in Eq. (15). The reconstructed backgrounds 𝒬_λX and ℛ_λY are computed using the iterative algorithm described in Sect. 3.2.3, with the largest spatial scale set to 2.5{X|Y}_λ.

When the background iterations converge, numerous sharp craters remain in the derived backgrounds {𝒬|ℛ}_λ{X|Y} that could create spurious structures if the images were used to flatten the structural components. To avoid this, the final flattening images 𝒬_λ and ℛ_λ (Figs. 9 and 10) are obtained by median filtering the background in circular sliding windows of radii 2O_λ and 5O_λ, respectively,

$$\begin{array}{c}\hfill {\{\mathcal{Q}|\mathcal{R}\}}_{\lambda }={\mathrm{mf}}_{\{2|5\}{O}_{\lambda }}\left({\{\mathcal{Q}|\mathcal{R}\}}_{\lambda \{X|Y\}}\right).\end{array}$$(26)

This important step ensures that flattening would never produce spurious structures in the detection images.

3.3.3. Flattening of the detection images

The detection images of the separated structural components are used to identify peaks of the sources and skeletons of the filaments, respectively (Sect. 3.4). Both source- and filament-detection images are flattened, that is, divided by the flattening images,

$$\begin{array}{c}\hfill {\{\mathcal{S}|\mathcal{F}\}}_{\lambda \mathrm{D}}=\frac{{\{\mathcal{S}|\mathcal{F}\}}_{\lambda }}{{\{\mathcal{Q}|\mathcal{R}\}}_{\lambda }}.\end{array}$$(27)

The standard deviations ${\mathrm{sd}}_{O_{\lambda }}({\{\mathcal{S}|\mathcal{F}\}}_{\lambda \mathrm{R}}{\{\mathcal{Q}|\mathcal{R}\}}_{\lambda }^{-1})$ in the regularized flattened components demonstrate that the detection images are remarkably flat outside the structures, as shown in Figs. 9 and 10. This ensures an accurate separation of significant structures from the fainter background and noise fluctuations during the subsequent extraction of sources and filaments.

3.4. Extraction of the structural components

To extract sources and filaments means to detect them and measure their properties. The background subtraction and flattening algorithms presented in Sects. 3.2 and 3.3 radically simplify the originals ℐ_​λ, separating two distinct structural components and creating the independent flat detection images {𝒮|ℱ}_λD. In contrast to the originals, the flat images are suitable for the detection techniques that apply a threshold value for the entire image.

3.4.1. Decomposition of the detection images

In order to accurately extract various structures that widely range in brightness and size, it is essential to use the benefits offered by the single-scale spatial decomposition (Appendix B). Following its general approach (Fig. 1, Sects. 3.2.1 and 3.3.1), getsf decomposes the detection images 𝒮_λD and ℱ_λD into single scales {𝒮|ℱ}_λDj and estimates the corresponding standard deviations σ_​λSj and σ_​λFj (Appendix B) that are necessary for separating significant structures from all other fluctuations. The decomposed components 𝒮_λDj and ℱ_λDj are shown in Figs. 11 and 12 after they were cleaned and combined over wavebands.

3.4.2. Cleaning of the single-scale detection images

Cleaning is the removal of insignificant background and noise fluctuations from detection images that needs to be done before combining them over wavebands (Sect. 3.4.3). The clean images of the structures are obtained by preserving only the pixels with values above the cleaning thresholds ϖ_λ{S|F}j and by setting all fainter pixels to zero,

$$\begin{array}{c}\hfill {\{\mathcal{S}|\mathcal{F}\}}_{\lambda \mathrm{D}j\mathrm{C}}=max({\{\mathcal{S}|\mathcal{F}\}}_{\lambda \mathrm{D}j},{\varpi }_{\lambda \{\mathrm{S}|\mathrm{F}\}j}),\end{array}$$(28)

where ϖ_λSj = 5σ_​λSj and ϖ_λFj = 2σ_​λFj. The filament threshold is significantly lower than that for sources because getsf additionally cleans ℱ_λDjC of the residual source-like clusters of connected pixels according to their definition in Eq. (15).

The resulting clean images {𝒮|ℱ}_λDjC (Figs. 11 and 12) are deemed to have signals only from the sources and filaments, respectively. In practice, some of them may have several faint spurious peaks and other structures that are discarded during the subsequent detection and measurement steps.

3.4.3. Combination of the clean single scales over λ

All previous image processing was done independently for each wavelength. It is recommended to always process the input images in parallel, independent getsf runs to reduce the total extraction time approximately by a factor N_W, the number of wavelengths. Now, getsf accumulates the clean single-scale images 𝒮_λDjC and ℱ_λDjC over the wavebands in order to use the independent information from all images and enhance the signal of the significant structures. This procedure follows the getold approach (Paper I), with the important improvement that filaments are handled in the same way as sources. When decomposed images on each scale are combined, differences in the angular resolutions between the wavebands are much less important because the single-scale images select and enhance the structures with widths similar to the scale size S_​j, not the resolution O_λ.

The clean detection images are normalized before their accumulation over wavelengths to make all cleaning thresholds equal to unity in all bands. The combination process is described by the following expression:

$$\begin{array}{c}\hfill {\{\mathcal{S}|\mathcal{F}\}}_{\mathrm{D}j\mathrm{C}}=N_{\{\mathrm{S}|\mathrm{F}\}}^{-1}{\displaystyle \sum _{\lambda }}{f}_{\lambda j}max({\{\mathcal{S}|\mathcal{F}\}}_{\lambda \mathrm{D}j\mathrm{C}},{\mathcal{Z}}_{\lambda \{\mathrm{S}|\mathrm{F}\}j}){\varpi }_{\lambda \{\mathrm{S}|\mathrm{F}\}j}^{-1},\end{array}$$(29)

where N{S|F} ≤ N_W is the number of the wavebands chosen to be used in the combination, 𝒵_λ{S|F}j is the threshold image (equal to ϖ_λ{S|F}jŌ in all pixels), and f_λj is a factor that gradually turns the smallest scales on,

$$\begin{array}{c}\hfill f_{\lambda j}=min({\left(S_j{O}_{\lambda }^{-1}\right)}^3,1).\end{array}$$(30)

This factor ensures that the small-scale noise or artifacts appearing on top of the resolved structures do not produce spurious detections in the combined images on scales S_​j <  O_λ. Sufficiently bright unresolved structures still contribute to {𝒮|ℱ}_DjC on the smallest scales below O_λ. This super-resolution is useful to detect blended unresolved peaks. Selected combined images of the two structural components are shown in Figs. 11 and 12.

The normalization to a common threshold in Eq. (29) is a natural way of maximizing sensitivity of the combined images. This procedure modifies the original dependence of the source brightness on spatial scales, however, which is analyzed by the detection algorithm (Sect. 3.4.4) to determine the characteristic size for each source. Therefore a second set of combined images is defined for the component of sources, normalized to the smallest scale in each waveband,

$$\begin{array}{c}\hfill {\stackrel{\sim }{\mathcal{S}}}_{\mathrm{D}j\mathrm{C}}={\displaystyle \sum _{\lambda }}\frac{w_{\lambda }}{{\mathcal{S}}_{\lambda \mathrm{D}1\mathrm{C}}}{\mathcal{S}}_{\lambda \mathrm{D}j\mathrm{C}},\end{array}$$(31)

where w_λ is the weight that enhances the contribution of the images with higher angular resolutions,

$$\begin{array}{c}\hfill w_{\lambda }={\left(\frac{\overline{O}}{O_{\lambda }}\right)}^7,\overline{O}=N_{\mathrm{W}}^{-1}{\displaystyle \sum _{\lambda }}{O}_{\lambda },\end{array}$$(32)

where Ō is the average resolution, and the power of 7 ensures complete separation of the contributions of different wavebands in Eq. (31). After the weighting, the summation of 𝒮_λDjC preserves the individual dependence of the peak intensity of each source on spatial scales, which provides an initial estimate of its size during detection before the actual measurements.

3.4.4. Detection of sources in the combined images

Sources are detected in 𝒮_DjC with almost the same algorithm that was used by getold (Sect. 2.5 of Paper I), which is briefly summarized here for completeness. An inspection of the entire set of single-scale images 𝒮_DjC shows that sources appear on relatively small scales become the brightest on scales roughly equal to their size and vanish on significantly larger scales (cf. Fig. B.1). All detectable sources appear isolated on small scales and become blended with other nearby sources on larger scales. The getsf source detection scheme identifies the sources in 𝒮_DjC and tracks their evolution from small to large scales, until they disappear or merge with a nearby brighter source.

To detect sources, getsf slices 𝒮_DjC by a number N_L of intensity levels I_jl, spaced by δlnI_j = 0.01, from the image maximum down to the lowest non-zero value. Each slice l cuts through all peaks brighter than I_jl, producing a set of partial images,

$$\begin{array}{c}\hfill {\mathcal{S}}_{\mathrm{D}j\mathrm{C}l}=max({\mathcal{S}}_{\mathrm{D}j\mathrm{C}},I_{\mathit{jl}}),l=1,2,\cdots ,N_{\mathrm{L}}.\end{array}$$(33)

The source detection algorithm works on the sequence of partial images, creating and updating source segmentation masks (for each j and l). This is done with the same tintfill algorithm that was applied in Sect. 3.2.2 to remove the source- and filament-like shapes. The resulting single-scale segmentation images of sources set all pixels belonging to a source to its number.

The scale j_F on which a source n becomes the brightest is referred to as the footprinting scale. It provides an initial estimate for its half-maximum size H_n = S_​jF (cf. Appendix B), which defines the initial footprint, that is, the entire area of all pixels making non-negligible contributions to the total flux of the source. From a practical point of view, getsf defines the initial footprint diameter of a circular source as ϕ_nH_n, where the footprint factor ϕ_n = 3. For the Gaussian sources (e.g., Fig. B.1), these footprints lead to the total fluxes that are underestimated by only 1.6%, well within the usual measurement uncertainties. Having detected the sources, getsf creates their initial footprints with the diameters {A, B}_Fn = ϕ_nH_n. The footprints become elliptically shaped in the wavelength-dependent measurements, reflecting the elongation of sources that is obtained from intensity moments. During the measurement iterations (Sect. 3.4.6), getsf changes ϕ_n to expand or shrink the footprints for those sources that are bright enough and whose intensity distributions indicate that their initial footprints are not optimal.

3.4.5. Detection of filaments in the combined images

Filaments are detected in ℱ_DjC with a completely new approach. In the getold algorithm (Sect. 2.4.4 of Paper II), intensity profiles at each pixel of the component of filaments are measured in four directions, and the pixel is deemed to belong to the crest (marks a skeleton point) if it has the highest value for each of the profiles. In practice, this simple approach sometimes creates artifacts at the filament end points, where the skeletons sometimes appear forked like a snake tongue. An important limitation of the getold skeletons is that they trace crests of the images of filaments without any dependence on the spatial scales.

The Herschel observations of nearby star-forming molecular clouds demonstrated that filaments are very complex, multiscale structures (e.g., Men’shchikov et al. 2010), unlike the simple case of the relatively round sources, whose intensities rapidly decrease in all directions from their peaks. Resolved sources are produced by the emission of dense, compact objects and may be reasonably well characterized by a single value of their half-maximum size (or spatial scale). In contrast to the sources, detection of filaments is fundamentally a scale-dependent problem, and a single skeleton that may be appropriate for a certain spatial scale cannot fully describe the complexity of the observed multiscale, profoundly substructured filaments. Resolved filaments often appear to be composed of thinner filaments on smaller scales, down to the angular resolution, and their widths, profiles, and crest intensities are quite variable along their skeletons.

The strong dependence of the observed filaments on spatial scales is illustrated in Fig. 13, which shows the surface densities of the filaments in three well-studied star-forming regions: Taurus, Aquila, and IC 5146. The observed images of the regions were downloaded from the Herschel Gould Belt Survey (André et al. 2010) archive⁴, and the hires surface densities 𝒟_13″ of the regions were computed from Eq. (8). Figure 13 shows the images of the spatially decomposed filaments on three selected scales: small, intermediate, and large. The images demonstrate that the observed filaments are highly substructured in the regions, and their shapes as traced by the skeletons are very different on various spatial scales. The skeletons, obtained on the small scales, are completely incompatible with the shapes and crests of the filaments on larger scales. The detected small-scale skeletons are often very curved, meandering back and forth even at the right angles, which implies a high degree of self-blending and leads to significant inaccuracies in the measured profiles and other derived properties of filaments. Therefore it is necessary to detect their skeletons on the scales that correspond to the widths of the structures being studied. Moreover, the small-scale substructures of the larger-scale filaments may even be the key to understanding the filament properties, the physical processes taking place inside them, and the formation of stars.

Instead of tracing the original image intensity profiles, getsf employs the Hilditch algorithm (Hilditch 1969), which skeletonizes two-dimensional shapes by erasing their outer pixels until the thinnest representation of the shapes is found. The original Hilditch algorithm has a deficiency in that the shapes oriented along the two main diagonals become completely erased during the iterations. To enable its application in getsf, the algorithm has been improved to preserve the diagonal skeletons.

Through the multiscale decomposition, getsf allows finding crests without any explicit analysis of the filament intensities. The single-scale images ℱ_DjC not only enhance the structures of the widths W ≈ S_​j, but also cause these filtered intensity distributions to become well centered on their zero-level footprints. The crests of the isolated decomposed filaments always approximate the medial axes of their footprints (cf. Figs. 12 and 13). If a single-scale filament blends with other filaments (or with itself), there is always a smaller scale on which it is isolated. This allows determining the scale-dependent skeletons as the medial axes of the zero-level filament masks.

The single-scale skeletons 𝒦_j are created using the Hilditch algorithm, with a width of three pixels to tolerate one-pixel displacements in the skeleton coordinates between scales. They are further accumulated over a limited range of scales to produce a set of N_K skeletons tracing the filamentary structures of various widths,

$$\begin{array}{c}\hfill {\mathcal{K}}_k={\displaystyle \sum _{j=J^{-}}^{J^{+}}}{\mathcal{K}}_j,k=1,2,\cdots ,N_{\mathrm{K}},\end{array}$$(34)

where J⁻ and J⁺ are the numbers of the smallest and the largest scales, S_​J⁻ = 2^−1/2S_​k and S_​J⁺ = 2^+1/2S_​k, in the accumulated skeleton 𝒦_k. The scale-dependent skeletons 𝒦_k sample the following scales:

$$\begin{array}{c}\hfill S_k=2^{1/2}{S}_{k-1},k=2,3,\cdots ,N_{\mathrm{K}},\end{array}$$(35)

where the scale S₁ = Ō is defined by Eq. (32) as the average angular resolution over the wavebands combined in ℱ_DjC (Sect. 3.4.3), and S_​NK = 4max_λ(Y_λ) is the largest spatial scale for the filament detection.

Each pixel of the accumulated skeleton 𝒦_k in Eq. (34) contains information on the filament detection significance ξ, defined as the number of scales between J⁻ and J⁺, on which the single-scale skeleton 𝒦_j contributes to 𝒦_k in that pixel. Depending on the filament intensity at the skeleton pixel, the significance range is 1 ≤ ξ ≲ ln2 (lnf)⁻¹ (≈14, assuming f ≈ 1.05, Appendix B). The algorithm automatically creates the final one-pixel-wide skeletons by thresholding: 𝒦_kξ = max(𝒦_k,ξ) with a default ξ = 2, and applying the Hilditch algorithm to the resulting shapes. Segmentation images of the skeletons 𝒦_kξ are computed using the tintfill algorithm, which sets all pixels belonging to a filament to its number.

3.4.6. Measurement of the sources

The source-measurement algorithm is an improved version of the one employed by getold (Sect. 2.6 of Paper I). Sources cannot be measured in their component 𝒮_λX (Sect. 3.2.3) because the subtracted background ℬ_λX contains substantial source residuals at low intensity levels (Fig. 8). The background ℬ_λX is derived specifically for the most complete and reliable source detection, not for accurate measurements. The sources are measured by getsf in the original ℐ_​λ after subtracting their backgrounds and deblending them from overlapping sources, which entails iterations. The background determination and deblending are more accurate for the sources with relatively small footprints. However, in crowded regions with larger areas of overlapping footprints and strongly fluctuating backgrounds, they may become very inaccurate.

The background ℬ_Fλ of each source is determined by a linear interpolation of ℐ_​λ across its footprint. The interpolation along two image axes and two diagonals is based on the adjacent pixels (not belonging to any source) outside the footprint, as was done by getold. To improve the background estimate in the presence of overlapping footprints, getsf evaluates the background only along those of the four directions for which the distances between the outside points being interpolated are within a factor of two of the smallest distance. For each pixel of the source footprint area, the background value is averaged between the directions used in the interpolation. The background ℬ_Fλ is median filtered using a sliding window with a radius O_λ, and the background-subtracted image of a source is then obtained as ℐ_​Sλ = ℐ_​λ − ℬ_Fλ.

In the measurements, the source coordinates x_n, y_n are known from the detection step and are kept unchanged. For the first measurement iteration, it uses the initial characteristic size H_n = S_​jF, provided by the detection algorithm (Sect. 3.4.4). The corresponding initial footprint {A, B}_Fn = ϕ_nH_n is a good approximation for only Gaussian sources, when H_n is close to the actual widths {A, B}_λn. However, the initial factor ϕ_n = 3 may strongly underestimate the footprints of the resolved power-law sources and overestimate those of the resolved flat-topped sources (see below). In the subsequent measurement iterations, the sizes and orientation {A, B, ω}_λn from the previous iteration are used.

The size derivation algorithm in getsf has become more accurate, hence it requires some clarifications. The half-maximum sizes were computed by getold using intensity moments (cf. Appendix F in Paper I) that give accurate sizes only for the sources with Gaussian shapes. In real-life observations, however, some sources are markedly non-Gaussian and their intensity moments give either over- or underestimated sizes, corresponding to the levels well below or above the half-maximum intensity.

The inaccuracies of the moment sizes become very large for the resolved sources with power-law intensity distributions. The simulated image of such a source shown in Fig. 14 has a half-maximum size of 10″. However, according to the intensity moments (over the entire image), the model source has a diameter of 76″. It is easy to see that this value corresponds to a level that is lower by an order of magnitude than the half-maximum intensity. The source size depends on the adopted footprint. Within the two footprints shown in the middle and right panels of Fig. 14, the moment sizes are 10.2 and 22.5″. The source flux is also underestimated by correspondingly large factors of 5.2 and 1.8.

Large inaccuracies of the half-maximum sizes also occur for the resolved starless cores that tend to have flat-topped shapes at short wavelengths (λ ≲ 250 μm, cf. Fig. 3), where the emission of their low-temperature interiors fades away. A simulated image of such a source is shown in Fig. 15, with the model half-maximum size of 49″. However, the intensity moments (over the entire image) indicate that its diameter is 31″, which corresponds to a level by a factor of 2 above the half-maximum intensity. In the simple example in Fig. 15, the source size and flux do not depend on the footprint size because the intensity profile in its outer parts is steep and the background is flat (zero).

The above examples demonstrate that the intensity moments do not provide accurate estimates of the half-maximum source sizes in the general case of arbitrary non-Gaussian intensity profiles. Therefore getsf determines accurate half-maximum sizes by the direct Gaussian interpolation of the source intensity distribution at its half-maximum and averaging the resulting distances from the source peak, thereby estimating an average radius h_λn. The source elongation E_Mλn = A_Mλn/B_Mλn and position angle ω_Mλn are computed independently from the intensity moments above the 10% level of the peak, excluding the low-intensity pixels that may be affected by the noise and background fluctuations. The major and minor half-maximum axes of the source are then computed from

$$\begin{array}{c}\hfill B_{\lambda n}={\mathrm{e}}^{-0.05(E_{\lambda n}-1)}{h}_{\lambda n}{E}_{\mathrm{M}\lambda n}^{-1/2},A_{\lambda n}=B_{\lambda n}{E}_{\mathrm{M}\lambda n},\end{array}$$(36)

where the (empirical) exponential factor converts the average radius h_λn into the equivalent-area radius (A_λnB_λn)^1/2 of an ellipse. The FWHM ellipse {A, B, ω}_λn from Eq. (36) is guaranteed to correspond to the source half-maximum intensity, in contrast to the ellipse estimated from the intensity moments. The moment sizes {A, B, ω}_Mλn are also computed by getsf because they contain independent information that can be useful for the analysis of the extracted sources.

During the measurement iterations (Sect. 2.6 of Paper I), getsf employs a footprint expansion and shrinkage algorithm to correct the footprint areas of those sources that need such adjustments. It is based on a simple observation that when a footprint area is too small, the source background contains a pedestal of the residual intensity distribution of the source (Fig. 14); when the source pedestal does not exist or is negative (Fig. 15), the footprint may be accurate or too large. The analysis is made in the regularized component 𝒮_λR from Eq. (25) without contribution from the complex background and filaments.

The presence of the background pedestal is indicated by the positive difference between the background values below the source and those in an external annulus just outside the source,

$$\begin{array}{c}\hfill 1.1B_{\mathrm{\Phi }\lambda n}>B_{\mathrm{\Psi }\lambda n}+D_{\mathrm{\Psi }\lambda n},\end{array}$$(37)

where B_Φλn is the median value within the footprint and B_Ψλn and D_Ψλn are the mean and the standard deviation inside the annulus. When the condition of Eq. (37) is fulfilled and the source is not too elongated (A_λn <  1.3B_λn) and bright enough (Ξ_λn >  50 and Ω_λn >  15, see Eq. (41)), getsf increases the factor ϕ_n by 5% before proceeding to the next measurement iteration. The footprint expansion terminates when the residual background pedestals (Fig. 14) are reduced as much as possible and the condition in Eq. (37) becomes false. As a final adjustment, the footprint is expanded once more by the factor 0.9 + 0.1(ϕ_n/3) to reduce the residual pedestal.

The footprints of the sources that do not need any expansion are attempted to be reduced in size. It is important to confine the footprints to the most local area occupied by the sources because oversized footprints may strongly decrease the accuracy of background subtraction and flux measurement for sources on complex (filamentary) backgrounds and in crowded areas. The need to shrink a source footprint is indicated by a negative difference between the background values below the source and in an external annulus just outside the source,

$$\begin{array}{c}\hfill 1.1B_{\mathrm{\Phi }\lambda n}<B_{\mathrm{\Psi }\lambda n}+D_{\mathrm{\Psi }\lambda n},\end{array}$$(38)

where the quantities are the same as in Eq. (37). When the condition of Eq. (38) is fulfilled, getsf decreases the factor ϕ_n by 2% before proceeding to the subsequent measurement iterations. The footprint shrinkage is completed when the condition in Eq. (38) becomes false, that is, when the reduced footprint causes a small residual background pedestal. In a final adjustment, the footprint is expanded by a factor of 1.1 to eliminate the pedestal created in the process (Fig. 15).

Extensive testing has shown that the simple footprint expansion and shrinkage algorithm performs well for most sources in complicated environments and backgrounds in both benchmarks and real-life observations. Despite the footprint expansion, total fluxes of the power-law sources may still remain underestimated by large factors because the faint outer areas of these sources, with an unknown full extent, vanish into the fluctuating backgrounds and noise.

After computing the background-subtracted images ℐ_Sλ, getsf deblends overlapping sources, calculating the peak intensities F_Pλn and the total fluxes F_Tλn for each source n. The iterative deblending algorithm employs the Gaussian shapes G_λn(x, y) defined by the source ellipse {A, B, ω}_λn and peak intensity F_Pλn. The intensity I_Sλ(x, y) is split between the source n and all overlapping sources n′ according to a fraction of the shape intensities,

$$\begin{array}{c}\hfill I_{\lambda n}(x,y)=\frac{G_{\lambda n}(x-x_n,y-y_n)}{{\displaystyle \sum _{n^{\prime }}}{G}_{\lambda n^{\prime }}(x-x_{n^{\prime }},y-y_{n^{\prime }})}{I}_{\mathrm{S}\lambda }(x,y),\end{array}$$(39)

where the summation is done over all surrounding sources whose footprints cover the pixel (x, y). The iterative deblending of the peak intensities starts with the original image values I_Sλ(x_n, y_n) of each source and proceeds with the splitting of the pixel values until I_λn(x_n, y_n) converges to the deblended peak intensity F_Pλn. After obtaining F_Pλn for all sources, getsf computes the deblended intensities I_λn(x, y) of all pixels within their footprints, estimates the ellipses {A, B, ω}_λn and {A, B, ω}_Mλn, and integrates the total fluxes F_Tλn. It also computes an independent flux estimate F_Gλn by integrating G_λn(x, y), which may only be accurate when a source shape resembles the two-dimensional Gaussian.

Uncertainties of the peak intensities F_Pλn are estimated by getsf as the standard deviations σ_Pλn, evaluated in the original image ℐ_​λ, in an elliptical annulus around each source n just outside its footprint. In heavily crowded fields, no local source-free annulus can be found near the sources, in which case the uncertainties are estimated from the more distant source-free pixels. The uncertainties σ_Tλn of the total fluxes F_Tλn are computed with the same assumptions as in getold (Sect. 2.6 of Paper I),

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{T}\lambda n}={\sigma }_{\mathrm{P}\lambda n}\frac{{(A_{\mathrm{F}\lambda n}{B}_{\mathrm{F}\lambda n})}^{1/2}}{{\varphi }_n{O}_{\lambda }},\end{array}$$(40)

where A_Fλn and B_Fλn are the major and minor axes of the source footprints.

It is convenient to define the detection significance Ξ_λn and the signal-to-noise (S/N) Ω_λn and Ψ_​λn, describing the detection and measurement properties of each extracted source,

$$\begin{array}{c}\hfill {\mathrm{\Xi }}_{\lambda n}=\frac{S_{\lambda \mathrm{D}{j}_{\mathrm{F}}n}}{{\sigma }_{\lambda \mathrm{S}{j}_{\mathrm{F}}}},{\mathrm{\Omega }}_{\lambda n}=\frac{F_{\mathrm{P}\lambda n}}{{\sigma }_{\mathrm{P}\lambda n}},{\mathrm{\Psi }}_{\lambda n}=\frac{F_{\mathrm{T}\lambda n}}{{\sigma }_{\mathrm{T}\lambda n}},\end{array}$$(41)

where j_F is the footprinting scale (Sect. 3.4.4) and S_​λDjFn is the intensity at the source position in 𝒮_λDjF (Sect. 3.4.1). The above quantities can be combined together to characterize the overall “goodness” of a source,

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{\lambda n}=\frac{{\mathrm{\Xi }}_{\lambda n}}{5}\frac{{\left({\mathrm{\Omega }}_{\lambda n}{\mathrm{\Psi }}_{\lambda n}\right)}^{1/2}}{2}\frac{B_{\lambda n}}{A_{\lambda n}},\end{array}$$(42)

normalized such that all acceptable sources in the extraction catalogs have Γ_λn ≳ 1. The sources with Γ_λn ≲ 1 may have quite unreliable measurements in waveband λ. The corresponding global quantities Ξ_n and Γ_n describe the source detection significance and goodness, respectively, in all wavebands,

$$\begin{array}{c}\hfill {\mathrm{\Xi }}_n={\left({\displaystyle \sum _{\lambda }}{\mathrm{\Xi }}_{\lambda n}^2\right)}^{1/2},{\mathrm{\Gamma }}_n={\left({\displaystyle \sum _{\lambda }}{\mathrm{\Gamma }}_{\lambda n}^2\right)}^{1/2}.\end{array}$$(43)

The getsf source extraction catalogs contain detailed headers, documenting the extraction parameters and explaining the tabulated quantities. Each data line presents the source number n, coordinates x_n, y_n (in pixels), world coordinates α_n, δ_n (computed with the xy2sky utility, Mink et al. 2002), global flag f_n, significance Ξ_n, and goodness Γ_n,

$$\begin{array}{c}\hfill n{x}_n{y}_n{\alpha }_n{\delta }_n{f}_n{\mathrm{\Xi }}_n{\mathrm{\Gamma }}_n,\end{array}$$(44)

followed (in the same line) by the measured quantities in each of the N_W wavebands,

$$\begin{array}{c}\hfill {\left(f_{\lambda n}{\mathrm{\Xi }}_{\lambda n}{\mathrm{\Gamma }}_{\lambda n}{F}_{\mathrm{P}\lambda n}{\sigma }_{\mathrm{P}\lambda n}{F}_{\mathrm{T}\lambda n}{\sigma }_{\mathrm{T}\lambda n}{A}_{\lambda n}{B}_{\lambda n}{A}_{\mathrm{M}\lambda n}{B}_{\mathrm{M}\lambda n}{\omega }_{\lambda n}\right)}_{N_{\mathrm{W}}},\end{array}$$(45)

where f_λn is a wavelength-dependent flag. In addition to this information, an expanded version of the catalog adds (to the same line) the Gaussian flux F_Gλn, characteristic size S_​jF, footprint factor ϕ_n, and footprint axes A_Fλn, B_Fλn. For surface density images, the F_Gλn column is replaced with source mass M_ƛn.

It is necessary to emphasize that sources from extraction catalogs must always be carefully selected (for each waveband separately) to ensure that only sufficiently good and accurately measurable sources are used in further analysis. This is especially important for the multiwavelength extraction catalogs, where sources can be prominent in one waveband and completely undetectable or not measurable in another. The getsf catalogs provide various quantities to enable the evaluation and selection of only acceptable sources and recommended the following selection criteria:

$$\begin{array}{l}{\Xi }_{\lambda n}>1\wedge {\Gamma }_{\lambda n}>1\wedge {\Omega }_{\lambda n}>2\wedge {\Psi }_{\text{}\lambda n}>2\hfill \\ \wedge A_{\lambda n}<2B_{\lambda n}\wedge A_{\text{F}\lambda n}>1.15A_{\lambda n}.\hfill \end{array}$$(46)

These empirical conditions, based on numerous test results obtained in various benchmarks (Pouteau et al., in prep.; Men’shchikov 2021), and verified in applications to a variety of observed images (e.g., Sect. 4), ensure that the selected subset of sources is reliable (not contaminated by significant numbers of spurious sources) and that selected sources have acceptably accurate measurements.

3.4.7. Measurements of the filaments

Filaments are measured in their background-subtracted ℱ_λY, derived in Sect. 3.2.3. When the maximum size Y_λ of the filaments of interest is estimated sufficiently accurately (Sect. 3.1.3), their background ℬ_λY does not reveal any filamentary residuals. Nevertheless, the background may well have substantial inaccuracies, especially when the filaments are wide and blended (Fig. 8). Observed filaments are the two-dimensional projections of the complex three-dimensional structures, which are much more difficult to disentangle, deblend, measure, and analyze than sources with their well-defined round shapes and compact footprints.

Sources can be represented by their peak intensity and half-maximum size, but filaments are extremely complicated in their shapes and widths, often interconnected with each other and with various nearby branches, and have variable intensity along their crests. It is quite clear that blending of the structures is a major source of large inaccuracies in the measured quantities of general interest (widths, fluxes, masses, profiles) and in other properties, derived from the measurements.

Another difficulty in understanding filaments (distinct physical structures) is that the filament length cannot be determined objectively. In most cases, it is quite unclear where a physical filament starts, where it ends, and which branches of the three-dimensional filamentary network belong to that filament. Fortunately, the global properties of the entire filaments (even if the latter could be clearly defined) are not as important for studying star formation as the local properties of their relatively short segments that develop appropriate physical conditions for the formation of prestellar cores.

The approach that is adopted in getsf is to simplify the very complex problem by separating all branches of the skeleton network, converting the latter into the simple, non-branching skeletons. The set of non-branching skeletons is derived during the segmentation of the skeletons 𝒦_kξ, the last step of the filament detection process (Sect. 3.4.5). The simplified skeletons enable an easy selection and better measurements of only the well-behaving, preferably isolated (not blended), and relatively straight parts of the filaments. No attempt is made by getsf to deblend filaments because a general algorithm for accurately deblending them is not available.

The segmentation image of all skeletons is scanned to trace each skeleton n and find coordinates of all its pixels; to smooth the skeletons, the integer coordinates of their pixels are averaged within a seven-pixel window. The resulting high-resolution coordinates x_n(i),y_n(i) of each skeleton point i are cataloged, together with the local position angles ϑ_n(i) of the skeleton direction and α_n(i),β_n(i) of the left and right normals. A normal is called left (α) or right (β) depending on which side it is located from the first skeleton point to the last. With an adopted distance to the observed region, getsf converts the angular units of the pixels into parsecs and measures each filament as a function of the length l along its skeleton and the distance r along its normals. If the distance is unknown or unspecified, a default distance of 100 pc is used; the measurements can be scaled to another distance after the extraction.

The observed filaments usually meander, hence their skeleton normals diverge from each other on one side and intersect with each other on the other side. In the absence of deblending, more accurate measurements for them are usually obtained from the one-sided quantities that correspond to the side on which the filament is the least affected by blending with itself and with other nearby structures. The filament surface density (or intensity) profiles and their full half-maximum widths are cataloged as the one-sided quantities D_{α|β}n(l, r) and W_{α|β}n(l) and as the average quantities D_n(l, r) and W_n(l) between the two sides. Also cataloged are the corresponding average profiles D_{α|β}n(r) and D_n(r) along the skeleton with their standard deviations ς_{α|β}n(r) and ς_n(r), as well as the median widths W_n and the slopes γ(r) of the filament profiles.

Although the total length L_n of a skeleton and mass M_n of a filament may not always be objective and physically meaningful quantities (see the discussion above), getsf derives the mass by direct integration of ℱ_λY within a filament footprint, assuming that the image is obtained from surface densities,

$$\begin{array}{c}\hfill M_{\{\alpha |\beta \}n}=2\mathrm{\mu }{m}_{\mathrm{H}}\int {\displaystyle \underset{{\mathrm{{\rm Y}}}_{\{\alpha |\beta \}n}}{\int }}{\mathcal{F}}_{\lambda Yn}(x,y)\mathrm{d}x\mathrm{d}y,\end{array}$$(47)

where M_{α|β}n are the one-sided mass estimates, from which the average mass M_n between the two sides is obtained. The one-sided footprints Υ_{α|β}n used in Eq. (47) are defined as the areas between the skeleton and the maximum extent of the filament on either side. In practice, a filament footprint Υ_n is the set of all pixels whose shortest distances from the skeleton are smaller than the filament normals.

When the filament mass M_{α|β}n and length L_n are known, the one-sided estimates of the average linear density⁵ of the entire filament are readily obtained,

$$\begin{array}{c}\hfill {\overline{\mathrm{\Lambda }}}_{\{\alpha |\beta \}n}=M_{\{\alpha |\beta \}n}{L}_n^{-1},\end{array}$$(48)

together with the average linear density ${\overline{\mathrm{\Lambda }}}_n$ between the two sides. The linear density of filaments is also computed by getsf as a function of the coordinate l along their skeletons,

$$\begin{array}{c}\hfill {\mathrm{\Lambda }}_{\{\alpha |\beta \}n}(l)=2\mathrm{\mu }{m}_{\mathrm{H}}{\displaystyle \underset{0}{\overset{R_{\{\alpha |\beta \}n}(l)}{\int }}}{\mathcal{F}}_{\lambda Yn}(l,r)\mathrm{d}r,\end{array}$$(49)

where the integration limits R_{α|β}n(l) along the left and right normals are chosen at zero surface density values or at a radial distance of the profile minimum at which the filament becomes blended with another structure. The median one-sided linear densities Λ_{α|β}n for the entire length L_n of the filament and its average linear density Λ_n are also computed and cataloged. The linear density values from Eqs. (48) and (49) are expected to be similar to each other for the well-behaved filaments.

4. Applications to observed regions

The multiscale, multiwavelength source- and filament-extraction method presented in Sect. 3 was very extensively tested using ∼40 images that were observed with different instruments and both ground-based and orbital telescopes during the past two decades. Multiwaveband observations of star-forming regions obtained in the Herschel Gould Belt Survey (André et al. 2010) and HOBYS (Motte et al. 2010) key projects, as well as the most recent interferometric images observed in the ALMA-IMF program (Motte et al., in prep.), played an important role in validating getsf.

The new extraction method has demonstrated very good results in ALMA benchmarks (Pouteau et al., in prep.) on images, created from a magnetohydrodynamic (MHD) simulation of a star-forming region (Ntormousi & Hennebelle 2019) that was populated with model cores and processed by the casa task simobs (McMullin et al. 2007) to resemble the real ALMA observations (Louvet et al., in prep.). Furthermore, getsf has been applied to source extraction in 15 regions of the ALMA-IMF program and 12 infrared dark clouds of the ASHES survey (Sanhueza et al. 2019; Li et al. 2020). However, the most significant and definitive testing and validation of extraction tools is achieved with simulated benchmarks for which everything is fully known about their components. The second paper (Men’shchikov 2021) presents a quantitative analysis of getsf extractions using several variants of the new benchmark (Sect. 2) and old benchmark (Papers I and III).

Sections 4.1–4.8 illustrate the performance of getsf on nine images obtained with different telescopes: XMM-Newton, the Galaxy Evolution Explorer (GALEX), Hubble, Spitzer, Herschel, the Atacama Pathfinder Experiment (APEX), the James Clerk Maxwell Telescope (JCMT), and ALMA from the X-ray domain to the millimeter wavelengths. These examples are presented to demonstrate that the method is applicable to a wide variety of observed images, visualizing the effects of the separation of structural components and flattening of detection images. Scientific analyses and discussions of these results, as well as their comparisons with previous studies, are beyond the scope of this paper. This can be accomplished using the corresponding extraction catalogs that are available on the getsf website⁶.

4.1. Supernova remnant RXJ 1713.7−3946

RXJ 1713.7−3946 was observed with XMM-Newton (EPIC camera) in the X-ray waveband (0.6−6 keV) centered at 0.0024 μm. The 0.5° ×0.5° image in Fig. 16 is a mosaic of multiple observations⁷, first presented in Acero et al. (2017). With an average angular resolution of 7″, it reveals the southeast segment of the supernova remnant shell that may have been created by the explosion of the historical supernova SN 393, whose center of explosion is located beyond the upper right image corner. For this source and filament extraction with getsf, maximum sizes {X|Y}_λ = {15, 25}″ were adopted (Sect. 3.1.3).

The observed X-ray image (Fig. 16) has relatively low counts of the detected photons per pixel and high levels of Poisson noise. The image is contaminated by linear artifacts and several spurious single-pixel spikes. The latter may appear in these images when just one or several photons are detected at an edge of the mapped area.

The image features several elongated shock fronts created by the expanding supernova shell, and a number of faint and bright point sources, all of them well isolated. The getsf extraction greatly simplified the image by separating the components of sources 𝒮_λ, filaments ℱ_λ, and their backgrounds ℬ_λ{X|Y}. The small-scale fluctuation levels across the observed image are only within a factor of two, therefore the improvement caused by the flattening is not clearly discernible in 𝒮_λ. However, the images of standard deviations show that the flat source detection image 𝒮_λD has uniform fluctuations over the entire image, which is beneficial for source detection.

The extraction catalog contains measurements of 41 sources, all of them selected as acceptably good by Eq. (46). Although the spurious one-pixel spikes were not removed before the extraction, getsf identified them as such (red squares in Fig. 16) and eliminated them from the catalog. Despite the faintness of the observed X-ray image and the Poisson noise, the three prominent shocks of the supernova shell become clearly visible and are extracted in the filament component.

4.2. Star-forming galaxy NGC 6744

NGC 6744 was observed with GALEX in a far-ultraviolet (FUV) waveband (1350−1750 Å) centered at 0.15 μm (Lee et al. 2011). The 0.4° × 0.4° image⁸ in Fig. 17 with an angular resolution of 4″ shows the spiral galaxy, which is considered to be similar to our own Galaxy. Despite noisiness of the FUV image, it displays the spiral arms with many hundreds of unresolved emission sources. These are the regions of ongoing star formation, heated by the embedded young massive stars. For this source and filament extraction, maximum sizes {X|Y}_λ = 20″ were adopted (Sect. 3.1.3).

Separation of the structural components by getsf provided independent images of sources 𝒮_λ, filaments ℱ_λ, and their backgrounds ℬ_λ{X|Y} (Fig. 17). Fluctuation levels in the observed image vary within a factor of two, largely in the central, brighter part of the galaxy. Flattening of the 𝒮_λ component effectively equalized the fluctuations across the detection image 𝒮_λD, improving the extraction results.

The source catalog contains measurements of 1169 sources, 1130 of which are selected as acceptably good by Eq. (46). Most of the sources likely correspond to the star-forming regions along the galactic spiral arms; many of them overlap with each other, hence they required deblending for accurate measurements of their fluxes. The filaments extracted in the ℱ_λ component represent the spiral arms and their branches. The 147 skeletons trace the simple, non-branching segments of the filamentary network (Sect. 3.4.5).

4.3. Supernova remnant NGC 6960

NGC 6960 was observed with Hubble in five UVIS wavebands (0.5−0.8 nm) centered at 0.6 μm, within the frame of the Hubble Heritage project (Mack et al. 2015, PI: Z. Levay). The small 73″ × 73″ image⁹ in Fig. 18 with an angular resolution of 0.2″ represents a small fragment of the Veil Nebula, which is a segment of the Cygnus Loop, the large expanding shell of a supernova remnant (Fesen et al. 2018). For this source and filament extraction with getsf, maximum sizes {X|Y}_λ = {0.5, 2}″ were adopted (Sect. 3.1.3).

The observed image (Fig. 18) is dominated by impressive fine filamentary structure of the nebula, seen in emission of a number of atomic lines. Many unresolved intensity peaks of sources are less prominent on this bright backdrop. The structural components were separated by getsf in the independent images of sources 𝒮_λ, filaments ℱ_λ, and backgrounds ℬ_λ{X|Y}; together with the flattening of detection images, this greatly facilitates their extraction and analysis.

The source catalog contains measurements of 786 sources, 690 of which are selected as acceptably good by Eq. (46). The strings of sources that run up through the middle of the image are the spurious peaks created by the linear artifacts. The spurious spikes were not cut out of the image before this extraction to illustrate that they need to be removed to avoid contamination of the source catalogs. The finely structured filamentary network of the nebula that is extracted by getsf in the ℱ_λ component comprises 100 skeletons representing its simple, non-branching segments (Sect. 3.4.5).

4.4. Star-forming cloud L 1688

L 1688 was observed with Spitzer in the IRAC 8 μm waveband (Evans et al. 2009). The 1° × 1° image¹⁰ in Fig. 19 with an angular resolution of 6″ shows a complex intensity distribution in this well-known star-forming region, with the background varying by almost two orders of magnitude and many sources situated in both faint and bright background areas. For this source and filament extraction with getsf, maximum sizes {X|Y}_λ = 30″ were adopted (Sect. 3.1.3).

The clean separation of the components of sources 𝒮_λ and filaments ℱ_λ from their backgrounds ℬ_λ{X|Y} provided by getsf (Fig. 19) represents an obvious improvement over the results obtained with getimages (Fig. 6 in Paper III). The old method of background derivation was indiscriminate with respect to the shapes of the components, hence the background-subtracted image also contained some filamentary structures on small scales. In contrast to getimages, which produced a single background, getsf derived and subtracted individual backgrounds for 𝒮_λ and ℱ_λ. The component 𝒮_λ of sources (Fig. 19) is completely free of the elongated structures. The standard deviations 𝒰_λ reveal that the small-scale background fluctuation levels vary by roughly three orders of magnitude across the image. If not equalized, the fluctuations would be extracted as numerous spurious sources and contaminate the source catalog. The very effective flattening of the detection image 𝒮_λD leads to a much more reliable extraction.

Several bright unresolved sources in the lower part of the observed image have very wide power-law wings and cross-like artifacts that are induced by the complex PSF of the Spitzer IRAC camera at 8 μm. Their intensity profiles are markedly non-Gaussian, and for a proper measurement of their integrated fluxes, getsf expanded their footprints by factors ∼4–15 using the footprint expansion algorithm (Fig. 14). The cross-like artifacts from the PSF were interpreted by getsf as filaments and were moved to the filament component, thereby improving 𝒮_λD for source detection. In addition to the cross shape, the complex PSF has ∼20 faint peaks that surround the main beam. They were extracted as several spurious sources, surrounding the brightest peaks; they must be eliminated in a post-extraction analysis.

The source catalog gives measurements of 1474 sources, 1162 of which are selected as acceptably good using Eq. (46). The filament component produced by getimages (Fig. 7 in Paper III) contains only the brightest parts of the filaments, their fainter intensities are missing. In contrast, getsf determines the intensity distributions down to very low intensity levels, with 286 skeletons tracing the simple, non-branching segments of the filaments (Sect. 3.4.5).

4.5. Embedded starless core L 1689B

L 1689B, one of the nearest well-resolved starless cores (a distance of 140 pc) embedded in a resolved filament, was observed with Herschel in five PACS and SPIRE wavebands (Ladjelate et al. 2020). The 160−500 μm images¹¹ and Eq. (8) were used to compute a 1.1° × 1.1° surface density image 𝒟_13″ in Fig. 20 with a resolution of 13.5″ to illustrate the new extraction method on a single image. In addition to the reduction of the number of images, the use of surface densities allows getsf to catalog physical parameters of the core and filament. For this extraction, maximum sizes {X|Y}_ƛ = {90,180}″ were adopted (Sect. 3.1.3).

The 𝒟_13″ image in Fig. 20 presents L 1689B in the wide filamentary structure near the edge of a diffuse cloud, all components are blended. The filament surface density is a factor of ∼5 below the peak surface density N_H2 = 3.7 × 10²² cm⁻², whereas at the values, lower by just a factor of 2, a round shape of the source becomes distorted by its complex environment. Separation of the components by getsf greatly simplifies the image, isolating the structures in their individual images 𝒮_ƛ, ℱ_ƛ, and ℬ_ƛ{X|Y} (Fig. 20). Subsequent flattening of the small-scale fluctuation levels allowed a reliable identification of the filament and sources in both low- and high-background areas of the observed image. The extraction catalog contains measurements of 20 sources, 12 of which are selected as acceptably good by Eq. (46). The single skeleton was obtained on spatial scales of ∼200″, corresponding to the maximum width Y_ƛ.

The main physical parameters of the starless core L 1689B, M = 0.6 M_⊙ and N_H2 = 10²² cm⁻², are underestimated because of the inaccuracies (Appendix A) of the standard surface density derivation approach (Sect. 3.1.2). The errors and correction factors can be found using the benchmark models from Sect. 2.2. A model of the critical Bonnor–Ebert sphere with T_BE = 14 K, M = 1 M_⊙, and N_H2 = 2.5 × 10²² cm⁻² has an FWHM size of 57″, almost the same as the size A = 58″ of L 1689B, measured by getsf. However, in the derived 𝒟_13″ image, the same model has M = 0.66 M_⊙ and N_H2 = 10²² cm⁻², implying correction factors of 1.5 and 2.5 for the mass and peak surface density, correspondingly. After correction, the measured mass of L 1689B becomes M ≈ 0.9 M_⊙; masses of the other sources in the image are lower by (at least) a factor of ∼10 . The filament measurements (Sect. 3.4.7) give its median value N₀ = 3.8 × 10²¹ cm⁻², length L = 0.8 pc, half-maximum width W = 0.14 pc (205″), mass M = 15 M_⊙, and linear density Λ = 14 M_⊙ pc⁻¹; the values are little affected by the fitting inaccuracies.

4.6. Star-forming cloud NGC 6334

NGC 6334 was observed with APEX at 350 μm, equipped with the ArTéMiS camera (André et al. 2016). The 0.5° ×0.5° image¹² in Fig. 21 with an angular resolution of 8″ represents an improvement by a factor of 3 with respect to the Herschel images at 350 μm. Subtraction of the correlated sky noise resulted in an image without signals on spatial scales above 120″ (André et al. 2016). Therefore the large-scale background and the zero level of the image are not known, and the structures in the image are smaller than the largest scale. Fortunately, these observational problems are entirely unimportant for getsf. For this source and filament extraction, maximum sizes {X|Y}_λ = 30″ were adopted (Sect. 3.1.3).

The observed image of NGC 6334 displays complex blended structures of various shapes and intensities (Fig. 21), including substantial numbers of negative areas and artifacts from the data reduction and map-making algorithms. The separated 𝒮_λ component shows all source-like peaks very clearly, even those that are hardly visible in the original image, because getsf is able to distinguish sources from the elongated filamentary shapes. Many of the sources overlap each other, therefore they require deblending for accurate measurements. The background ℬ_λY of filaments is fairly low, hence its subtraction enhanced the visibility of filaments in ℱ_λ only little. Nonuniform small-scale fluctuations in 𝒮_λ were effectively equalized in the detection image 𝒮_λD by the flattening algorithm.

The source catalog contains measurements of 124 sources, 91 of which are selected as acceptably good by Eq. (46). In the component of filaments, getsf identified 26 skeletons, tracing the simple, non-branching segments of the filaments (Sect. 3.4.5) on spatial scales of ∼30″, corresponding to the maximum width Y_λ.

4.7. Star-forming cloud Orion A

Orion A was observed with JCMT at 450 and 850 μm with the SCUBA-2 camera (Lane et al. 2016) with angular resolutions of 9.8 and 14.6″, respectively. The 0.86° ×0.86° image¹³ at 850 μm in Fig. 22 displays the northern part of the integral-shaped filament (ISF). Like with other ground-based submillimeter observations that must subtract sky background, large-scale emission in the images has been filtered out (Kirk et al. 2018). A visual estimate suggests that the image contains substantial signal on spatial scales of up to ∼100″. For the two-wavelength getsf extraction, employing both 450 and 850 μm images, maximum sizes {X|Y}_{450|850}  =  {20, 30, 30, 45}″ were adopted (Sect. 3.1.3).

The 850 μm image of the ISF in Fig. 22 reveals the small-scale structures of the area most clearly because of the spatial filtering effect of the observational technique. However, the central bright part of the ISF remains blended, and the spatial decomposition by getsf helps isolate the sources 𝒮_λ in that area from the filaments ℱ_λ and their backgrounds ℬ_λ{X|Y}. The background of filaments is found to be slightly negative, except in its central bright round area. In comparison with an average value of small-scale fluctuations in 𝒮_λ, they are larger by a factor of 2.7 in the central zone and lower by a factor of 1.7 in the lower right corner. The standard deviations 𝒰_λ reveal imprints of the five overlapping scans from the observations. The flattening algorithm of getsf effectively equalizes them and creates the flat detection images {𝒮|ℱ}_λD of sources and filaments, improving their detection reliability.

The two-band source extraction in ISF with getsf cataloged 344 sources, detected and measured in both wavebands simultaneously. Only 257 and 212 sources at 450 and 850 μm, respectively, are selected as acceptably good by Eq. (46); the S/N for the remaining detections is too low or they have other defects that are identified by the measurements. Two additional getsf extractions, done on each image independently, resulted in catalogs with 319 and 283 sources at 450 and 850 μm, respectively. Independent extractions ignore the valuable information from the other image, hence there are higher chances of spurious sources. With the additional condition that cataloged sources must be detected in both images, the combined extraction catalog contains 223 sources; 196 and 183 of these sources at 450 and 850 μm, respectively, are acceptably good. They represent the most reliable sources in the images, hence it is highly unlikely that there are spurious sources among them.

The missing large-scale emission of the SCUBA-2 image helped getsf expose the many relatively faint, narrow filaments within the wide, massive ISF. In the ℱ_λ component, getsf identified 267 and 199 simple, non-branching segments of the filaments (Sect. 3.4.5) at 450 and 850 μm, respectively, on transverse scales of 28 and 39″. This is similar to the existence of narrow sub-filaments on small scales within the resolved Taurus, Aquila, and IC 5146 filaments (Fig. 13) and consistent with the recent ALMA observations of ISF (Hacar et al. 2018).

4.8. Star-forming cloud W 43-MM1

W 43-MM1 was observed with the 12 m array of the ALMA interferometer (baselines 13−1045 m) in the 233 GHz band centered at 1300 μm (Motte et al. 2018; Nony et al. 2020). The small 68″ × 68″ image in Fig. 23 with an angular resolution of 0.44″ contains spatial scales of up to 12″, beyond which the interferometer was insensitive to the emission. For this source and filament extraction with getsf, the maximum size {X|Y}_λ = {0.8, 1.3}″ was adopted (Sect. 3.1.3).

The interferometric image of W 43-MM1 (Fig. 23) shows a cluster of relatively bright sources, some of them blended, and three faint filamentary structures that appear to connect them. Separation of the components of sources 𝒮_λ and filaments ℱ_λ confirms that most sources are concentrated on (or near) the faint continuous filaments. Almost the entire background ℬ_λY of the filamentary structures is negative, which is caused by the missing large scales in the observed images.

The small-scale fluctuation levels steeply increase toward the image center by more than an order of magnitude (Fig. 23), as evidenced by the standard deviations 𝒰_λ. The small-scale structured noise from the interferometric observations have both round or irregular, elongated shapes. Consequently, the separation of structural components by getsf leads to many faint peaks in 𝒮_λ and ℱ_λ. The flattening algorithm equalizes the fluctuation levels very effectively, providing reliable detection of sources in the flat {𝒮|ℱ}_λD. If not suppressed, such highly variable structured noise would produce many spurious sources and filaments in the central area of the image.

This ALMA image of W43-MM1 contains only moderate numbers of sources and filaments. The extraction catalog contains measurements of 44 sources, and all of them are selected as acceptably good by Eq. (46). This simple field allows a visual verification that they all are the true sources and are not contaminated by the noise fluctuations. In the filament component, getsf identified 15 skeletons, tracing the simple, non-branching segments of the filaments (Sect. 3.4.5) on spatial scales of ∼2″, similar to the maximum width Y_λ adopted for the extraction.

5. Strengths and limitations

5.1. Strengths

In contrast to the other methods, getsf extracts sources and filaments simultaneously by combining available information from all wavebands. Its flexible multiwavelength design enables handling of up to 99 images, not necessarily all of them observed in different wavebands. The maximum number of images is arbitrary, representing the largest two-digit integer number used in the output file names; the code can be updated to use a higher value if required for some applications. Any subset of the input images that is deemed beneficial for the detection purposes can be used to detect the sources and filaments, whereas measurements of the identified structures are provided for all input images. In a nonstandard application, the method can also be employed with the position-velocity cubes if they are split into separate images along the velocity axis (getold was used in this way by Shimajiri et al. 2019).

The images that are selected for detection are spatially decomposed to isolate the contributions of similar scales (Appendix B) and are then combined in a wavelength-independent set of single-scale detection images (Sect. 3.4.3). This eliminates the necessity of associating independent detections across wavelengths in images with greatly different angular resolutions and improves the detection and measurement accuracy. For example, positional association of nearby sources detected at 160 μm and completely blended into a single clump at 500 μm does not make sense.

Separation of structural components in the images of highly structured observed regions in space provides independent images of sources, filaments, and their backgrounds (Sect. 3.2), which is highly beneficial for the analysis and interpretation of observations. Flattening of detection images equalizes the (nonuniform) small-scale background and noise fluctuations (Sect. 3.3). This greatly simplifies the images and allows reliable detection of sources and filaments in decomposed single-scale images using a constant threshold, with a very low rate of spurious sources.

Sources and filaments of any size and width can be extracted by getsf provided that they are significantly smaller than the image. Only the maximum size of the structures of interest must be specified for each image in order to limit the range of spatial scales considered and the size of the structures to be measured and cataloged. The single parameter of the observed images that getsf needs to know is the maximum size, which is determined from the images by users (Sect. 3.1.3) on the basis of their research interests. This single constrained parameter reduces the dependence of the extraction results on the human factor to a minimum and makes their analysis and derived conclusions as objective as possible.

The numerical code is designed to be user-friendly and easy to run, providing diagnostics to help users avoid common problems. It verifies the getsf configuration, input images, and masks for consistency, and it suggests solutions in various circumstances during extractions. The software includes 21 utilities and scripts (Appendix C), providing all kinds of image processing necessary for getsf to run and more. They include the fitfluxes utility for spectral energy distribution fitting of source fluxes or image pixels (and mass derivation) and the hires script that computes the high-resolution surface density images (Sect. 3.1.2). Most of the utilities are very useful for command-line image manipulations, even without source and filament extractions.

5.2. Limitations

The method is designed and expected to work for the images that are not very sparse: most pixels must contain detectable signals (measurable data). Examples of the images for which getsf might not produce reliable results are some extremely faint X-ray or UV low-count images with isolated spiking pixels that are surrounded by large areas of pixels that were not assigned any detectable signal. For such nonstandard images, getsf would still work and complete extractions, but its results might not be reliable because the method relies on the standard deviations of the background or noise fluctuations outside structures, whose values may not correctly represent the statistics of the observed data in these images. On the other hand, the images for getsf extractions must not be extremely smooth: they must have some variations on scales of about the angular resolution. However, such smooth images can easily be made perfectly suitable for getsf just by adding Gaussian noise at some faint level that does not alter the structures of interest.

Separation of sources from filaments is not (and cannot be) perfect. It leaves very faint residuals of sources that end up in the filament component. In practice, this is not important because most of the residuals are too faint (Fig. 12) to affect the filament properties. The background of very wide and/or overlapping filaments is likely to be derived less accurately than that of the narrower and/or isolated filaments because the filaments are separated from the wider background areas. Filaments that are separated from wider background peaks of comparable widths are likely to receive some contribution from the background (Fig. 8). In very rare cases, the footprint of a bright power-law peak might not be sufficiently expanded, which leads to an underestimated flux.

The method takes quite considerable time to complete extractions, although getsf was optimized to run as fast as possible. The aim of its design was to produce extraction results that are as reliable as possible because completeness and accuracy, not speed, are of prime importance in astrophysical research. The runtime for the getsf applications presented in Sect. 4 is in the range of three hours to a week (the images with 430² to 2000² pixels and file sizes of 800 KB to 16 MB). The two-wavelength extraction of sources and filaments for the subfield of Orion A described in Sect. 4.7 took 43 h and required ∼10 GB of disk space. The total processing time with getsf depends on the numbers of pixels, wavelengths, iterations, detected sources and filaments, and on the processor and file system speed and load. A source extraction run on eight large images, each with 4800² pixels (92 MB file size), that detects and measures ∼3000 sources, may need about three weeks and ∼200 GB of disk space. Most of the time getsf spends in the iterative separation of structural components: the actual extraction of sources and filaments takes less than 10% of the runtime. For the source extraction alone, the execution time is halved. In a properly planned research, the processing time is almost never a limiting factor: much more time is usually spent on the analysis and interpretation of the information delivered by the extraction and on describing the findings in a paper.

Many intermediate images are produced in the getsf extractions at each wavelength (for spatial decomposition, iterations, etc.), hence they require large storage space. Between hundreds of MB and GB may be necessary for an extraction, depending on the image size and the numbers of wavebands and iterations. It is necessary to keep many images until the end of the extraction process; however, most of them may be deleted by getsf after the extraction has finished. The extraction results themselves represent only ∼20% of the total size of the extraction directory. Computers with sufficiently large random access memory are required to run getsf extractions on very large images. For the above range of image sizes, between 8 and 64 GB may be necessary (the more memory, the better). The actual memory usage strongly depends on the number of sources being detected and measured. Numbers of sources up to ∼15 000 do not pose any problems to getsf, but substantially larger numbers of detected sources require very large memory and long execution time.

6. Conclusions

This paper presented getsf, the new multiscale method for extracting both sources and filaments in astronomical images using separation of their structural components. It is specifically designed to handle multiwavelength sets of images and extremely complex filamentary backgrounds, but it is perfectly applicable to a single image or very simple backgrounds. The new code is freely downloadable from its website¹⁴, from the Astrophysics Source Code Library¹⁵, and also available from the author.

The main processing steps of getsf include (1) preparation of a complete set of images and derivation of high-resolution surface densities, (2) spatial decomposition of the original images and separation of the structural components of sources and filaments from each other and from their backgrounds, (3) flattening of the residual noise and background fluctuations in the separate images of sources and filaments, (4) spatial decomposition of the flattened components of sources and filaments and their combination of the over wavelengths, (5) detection of sources (positions) and filaments (skeletons) in the combined images of the components, and (6) measurements of the properties of the detected sources and filaments and creation of the output catalogs and images. Like its predecessor (getold, Papers I–III), getsf has a single user-definable parameter (per wavelength), the maximum size of the structures of interest to extract. All internal parameters of getsf have been calibrated and verified in numerous tests using various images from simulations and observations to ensure that the method works well in all cases.

This paper formulated hires, the algorithm for the derivation of high-resolution surface densities and temperatures from the diffraction-limited multiwavelength far-infrared and submillimeter continuum observations, such as those obtained with Herschel. A substantial improvement over the original algorithm (Palmeirim et al. 2013) is the angular resolution of the derived surface densities that may become as high as that of the shortest-wavelength image of a sufficient quality. In the case of the Herschel observations, the resolution may be as high as 5.6″ for the slow scanning speed (20″ s⁻¹) or 8.4″ for the fast parallel mode (60″ s⁻¹). If the 70 μm image appears too noisy, excessively contaminated by the emission of polycyclic aromatic hydrocarbons or transiently heated very small dust grains, or if it cannot be used for other reasons, then the highest resolution of surface densities is limited to that of the 100 or 160 μm images, that is, to 6.8−11.3″ or 8.4−13.5″, for the slow- or fast-scanning modes, respectively. These high-resolution surface density images are especially useful for the detailed studies of the highly complex structural diversity in the observed images and for deeper understanding of the physical processes within the heavily substructured filaments and their relation to the formation of stars.

This paper described the set of simulated multiwavelength benchmark images for testing and comparing the source and filament extraction methods to allow the researchers who need to perform such extractions to choose the most accurate algorithm for their projects. Although the benchmark was designed to resemble the Herschel observations of star-forming regions, the images are suitable for testing and evaluating extraction methods for any astronomical projects and applications. It consists of the complex fluctuating background cloud, the long dense filament, and many starless and protostellar cores with wide ranges of sizes, masses, and intensity profiles, computed with a radiative transfer code. A separate paper (Men’shchikov 2021) presents a series of the multiwavelength source extractions with getsf using three variants of the new benchmark with increasing complexity levels and compares their results with those produced by getold. All benchmark images, the truth catalogs containing the model parameters, and the reference extraction catalogs obtained by the author with getsf are available on its website.

The new extraction method can be used to conduct consistent and comparable studies of sources and filaments in various projects: getsf is designed to work for all images with nonzero background or noise fluctuations, where most pixels carry nonzero measured signal. The method is not limited to any particular area of astronomical research nor to the type of the telescopes or instruments used, as demonstrated by its applications to the images obtained with XMM-Newton, GALEX, Hubble, Spitzer, Herschel, APEX, and ALMA. Although no finite numbers of specific examples can prove that getsf is universally applicable, they confirm a remarkably wide applicability of the method.

---

Acknowledgments

This study used the cfitsio library (Pence 1999), developed at HEASARC NASA (USA), saoimage ds9 (Joye & Mandel 2003) and wcstools (Mink et al. 2002), developed at the Smithsonian Astrophysical Observatory (USA), and the stilts library (by Mark Taylor), developed at Bristol University (UK). The plot utility and ps12d library, used in this work to draw figures directly in the PostScript language, were written by the author using the psplot library (by Kevin E. Kohler), developed at Nova Southeastern University Oceanographic Center (USA), and the plotting subroutines from the MHD code azeus (Ramsey et al. 2012), developed by David Clarke and the author at Saint Mary’s University (Canada). This work used observations obtained with XMM-Newton, an ESA science mission with instruments and contributions directly funded by ESA Member States and NASA. This work used observations made with the Spitzer Space Telescope, which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under a contract with NASA. This work used observations made with the NASA/ESA Hubble Space Telescope, and obtained from the Hubble Legacy Archive, which is a collaboration between the Space Telescope Science Institute (STScI/NASA), the Space Telescope European Coordinating Facility (ST-ECF/ESA) and the Canadian Astronomy Data Centre (CADC/NRC/CSA). This paper used the SCUBA-2 data obtained at JCMT under program MJLSG31. The James Clerk Maxwell Telescope is operated by the East Asian Observatory on behalf of The National Astronomical Observatory of Japan; Academia Sinica Institute of Astronomy and Astrophysics; the Korea Astronomy and Space Science Institute; Center for Astronomical Mega-Science (as well as the National Key R&D Program of China with No. 2017YFA0402700). Additional funding support is provided by the Science and Technology Facilities Council of the United Kingdom and participating universities and organizations in the United Kingdom and Canada. Additional funds for the construction of SCUBA-2 were provided by the Canada Foundation for Innovation. This paper used the following ALMA data: ADS/JAO.ALMA#2013.1.01365.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), NSC and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. The simulated surface density background ℬ was derived from a synthetic scale-free background image created by Ph. André. A large set of images, used for testing and validation of getsf, includes those obtained in the Herschel Gould Belt Survey (http://gouldbelt-herschel.cea.fr) (HGBS, PI Ph. André), HOBYS (http://hobys-herschel.cea.fr) (PIs F. Motte, A. Zavagno, S. Bontemps), and ALMA-IMF (PIs F. Motte, A. Ginsburg, F. Louvet, P. Sanhoueza). HGBS and HOBYS are the Herschel Key Projects jointly carried out by SPIRE Specialist Astronomy Group 3 (SAG3), scientists of several institutes in the PACS Consortium (e.g., CEA Saclay, INAF-IAPS Rome, LAM/OAMP Marseille), and scientists of the Herschel Science Center (HSC). The author appreciates the valuable feedback, received from G. Zhang, F. Louvet, and N. Kumar, on the getsf extractions in the X-shaped nebula, MHD simulations, and Mon R2, respectively. The author is grateful to A. Zavagno, T. Nony, Y. Shimajiri, Ph. André, D. Arzoumanian, and P. Palmeirim for their comments on the manuscript.

References

Appendix A: Inaccuracies of the derived surface densities and temperatures

The algorithms described in Sect. 3.1.2 imply that the 160, 250, 350, and 500 μm images have an accurate (consistent) intensity calibration. When we assume that the calibration inaccuracies can be described by constant wavelength-dependent offsets, simple consistency checks and corrections can be made. Three independent estimates of low-resolution temperatures (𝒯_L1, 𝒯_L2, 𝒯_L3) are readily available from fitting the images in three pairs of wavebands (160−250, 250−350, and 350−500 μm) with a low resolution of O₅₀₀. If the median values of the three temperature maps differ by more than several percent, it would be necessary to adjust some of the offsets and estimate 𝒯_L{1|2|3} again. This iterative process is stopped when the three temperatures become consistent.

The higher-resolution images are obtained at the cost of significantly stronger noise and greater chances of distortions and spurious peaks. The quality of the resulting {𝒟|𝒯}_P from Eq. (5) strongly depends on the quality of the original short-wavelength images. Higher levels of noise or map-making artifacts in the 250 and 160 μm images would be amplified in the resulting maps in the process of fitting the spectral shapes Π_λ of pixels, which is likely to create significant small-scale distortions, predominantly in the pixels with strong line-of-sight temperature gradients that are usually located over the dense sources or filaments. The differential terms δ{𝒟|𝒯}_{3|2} that contain the higher-resolution information are increasingly less accurate because they are obtained from fitting of only three and two (noisier) images. It is very important to carefully inspect {𝒟|𝒯}_P to ensure that they are free of spurious small-scale structures before using them in any extraction. The hires images {𝒟|𝒯}_OH​ from Eqs. (8) and (11) are much less affected by the problems because they use the contributions δ{𝒟|𝒯}_{4|3|2} from all three variants of the fitted temperatures for each of the six resolutions of the original images.

The essential idea of the differential algorithm for improving the angular resolution of surface density was validated using the benchmark images (Sect. 2). The complete surface density 𝒟_C + 𝒟_S was first convolved to the resolutions of all Herschel wavebands (Sect. 2.3). The algorithm of Eq. (5), generalized to all six wavebands, was then applied to improve the lowest-resolution surface density using the unsharp masking of Eq. (6) and to successively recover each of the higher-resolution surface densities, all the way up to the highest adopted resolution O₇₀ = 8.4″, with a resulting maximum error below 0.5%. Although this is an excellent accuracy of the scheme, real-life applications of the method involve fitting of the spectral pixel shapes Π_λ, hence they inevitably suffer from larger inaccuracies (Fig. A.1).

The derived surface densities 𝒟_P and 𝒟_OH (Sect. 3.1.2) are not suitable for measuring dense structures, especially those with a central source of heating, because their inaccuracies in the pixels with strong line-of-sight temperature gradients are too large (e.g., Men’shchikov 2016). Comparisons with the true surface densities in Fig. A.1 show that the vast majority of the pixels outside bright sources are quite accurate, to better than 0.5%. However, the inaccuracies become much larger in the places that are occupied by the sources with steep gradients of the line-of-sight temperature. The starless cores and protostellar envelopes have markedly different radial temperature profiles, therefore the errors that are induced in the derived surface densities are also very dissimilar in both their sign and magnitude.

Appendix B: Single-scale spatial decomposition and standard deviations

Following the getold general approach, getsf employs successive unsharp masking to decompose the prepared original images ℐ_​λ (Sect. 3.1.1) into N_S single scales,

$$\begin{array}{c}\hfill {\mathcal{I}}_{\lambda j}={\mathcal{G}}_{j-1}\ast {\mathcal{I}}_{\lambda }-{\mathcal{G}}_j\ast {\mathcal{I}}_{\lambda },j=1,2,\cdots ,N_{\mathrm{S}},\end{array}$$(B.1)

where 𝒢_j are the circular Gaussian convolution kernels (𝒢₀ is to be regarded as the delta function) with progressively increasing half-maximum sizes,

$$\begin{array}{c}\hfill S_j=f{S}_{j-1},S_0=S_1{f}^{-1},S_{\mathrm{m}in}\le S_j\le S_{\mathrm{m}ax},\end{array}$$(B.2)

where f >  1 is the discretization factor (typically f ≈ 1.05) and the limiting scales of the decomposition range are

$$\begin{array}{c}\hfill \begin{array}{cc}{S}_{\mathrm{m}in}\hfill & =max(2\mathrm{\Delta },0.33{\min }_{\lambda }\left(O_{\lambda }\right)),\hfill \\ S_{\mathrm{m}ax}\hfill & ={\max }_{\lambda }(max(4X_{\lambda },4Y_{\lambda })),\hfill \end{array}\end{array}$$(B.3)

where Δ is the pixel size. The first image ℐ_​λ1 contains the contribution from all scales below S_min, whereas the last image ℐ_​λNS does not contain the signals from the scales above S_max, they are outside the range of scales being analyzed. The convolution is done with rescaling to conserve the total flux, hence the originals ℐ_​λ can be recovered by summation of the N_S scales and all remaining largest spatial scales,

$$\begin{array}{c}\hfill {\mathcal{I}}_{\lambda }={\displaystyle \sum _{j=1}^{N_{\mathrm{S}}}}{\mathcal{I}}_{\lambda j}+{\mathcal{G}}_{N_{\mathrm{S}}}\ast {\mathcal{I}}_{\lambda }.\end{array}$$(B.4)

The spatial decomposition is illustrated in Fig. B.1 using an example of a simple two-dimensional Gaussian shape 𝒫. As demonstrated in Papers I and II, the spatial decomposition has many useful properties. The filtered single-scale images contain signals from a relatively narrow range of spatial scales, and their properties resemble the Gaussian statistics much better than those of the originals, which are blends of all spatial scales. On the scales much smaller than the image size, the decomposed images are well described by the global value of the standard deviation σ_​λj. Significant departures from the Gaussian distribution in single scales above a certain threshold (e.g., I_λj ≳ 5σ_​λj) indicate the presence of the real structures. The decomposition highlights the structures of a specific width in the decomposed images on a matching scale. For example, a resolved isolated circular source with a half-maximum size H_λ has its maximum brightness in ℐ_​λj on the scale S_​j ≈ H_λ and a completely unresolved source produces the brightest signal on the smallest spatial scales S_​j ≲ O_λ.

Following the getold approach (Papers I and II), getsf employs an iterative algorithm to determine the single-scale σ_​λj over the entire usable area ℐ_​λjℳ_λ of the image to separate the real structures from other insignificant background or noise fluctuations. Before the iterations, the global σ_​λj0 and the threshold ϖ_λj0 = 3σ_​λj0 are computed over all pixels. At the first and all subsequent iterations (i = 1, 2, …, N_I), significant peaks and hollows with |I_λj|≥ϖ_{λji − 1} are masked. The absolute value is taken, because structures have both positive and negative counterparts in the decomposed images. Then getsf calculates a new (lower) σ_​λji value outside the masked areas and all structures with |I_λj|≥ϖ_λji are masked again. The iterations continue until the threshold converges to a stable value of ϖ_λji, with corrections δϖ_λji <  1%. The final single-scale standard deviation is obtained as σ_​λj = ϖ_λj/3 and its total value as ${\sigma }_{\lambda }^2={\sum }_j{\sigma }_{\lambda j}^2$. The constant 3, chosen empirically, provides both suitable values of the resulting σ_​λj values and good convergence of the iterations.

A notable difference with getold is that getsf does not need to correct the iterated thresholds using the higher-order statistical moments (skewness and kurtosis) because significant structures are detected in accurately flattened detection images (Sect. 3.4), which ensures that the majority of pixels resemble a normal distribution. Furthermore, precise σ_​λj values are of relatively minor importance for the separation of structural components because the separation is done in iterations and is based on the shapes that are removed from the single-scale slices (Sect. 3.2.2), not on the σ_​λj value itself.

Appendix C: Software suite

The method has been developed as a bash script getsf that executes a number of FORTRAN utilities, doing all numerical computations. Linux or macOS systems with the ifort or gfortran compilers can be used to install the code. For reading and writing images, getsf uses the cfitsio library (Pence 1999); for resampling and reprojecting images, it calls swarp (Bertin et al. 2002); for convolving images, it uses the fast Fourier transform routine rlft3 (Press et al. 1992); for determining the source coordinates α and δ, it applies xy2sky from wcstools (Mink et al. 2002); and for a colored screen output, it uses the highlight utility (by André Simon)¹⁶, if the latter is installed.

The following list of the getsf utilities and scripts explains their purpose and functions. They are quite useful for command-line image manipulations, even if there is no need in a source or filament extraction. Their usage information is displayed when a utility is run without any parameter. The utilities are sorted in the decreasing sequence of their general usability outside getsf.

modfits	modify an image or its header in various ways: math transformations; profiling an image along a line; image segmentation; filament skeletonization; removal of connected clusters of pixels; addition or removal of border areas; correction of saturated or bad pixel areas; conversion of intensity units; changes of the header keywords; etc.
operate	operate on two input images: addition, subtraction, multiplication, division; relative differencing; minimization or maximization; extension or expansion of masks; copying of an image header; computation of surface densities, temperatures, or intensities; etc.
imgstat	display and/or save image statistical quantities; produce mode-, mean-, or median-filtered images; compute images of standard deviations, skewness, kurtosis; etc.
fftconv	fast Fourier transform or convolve image with few predefined kernels or an external kernel image
fitfluxes	fit spectral shapes of source fluxes or image pixel intensities to derive masses or surface densities
convolve	convolve an image to a desired lower resolution
resample	resample and reproject an image with rotation
hires	high-resolution surface densities and temperatures
prepobs	convert observed images into the same pixel grid
installg	install getsf on a computer (macOS, Linux)
iospeed	test I/O speed of a hard drive for a specific image
readhead	display an image header or save selected keywords
cleanbg	interpolate background below source footprints
ellipses	overlay an image with ellipses of extracted sources
sfinder	detect sources in combined single-scale images
smeasure	measure and catalog properties of detected sources
fmeasure	measure and catalog properties of detected filaments
finalcat	produce the final catalogs of detected sources
expanda	expand masked areas of an image to its edges
extractx	extract all image extensions in separate images
splitcube	split a data cube into separate images

The code is automated, flexible, and user-friendly; it can be downloaded from the website¹⁷, the Astrophysics Source Code Library¹⁸, and it is also available from the author upon request. The website also contains a detailed User’s Guide and a complete validation extraction of sources and filaments in a small image for those who would like to verify that their installed getsf produces correct results.

All Figures