Multiscale, multiwavelength extraction of sources and filaments using separation of the structural components: getsf

High-quality astronomical images delivered by modern ground-based and space observatories demand adequate, reliable software for their analysis and accurate extraction of sources, filaments, and other structures, containing massive amounts of detailed information about the complex physical processes in space. The multiwavelength observations with highly variable angular resolutions across wavebands require extraction tools that preserve and use the invaluable high-resolution information. This paper presents getsf, a new method for extracting sources and filaments in astronomical images using separation of their structural components, designed to handle multiwavelength sets of images and very complex filamentary backgrounds. The method spatially decomposes the original images and separates the structural components of sources and filaments from each other and from their backgrounds, flattening their resulting images. It spatially decomposes the flattened components, combines them over wavelengths, detects the positions of sources and skeletons of filaments, and measures the detected sources and filaments. This paper presents a realistic multiwavelength set of simulated benchmark images that can serve as the standard benchmark problem to evaluate qualities of source- and filament-extraction methods. This paper describes hires, an improved algorithm for the derivation of high-resolution surface densities from the multiwavelength far-infrared Herschel images. The algorithm allows creating the surface densities with angular resolutions that reach 5.6 arcsec, when the 70 micron image is used. The codes getsf and hires are illustrated by their applications to a variety of images, from the X-ray domain to the millimeter wavelengths.


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 highquality data from space telescopes and large ground-based inter-Send offprint requests to: Alexander Men'shchikov ferometers, 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 ob-Article number, page 1 of 37 arXiv:2102.11565v2 [astro-ph.IM] 1 Mar 2021 A&A proofs: manuscript no. getsf served 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 starforming regions, which is used below to illustrate the method and in a separate paper (Men'shchikov 2021, submitted) 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 1 and a quality evaluation system (Men'shchikov 2021, submitted) 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 2013Men'shchikov , 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" 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 im-ages 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., A, B, C) and software names and numerical methods are typeset slanted (e.g., getsf ) to distinguish them from other emphasized words.

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 B λ , a long filament F λ , round sources S λ , and small-scale instrumental noise N λ : where C λ = B λ + F λ 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).

Simulated filamentary background
An image of the background surface density was computed from a purely synthetic scale-free background D A (cf. Paper I), with N H 2 ∼ 2.7 × 10 20 to 5 × 10 22 cm −2 that had uniform fluctuations across the entire image. To simulate complex astrophysical backgrounds with strongly nonuniform fluctuations (e.g., Könyves et al. 2015), D A was multiplied by a circular shape P 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 21 cm −2 was added to increase the minimum value. The surface densities of the resulting background cloud image D B (Fig. 2) are 1.5 × 10 21 to 4.8 × 10 22 cm −2 and the fluctuations differ by approximately two orders of magnitude. The total mass of the cloud is M B = 1.78 × 10 3 M . To simulate filamentary backgrounds, a long spiral filament was added to the background cloud D 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 D F has a crest value of N 0 = 10 23 cm −2 , 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. 2011Arzoumanian et al. , 2019, 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 H 2 (θ) ∝ θ −3 for θ Θ.
The filament mass M F = 3.04 × 10 3 M and length L F = 10.5 pc correspond to the linear density Λ F = 290 M pc −1 . The resulting surface densities D C = D B + D F of the filamentary cloud are in the range of 1.7 × 10 21 to 1.4 × 10 23 cm −2 (Fig. 2), and its total mass is M C = 4.82 × 10 3 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 The pixel values of the resulting temperature image T 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 C λ in all Herschel wavebands, assuming optically thin dust emission: 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 2 molecule, and m H is the hydrogen mass. The dust opacity was parameterized as a power law κ ν = κ 0 (ν/ν 0 ) β with κ 0 = 9.31 cm 2 g −1 (per gram of dust), λ 0 = 300 µm, and β = 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 6 cm −3 and coagulation time t = 10 5 yr. The opacity values at λ > 160 µm were replaced with a power law κ λ ∝ λ −2 , 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 −2 . 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(Motte et al. , 2001. Both types of cores were embedded in background spherical clouds with a uniform surface density of 3 × 10 21 cm −2 and outer radius of 1.4 × 10 5 AU (1000 or 0.68 pc). In an isotropic interstellar radiation field (Black 1994) with the strength parameter G 0 = 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 T 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.

Complete simulated images
Individual surface density images of the models of 828 starless and 91 protostellar cores were distributed in the dense areas (N H 2 ≥ 5 × 10 21 cm −2 ) of the filamentary cloud D C . They were added quasi-randomly, without overlapping, to the D C image at positions, where their peak surface density exceeded that of the cloud N H 2 value. An initial mass function (IMF)-like broken power-law mass function with a slope dN/dM of −1.3 for M ≤ 0.5 M and −2.3 for M > 0.5 M was used to determine the numbers of models per mass bin δlog 10 M ≈ 0.1 in each of the four populations. This resulted in the surface densities D S , the intensities S λ of sources (Fig. 3), and in the complete simulated images C λ + S λ .

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.

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.  Fig. 4. Images H λ of the simulated star-forming region, defined by Eq.
(1), shown at three selected wavelengths. The benchmark images are a superposition of four structural components: the background B λ , the filament F λ , the sources S λ , and the noise N λ . Two simpler variants of this benchmark are also available: without the filament and without the background. Square-root color mapping.

Original observed set of images
To prepare multiwavelength H λ for processing with getsf, it is necessary to convert them into the images I λ , 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 I λ 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 I λ coverage in the image processing, it is necessary to create masks M λ (with pixel values 1 or 0). With these masks, getsf can process only the good areas of I λ that have a mask value of 1. To facilitate the image preparation, getsf always creates default masks M λ = 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 2 that allows users to create a polygon over an image, convert the polygon into a mask, and save it in the FITS format.

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 powerlaw 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 lowerresolution 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 lowresolution images as differential terms, where the base surface density and temperature {D|T } 4 are derived by fitting the 160, 250, 350, and 500 µm images at the lowest resolution O 500 = 36.3 . The additional terms, containing the higher-resolution contributions, are produced by unsharp masking,  densities, with the assumptions and parameterizations of Eq. (4). It is required that the resolution of temperatures must not be higher than O λ , which excludes D O 350 2 and D O 500 {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 where λ H denotes the wavelength of the image I λ 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, where G O λ+ is the Gaussian kernel (regarded as the delta function at 500 µm), convolving D O λ {2|3|4} to a lower resolution of the next longer wavelength. For images at λ < 250 µm, only the positive values of δD 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 250 ) between I λ and the lower-resolution T {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 T {2|3|4} and overestimated values (within an order of magnitude) of peak surface densities at higher resolutions O λ < O 250 . This means that unsharp masking of the overestimated peaks could create negative annuli in δD O λ {2|3|4} and negative pixels in D O H . Fortunately, the surface densities are quite accurate outside the unresolved peaks (Appendix A).
A slight modification of Eq. (8) allows deriving the highresolution surface densities with an enhanced contrast of all un-resolved or slightly resolved structures, where the positive parts of the differential high-resolution terms from Eq. (9) are added to D O H . 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 T O H , consistent with the highresolution surface density D O H , is computed by numerically inverting the Planck function, with ν H = cλ −1 H , where c is the speed of light. The high-resolution images {D|T } 13 are shown in Fig. 5 along with the true simulated D C + D 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 500 , 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 −1 , 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 D O λ {2|3|4} that use all three temperatures T {2|3|4} with each original I λ .
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 T 2 at the resolution O 250 is removed from Eq. (7) and T {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 T 4 at the lowest resolution O 500 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 (Ap-pendix 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 highresolution surface densities and temperatures from Eqs. (8) and (11), it is straightforward to obtain such images: with the assumptions and parameterizations of Eq. (4). For example, the intensities J λ13 at 250, 350, and 500 µm would be sharper than I λ by the factors 1.3, 1.8, and 2.7, respectively. When the available original set of images I λ allows creation of D O H , it is advantageous to have it complement the original data set, handling it as an image I "observed" in a fictitious waveband . In the multiwavelength extractions with getsf, it may be recommended to use D O H 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 I λ . 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 I ≡ D 13 (Fig. 5), a total of N W = 7 wavelengths.

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 I λ 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 I λ 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 I ≡ D 13 (Fig. 5), the {X|Y} values are the same as those for the 250−500 µm images.

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 B λX is derived to separate the component of sources S λ , whereas the Y λ -scale background B λY is obtained to separate the component of filaments F λ . The backgrounds are collectively referred to as B λ{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 B λ{X|Y} are defined in getsf as the smooth intensity distributions on spatial scales S j larger than 4{X|Y} λ that remain in I λ 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.

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 I λ into a set of single-scale images I λ 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 I λ .

Separation of the structural components
The backgrounds B λ{X|Y} are computed by cutting small round peaks and elongated structures off the decomposed images I λ 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 I λ j by a number N L of intensity levels I λ jl , spaced by δ ln I λ 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 I λ j on that intensity level, producing various shapes of connected pixels, Relatively round source-like peaks in I λ j may be effectively distinguished from elongated structures by the number of connected pixels N λ jl that their shapes occupy in the slice I λ 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 I λ jl that are comparable to the area πS j 2 of the convolution kernel G j . In contrast to the round peaks, elongated shapes in I λ jl have greater lengths L λ than widths W λ ≈ S j , which means that the filamentary shapes in slices I λ jl extend over much larger areas than π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), Article number, page 9 of 37  where ∆ is the pixel size. Only simple and relatively straight filamentary shapes can be identified in I λ 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 I λ j using the three quantities described above. The shapes produced by sources in a slice I λ jl are not very elongated, not very sparse, and not very large. In contrast, the shapes produced by filaments in a slice I λ jl are elongated or sparse. Hence, these definitions for the source-like and filament-like shapes are written as 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 P in single scales P 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, 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) 3 , previously employed by getold to detect sources and filaments (Papers I and II). Deriving the background B λX of sources, getsf decomposes I λ and removes all source-like shapes from I λ jl , according to their definition in Eq. (15), in an iterative procedure (Sect. 3.2.3). Deriving the background B λY of filaments, getsf decomposes B λX and removes all filamentlike shapes from B λX jl , 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.

Reconstruction of the backgrounds
When we denote with B λ{X|Y} jlC either of the single-scale background slices B λX jlC or B λY jlC after the shape removal, the backgrounds on scale j are reassembled from the clipped slices as   (Fig. 2), convolved to the same resolution. The filament is heavily blended with itself in the central area, therefore its background is systematically underestimated there (lower right). Square-root color mapping, except in the right panels, which show linear mapping.
To properly reconstruct the complete backgrounds B λ{X|Y} from B λ{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, 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, where J λ{X|Y} is the number of the largest spatial scales 4{X|Y} λ for the backgrounds B λ{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, 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, 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, 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 The original images can be recovered by summing the three separated components: The positive parts of the small-scale background fluctuations and instrumental noise are contained in the component S λ , hence the component F λ appears fairly smooth (Fig. 8).

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 B λ{X|Y} greatly simplifies the original images, it does not reduce the strong variations of the smaller-scale fluctuation levels across {S|F } λ . 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 com-plete 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 Q λ and R λ that are derived by getsf from the images U λ and V λ of the standard deviations computed in the structural components with a circular sliding window of a radius O λ , where S λR and F λR are the regularized images S λ and F λ , 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 λ , This is done to improve the quality of {U|V} λ for further processing because the structural components from by Eq. (23) 27) and its standard deviations sd O λ (F λR R −1 λ ), which are much flatter (outside the filament) across the image. Square-root color mapping, except in the right panels, which show logarithmic mapping.

Decomposition of the standard deviations
The advantages of the spatial decomposition (Appendix B) apply also to the standard deviations {U|V} λ . The getsf method produces the single-scales {U|V} λ 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 I λ in Sect. 3.2.1.

Removal of the structural features
The {U|V} λ images sample local fluctuations and intensity gradients, revealing all sources and filaments present in I λ (Figs. 9, 10). To produce the corresponding flattening images, it is necessary to remove all such features from {U|V} λ , hence to determine their {X|Y} λ -scale backgrounds. Deriving the latter, getsf creates single-scale slices {U|V} λ jl , in a complete analogy with I λ jl in Sect. 3.2.2, and clips from them all source-and filamentlike shapes according to their definitions in Eq. (15). The reconstructed backgrounds Q λX and R λ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 {Q|R} λ{X|Y} that could create spurious structures if the images were used to flatten the structural components. To avoid this, the final flattening images Q λ and R λ (Figs. 9, 10) are obtained by median filtering the background in circular sliding windows of radii 2O λ and 5O λ , respectively, This important step ensures that flattening would never produce spurious structures in the detection images.

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 filamentdetection images are flattened, that is, divided by the flattening images, The standard deviations sd O λ ({S|F } λR {Q|R} −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.

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 I λ , separating two distinct structural components and creating the independent flat detection images {S|F } λD . In contrast to the originals, the flat images are suitable for the detection techniques that apply a threshold value for the entire image.

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 and Sects. 3.2.1, 3.3.1), getsf decomposes the detection images S λD and F λD into single scales {S|F } λD j and estimates the corresponding standard deviations σ λS j and σ λF j (Appendix B) that are necessary for separating significant structures from all other fluctuations. The decomposed components S λD j and F λD j are shown in Figs. 11 and 12 after they were cleaned and combined over wavebands.

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, where λS j = 5σ λS j and λF j = 2σ λF j . The filament threshold is significantly lower than that for sources because getsf additionally cleans F λD jC of the residual source-like clusters of connected pixels according to their definition in Eq. (15). The resulting clean images {S|F } λD jC (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.

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 S λD jC and F λD jC 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: where N {S|F} ≤ N W is the number of the wavebands chosen to be used in the combination, Z λ{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, 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 {S|F } D jC 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, where w λ is the weight that enhances the contribution of the images with higher angular resolutions, 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 S λD jC 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.

Detection of sources in the combined images
Sources are detected in S D jC 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 S D jC 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 S D jC and tracks their evolution from small to large scales, until they disappear or merge with a nearby brighter source.
To detect sources, getsf slices S D jC by a number N L of intensity levels I jl , spaced by δ ln I 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, S D jCl = max S D jC , I jl , l = 1, 2, . . . , N L .  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 filamentlike shapes. The resulting single-scale segmentation images of sources sets 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 j F (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 φ n H 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 = φ n H 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.

Detection of filaments in the combined images
Filaments are detected in F D jC 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 halfmaximum 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 ) archive 4 , and the hires surface densities D 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 understand-4 http://gouldbelt-herschel.cea.fr/archives ing 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 F D jC 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 K 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  The flattened components F D derived from the hires surface densities D 13 obtained from Eq. (8) using the Herschel 160, 250, 350, and 500 µm images are shown. The minimum scales of 36 (left column) correspond to 2.8 times the angular resolution, whereas the maximum scales (right column) correspond to 0.3 pc at the adopted distances of the regions (140, 260, and 460 pc, respectively). Intermediate scales between the two extremes are displayed in the middle column. The images were cleaned using the default threshold F j = 2σ F j . Overlaid on the filaments are their skeletons obtained from the images using the Hilditch algorithm (Sect. 3.4.5). The observed filaments are heavily substructured, and their appearance, detected skeletons, and measured properties depend strongly on spatial scales. Logarithmic color mapping.
Article number, page 17 of 37 A&A proofs: manuscript no. getsf further accumulated over a limited range of scales to produce a set of N K skeletons tracing the filamentary structures of various widths, where J − and J + are the numbers of the smallest and the largest scales, S J − = 2 −1/2 S k and S J + = 2 +1/2 S k , in the accumulated skeleton K k . The scale-dependent skeletons K k sample the following scales: where the scale S 1 =Ō is defined by Eq. (32) as the average angular resolution over the wavebands combined in F D jC (Sect. 3.4.3), and S N K = 4 max λ (Y λ ) is the largest spatial scale for the filament detection. Each pixel of the accumulated skeleton K 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 K j contributes to K k in that pixel. Depending on the filament intensity at the skeleton pixel, the significance range is 1 ≤ ξ < ∼ ln 2 (ln f ) −1 (≈ 14, assuming f ≈ 1.05, Appendix B). The algorithm automatically creates the final onepixel-wide skeletons by thresholding: K kξ = max (K k , ξ) with a default ξ = 2, and applying the Hilditch algorithm to the resulting shapes. Segmentation images of the skeletons K kξ are computed using the tintfill algorithm, which sets all pixels belonging to a filament to its number.

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 S λX (Sect. 3.2.3) because the subtracted background B λX contains substantial source residuals at low intensity levels (Fig. 8). The background B λ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 I λ 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 B Fλ of each source is determined by a linear interpolation of I λ 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 B Fλ is median filtered using a sliding window with a radius O λ , and the background-subtracted image of a source is then obtained as 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 j F , provided by the detection algorithm (Sect. 3.4.4). The corresponding initial footprint {A, B} Fn = φ n H 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 halfmaximum 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 halfmaximum 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 where the (empirical) exponential factor converts the average radius h λn into the equivalent-area radius (A λn B λ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 A. Men'shchikov: Multiscale, multiwavelength extraction of sources and filaments: getsf (2) with Θ = 5 and ζ = 1, transforming into a power law I ∝ θ −2 for θ Θ and filling up the entire image, its faint outer areas (I ∼ 0.2) are largely lost within the noise. The initial footprint factor φ n = 3 (Sect. 3.4.4) is too small for these power-law sources, hence background subtraction leaves a relatively bright pedestal containing a large amount of the source emission (middle). The footprint expansion algorithm (Sect. 3.4.6) enlarges φ n by a factor of 2.2 (right), which lowers the source background by a factor of 5, and as a result, increases the source flux by a factor of 2.7. The improved flux is still below the true value by a factor of 1.9 because the actual footprint is about three times larger. Square-root color mapping.
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 S λ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, 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, 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 I 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, 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 A&A proofs: manuscript no. getsf image I λ , 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), 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 ratios (S/Ns) Ω λn and Ψ λn , describing the detection and measurement properties of each extracted source, where j F is the footprinting scale (Sect. 3.4.4) and S λD j F n is the intensity at the source position in S λD j F (Sect. 3.4.1). The above quantities can be combined together to characterize the overall "goodness" of a source, 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, 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 2002), global flag f n , significance Ξ n , and goodness Γ n , n x n y n α n δ n f n Ξ n Γ n, followed (in the same line) by the measured quantities in each of the N W wavebands, 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 j F , 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: These empirical conditions, based on numerous test results obtained in various benchmarks (Pouteau et al., in prep.;Men'shchikov 2021, submitted), 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.

Measurements of the filaments
Filaments are measured in their background-subtracted F λY , derived in Sect. 3.2.3. When the maximum size Y λ of the filaments Article number, page 20 of 37 of interest is estimated sufficiently accurately (Sect. 3.1.3), their background B λ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 halfmaximum 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 threedimensional 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 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 F λY within a filament footprint, assuming that the image is obtained from surface densities, where M {α|β}n are the one-sided mass estimates, from which the average mass M n between the two sides is obtained. The onesided footprints Υ {α|β}n used in Eq. (45) 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 5 of the entire filament are readily obtained, together with the average linear densityΛ 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, 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 Eq. (46) and Eq. (47) are expected to be similar to each other for the well-behaved filaments.

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  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, submitted) 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 6 .

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 7 , 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 S λ , filaments F λ , and their backgrounds B λ{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 S λ . However, the images of standard deviations show that the flat source detection image S λ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. (44). 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.

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 8 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 S λ , filaments F λ , and their backgrounds B λ{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 S λ component effectively equalized the fluctuations across the detection image S λD , improving the extraction results.
The source catalog contains measurements of 1169 sources, 1130 of which are selected as acceptably good by Eq. (44). 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 F λ 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).

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 9 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 S λ , filaments F λ , and backgrounds B λ{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. (44). 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 F λ component comprises 100 skeletons representing its simple, non-branching segments (Sect. 3.4.5).

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 10 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 S λ and filaments F λ from their backgrounds B λ{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 S λ and F λ . The component S λ of sources (Fig. 19) is completely free of the elongated structures. The standard deviations U λ 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 S λ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 to 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 S λ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. (44). 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).

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 11 and Eq. (8) were used to compute a 1.1 • × 1.1 • surface density image D 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 D 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 H 2 = 3.7 × 10 22 cm −2 , 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 S , F , and B {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. (44). 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 H 2 = 10 22 cm −2 , 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 H 2 = 2.5 × 10 22 cm −2 has an FWHM size of 57 , almost the same as the size A = 58 of L 1689B, measured by getsf. However, in the derived D 13 image, the same model has M = 0.66 M and N H 2 = 10 22 cm −2 , 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 0 = 3.8 × 10 21 cm −2 , length L = 0.8 pc, half-maximum width W = 0.14 pc (205 ), mass M = 15 M , and linear density Λ = 14 M pc −1 ; the values are little affected by the fitting inaccuracies.

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 12 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 S λ 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 B λY of filaments is fairly low, hence its subtraction enhanced the visibility of filaments in F λ only little. Nonuniform small-scale fluc-tuations in S λ were effectively equalized in the detection image S λD by the flattening algorithm.
The source catalog contains measurements of 124 sources, 91 of which are selected as acceptably good by Eq. (44). 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 λ .

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 13 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 smallscale 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 S λ in that area from the filaments F λ and their backgrounds B λ{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 S λ , 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 U λ reveal imprints of the five overlapping scans from the observations. The flattening algorithm of getsf effectively equalizes them and creates the flat detection images {S|F } λ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. (44); 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 F λ 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 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 S λ and filaments F λ confirms that most sources are concentrated on (or near) the faint continuous filaments. Almost the entire background B λ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 U λ . 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 S λ and F λ . The flattening algorithm equalizes the fluctuation levels very effectively, providing reliable detection of sources in the flat {S|F } λ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 W 43-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. (44). 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.

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 singlescale 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.

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 Xray 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 2 to 2000 2 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 hours 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 2 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 ∼ 15000 do not pose any problems to getsf, but substantially larger numbers of detected sources require very large memory and long execution time.

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 14 , from the Astrophysics Source Code Library 15 , 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 sur-  Fig. 22. Application of getsf to the JCMT λ = 850 µm image (14.6 resolution) of the Orion A star-forming cloud, adopting {X|Y} λ = 30 . The top row shows the original image I λ and the backgrounds B λ{X|Y} of sources and filaments. The middle row shows the component S λ , the footprint ellipses of 212 acceptably good sources on S λD , and the component F λD with 199 skeletons K k2 corresponding to the scales S k ≈ 39 . The bottom row shows the standard deviations U λ in the regularized component S λR , the flattening image Q λ , and the standard deviations in the flattened component S λR Q −1 λ . Some skeletons may only appear to have branches because they were widened for this presentation. Intensities (in MJy sr −1 ) are limited in range with logarithmic color mapping, except for Q λ , which is shown with linear mapping.
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 (T L1 , T L2 , T 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 500 . 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 T 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 {D|T } 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 δ{D|T } {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 {D|T } P to ensure that they are free of spurious small-scale structures before using them in any extraction. The hires images {D|T } O H from Eqs. (8) and (11) are much less affected by the problems because they use the contributions δ{D|T } {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 D C + D 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 70 = 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 D P and D O H (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-ofsight 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 I λ (Sect. 3.1.1) into N S single scales, I λ j = G j−1 * I λ − G j * I λ , j = 1, 2, . . . , N S , (B.1) where G j are the circular Gaussian convolution kernels (G 0 is to be regarded as the delta function) with progressively increasing half-maximum sizes, where f > 1 is the discretization factor (typically f ≈ 1.05) and the limiting scales of the decomposition range are S min = max (2∆, 0.33 min λ (O λ )) , S max = max λ (max (4X λ , 4Y λ )) , where ∆ is the pixel size. The first image I λ1 contains the contribution from all scales below S min , whereas the last image I λN S 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 I λ can be recovered by summation of the N S scales and all remaining largest spatial scales, The spatial decomposition is illustrated in Fig. B.1 using an example of a simple two-dimensional Gaussian shape P. 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 I λ 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 I λ j M λ 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 σ λ 2 = j σ λ j 2 . The constant 3, chosen empirically, provides both suitable values of the resulting σ λ j values and good convergence of the iterations.  At the highest resolution of 8 , the derived images are the most accurate, with the exception of the unresolved protostellar peak surface densities (Fig. 3), which become strongly overestimated (up to a factor of ∼ 10) because the temperatures T {2|3|4} along the lines of sight with large temperature gradients are underestimated. The range of displayed values is reduced for better visibility. Linear color mapping.
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. 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 19 , the Astrophysics Source Code Library 20 , 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. 19 http://irfu.cea.fr/Pisp/alexander.menshchikov/ 20 https://ascl.net/2012.001