© ESO, 2017

1. Introduction

Recent high-resolution and high-sensitivity observations of large areas of Galactic star-forming regions with the Herschel Space Observatory (Pilbratt et al. 2010) revealed dense interstellar clouds as extremely complex. Multiwavelength images of the numerous regions observed within the Herschel Gould Belt (André et al. 2010) and HOBYS (Motte et al. 2010) surveys display bright and strongly variable backgrounds (on all spatial scales), blended with the omnipresent filamentary structures and many hundreds of (crowded) sources of different nature (e.g., Men’shchikov et al. 2010; André et al. 2014). The observed extreme complexity reflects that of the astrophysical reality, which is further exacerbated by its projection onto the sky plane. To analyze and understand the reality, it is necessary to separate the blended components: sources, filaments, and background (and/or noise).

The enormous complexity of observed images presents serious problems for the fully automated source extraction methods that apply global thresholds (e.g., sextractor, cutex, getsources, getfilaments; Bertin & Arnouts 1996; Molinari et al. 2011; Men’shchikov et al. 2012; Men’shchikov 2013) in order to separate sources from background and noise. The use of thresholds in terms of the background brightness or fluctuation levels that are computed for an entire image is based on the assumption that these thresholds are approximately uniform across the image. The observations of complex backgrounds demonstrate, however, that intensities and their fluctuations both often vary by several orders of magnitude in different parts of large fields. It is easy to see that any single threshold value would be inappropriate for such images, either producing many spurious sources or leaving many sources undetected. For a complete and reliable source detection, it would be highly desirable to simplify the complex observed images by removing signals from irrelevant spatial scales and equalizing the background fluctuation levels (outside sources) over the images. A solution of this non-trivial problem might require iterations, as the locations of sources are unknown before source extraction.

An attempt to solve this problem was made in two previous papers (Men’shchikov et al. 2012; Men’shchikov 2013, hereafter referred to as Paper I and Paper II), where the multiscale, multiwavelength source and filament extraction methods getsources and getfilaments were presented. The algorithms employ two iterations for a complete source extraction. An initial (somewhat deeper) extraction aims to detect all sources (possibly including spurious sources) and estimate their sizes. After extracting and removing the sources and filaments, the method attempts to equalize fluctuation levels of the remaining background. A second (and final) extraction uses the flattened detection images to identify sources more reliably, reducing the chances of detecting spurious sources. Although this approach works well for some types of images, it is not universal and therefore not fully satisfactory. In the most complex images, the initial extraction tends to produce large areas of many overlapping sources. Extremely complex observed backgrounds, in combination with large crowded areas of overlapping sources, could make the derived background inexact and hence render that approach inaccurate in some applications.

This paper presents a much better solution of the problem of background derivation and image flattening that does not require any prior source extraction or knowledge of the locations and sizes of sources. The new method, called getimages, is based on a simple and straightforward transformation of the observed images using median filtering.

A well-known technique (Tukey 1977), median filtering has been widely used during the past four decades in image processing, especially to reduce noise or small-scale artifacts. The idea of this image transformation is to replace the value of each pixel with a median value computed over all pixels in a window centered on the pixel that is being modified. The size of the sliding window is a free parameter of this technique. For the purpose of reducing sharp pixel-to-pixel intensity jumps, even a small 3 × 3 pixel window may be sufficient, whereas larger window sizes are necessary to remove wider features.

Two-dimensional median filtering has the very useful property of efficiently truncating any peaks – those produced by instrumental noise, intensity fluctuations of a background molecular cloud, filamentary structures, or astrophysical emission sources of any type. Finding a universal way of removing sources or other structures of all sizes of interest would essentially mean estimating their background, a major step toward designing a reliable method of flattening background-subtracted images to ensure uniform fluctuation levels in all their parts. These are the ideas that initiated the development of getimages.

As in Papers I and II, images are represented by capital calligraphic characters (e.g., ) to distinguish them from other parameters. The median filtering operator using a window of radius R is denoted as mf_R(ℐ), and the standard deviations operator using a circular n-pixels window is denoted as sd_n(ℐ). The new background derivation and image flattening method is described in Sect. 2 and discussed in Sect. 3, the conclusions are outlined in Sect. 4, and further details are presented in Appendices A–C.

2. Background estimation and flattening

2.1. Median filtering as a structure removal operator

The point-spread functions (PSFs, beams) of modern telescopes are usually represented by two-dimensional Gaussians down to percent levels, below which the beams become more complicated¹. Observations show mostly unresolved or slightly resolved structures with Gaussian profiles. However, some bright and well-resolved structures are better approximated by a Gaussian core with power-law wings. In addition to the simple Gaussian shapes, this paper therefore also considers profiles with the functional form used in Papers I and II: $\mathit{I}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{=}{\mathit{I}}_{\mathrm{P}}{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{f}\mathrm{(}\mathit{\zeta }\mathrm{)}\hspace{0.17em}\mathrm{(}\mathit{\theta }\mathit{/}\mathrm{\Theta }{\mathrm{)}}^{\mathrm{2}}\mathrm{)}}^{\mathrm{-}\mathit{\zeta }}\mathit{,}$(1)where I_P is the structure peak intensity, θ the angular distance, Θ the structure half-width at half-maximum, ζ a power-law exponent, and f(ζ) = (2^{1 /ζ}−1) adjusts the profile to have the same Θ for all values of ζ. This function has an almost Gaussian shape in its core, transforming into a power-law profile I(θ) ∝ θ^{− 2 ζ} for θ ≫ Θ. For simplicity, ζ = 1 is fixed in this section, and the term “power law” refers to the shapes with I(θ) ∝ θ^-2 at large θ. The Gaussian and power-law shapes can represent both starless cores and protostellar envelopes (cf. Appendix A), depending on their signal-to-noise ratio (S/N) and on the degree to which they are resolved.

Fig. 1 Median filtering applied to Gaussian sources (upper row), power-law sources (second row), Gaussian filaments ℱG (third row), and power-law filaments ℱP (bottom row) with sizes H (FWHM) of 8′′ and 16′′ (left) and to the same images after addition of random Gaussian noise with σ = 0.33 (middle) and background ℬ modeled as a large Gaussian with FWHM of 512′′ (right). The radii R of circular sliding windows given by the indices on the curve labels correspond to the relative radii W = R/H = { 1,2,4,8 }. Original intensity profiles of , , ℱG, and ℱP are shown in gray, filtered profiles of the sources and filaments with H = 8 and 16′′ are plotted in red and blue, whereas profiles of and ℬ are colored in cyan and green, respectively. Large sliding windows (W> 2) truncate Gaussian shapes so efficiently (cf. Table 1) that some of the annotated curves become invisible (fall entirely below the lower edge of the plots).

Efficiently reducing noise and various artifacts, median filtering also erases other types of structures, such as sources and filaments. The steeper their intensity distributions, the better they are removed with smaller sliding windows. This property of median filtering is demonstrated below using simulated images of sources and linear filaments with Gaussian and power-law intensity distributions: , , ℱ_G, and ℱ_P. The sources and filaments have sizes H = 8 and 16′′(full-width at half-maximum, FWHM) and peak intensities I_P = 100 (in arbitrary units), the filaments extend across the entire image. The simulated images have dimensions of 803 × 803 with 0.67″ pixels. Another variant of the images adds random Gaussian noise convolved to a 2′′resolution with a standard deviation σ = 0.33 (zero mean). The third variant includes a large-scale background ℬ modeled as a wide Gaussian with a peak value of 50 and a half-maximum width of 512′′.

Results of median filtering of all simulated structures are presented in Fig. 1 by the intensity profiles through their peaks (see Appendix B for the corresponding images). The resulting truncation factors, defined as the ratio f_T = I_P/I_F of the original and filtered peak intensities, are summarized in Table 1. Median filtering was parameterized by W = R/H, the window size in units of the structure half-maximum size, increasing by factor of two (1, 2, 4, etc.). The original structures of different sizes (H = 8 and 16′′) are truncated by the same factors depending only on W. For the same window, Gaussian shapes are erased much more efficiently than power-law shapes and sources are truncated much deeper than filaments. This is a simple consequence of the fact that the more extended or elongated structures fill sliding windows to a much higher degree. The presence of noise and background does not affect the truncation factors until they increase to f_T ~ 100, whereas for larger windows, noise and background tend to reduce the factors.

This demonstrates that median filtering efficiently removes Gaussian sources and filaments in small windows (W ≈ 2 and 4) with very low residuals of ≲ 1%. More extended intensity distributions of the power-law sources and filaments require larger windows (W ≈ 8 and 16) for their removal to the same high accuracy. Median filtering with a smaller window (W = 4) erases power-law sources and filaments to reasonably good levels of 3 and 10%, respectively (Table 1).

If the sources with size H are cleanly erased, the filtered image closely approximates background ℬ_H for those sources and ℐ − ℬ_H becomes their background-subtracted image. In the easiest case of isolated sources on a simple background, their accurate sizes and fluxes might be measured directly in ℐ − ℬ_H. The problem is that in real-life images sources have different sizes and backgrounds are not simple.

2.2. Background derivation by median filtering

Astrophysical backgrounds are highly filamentary and strongly variable on all spatial scales. Since median filtering truncates any peaks, in addition to sources, it also removes some intensity peaks that belong to the filaments and complex backgrounds. Backgrounds obtained by median filtering may therefore become underestimated in some places, hence the background-subtracted images may contain contributions from the filaments and backgrounds. Direct measurements of source sizes and fluxes in such images are not always feasible because they will be inaccurate for some sources. However, with drastically reduced contributions of filaments and bright fluctuating backgrounds, the background-subtracted images are always much simpler than the originals. There are great benefits in using them as detection images, in addition to the originals, which should be used as measurement images (cf. Papers I and II).

Images of star-forming regions may be considered as superpositions of sources, backgrounds, filaments, and noise. To illustrate the new algorithm, images of a simulated star-forming region were computed at Herschel wavelengths, for simplicity, consisting of only the sources and background: . The images are essentially identical to those described in Paper I (Appendix C) without instrumental noise; they have dimensions of 1900 × 1900 with 2′′pixels. The purely synthetic scale-free cirrus background images were made consistent with a planar temperature gradient between 15 and 20 K along one image diagonal (from lower right to upper left, as in Sect. 2.8 of Paper I). Examples of the simulated images ℐ_λ at wavelengths λ of 70, 160, and 350 μm, convolved to the corresponding angular resolutions O_λ of Herschel (8.4,13.5, and 24.9′′FWHM), are shown in Fig. 2 (upper panels).

Astronomical images observed with a limited resolution contain both unresolved and resolved sources whose sizes H_λ ≥ O_λ are unknown before source extraction. Determining the background of such images using median filtering with a single sliding window radius R_λ may not be optimal, except when the goal is to extract sources of similar sizes. For example, one might be interested in detecting unresolved or slightly resolved sources with H_λ ≅ O_λ and use a sliding window with R_λ ≈ 4 O_λ (cf. Sect. 2.1). When extracting sources with a wide range of sizes, it is conceivable to use a sliding window tailored to the largest source. This may not always work, however, as median filtering with large windows spreads bright emission down to the low-background areas, leading to an overestimated background and undetected faint sources in those areas of a background-subtracted image.

Fig. 2 Background derivation by median filtering. Shown are the images ℐλ of a simulated star-forming region (upper row), median-filtered background $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathrm{ℬ}}_{\mathit{\lambda }}\end{array}$ (middle row), and background-subtracted images $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathrm{𝒮}}_{\mathit{\lambda }}\end{array}$ (bottom row) at selected Herschel wavelengths (after M = 30 iterations). The maximum source sizes Xλ of 25, 100, and 100′′ for the getimages method were estimated directly from ℐλ (cf. Sect. 3.2). In panel a, small holes are starless cores seen in absorption at 70 μm, and white emission peaks are protostellar sources. Intensities (in MJy/sr) are limited in range, and their color scaling is linear.

Fig. 3 Quality of the derived background presented in Fig. 2. Shown are the relative accuracies $\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathrm{ℬ}}_{\mathit{\lambda }}\end{array}\mathit{/}{\mathrm{ℬ}}_{\mathit{\lambda }}\hspace{0.17em}\mathrm{-}\hspace{0.17em}\mathrm{1}$ of the background $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathrm{ℬ}}_{\mathit{\lambda }}\end{array}$ obtained by median filtering with respect to the true background ℬλ. Only the relevant pixels within the model sources are shown for values limited by the range [−0.3,0.3 ] with linear color coding. In panel a, the minima visible in the centers of some footprints are caused by the absorption of radiation in the central parts of starless cores; the two deepest minima are at levels of −0.7 and −0.6. In panels b and c, the maxima are at levels of 0.8 and 0.7, respectively.

The getimages method provides a simple and universal procedure to remove all structures of any width below an arbitrary maximum size X_λ (FWHM) and to derive background $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathrm{ℬ}}_{\mathit{\lambda }}\end{array}$ for the size range O_λ ≤ S_λ ≲ X_λ. Given the maximum size X_λ, the algorithm defines N sliding windows with radii R_{λ 1},R_{λ 2},...,R_λN, such that R_{λ 1} = 2 O_λ, R_λj = f_WR_λj−1, and R_λN = 4 X_λ, where the factor f_W> 1 must be small enough to sample the range of sizes. Median filtering of the original images ℐ_λ is repeated using N sliding windows and minimizing the resulting set of images: $\mathrm{ℬ̂}\mathit{\lambda }\mathrm{=}\begin{array}{c}\min \end{array}\mathrm{\{}\hspace{0.17em}{\mathrm{mf}}_{{\mathit{R}}_{\mathit{\lambda }\hspace{0.17em}\mathit{j}}}\mathrm{(}{\mathrm{ℐ}}_{\mathit{\lambda }}\mathrm{)}{\hspace{0.17em}}^{\mathrm{\}}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{\left(}\mathit{j}\mathrm{=}\mathrm{1}\mathit{,}\mathrm{2}\mathit{,}\mathit{...}\mathit{,N}\mathrm{\right)}\mathit{.}$(2)Finally, the median filtering with the largest window is repeated to improve the smoothness of the resulting image: $\begin{array}{c}\mathrm{˜}\\ {\mathrm{ℬ}}_{\mathit{\lambda }}\end{array}\mathrm{=}\begin{array}{c}\min \end{array}\left\{\hspace{0.17em}\mathrm{ℬ̂}\mathit{\lambda }\mathit{,}\hspace{0.17em}{\mathrm{mf}}_{{\mathit{R}}_{\mathit{\lambda }\hspace{0.17em}\mathit{N}}}\mathrm{\left(}\mathrm{ℬ̂}\mathit{\lambda }\mathrm{\right)}\hspace{0.17em}\right\}\mathit{.}$(3)The smallest window with R_{λ 1} = 2 O_λ completely removes Gaussian sources, and the largest window with R_λN = 4 X_λ truncates Gaussian and power-law sources and filaments sufficiently well (cf. Sect. 2.1, Table 1). Practical details of the definition of X_λ are discussed in Sect. 3.2.

A reasonably low value f_W = 2^1/2 of the discretization factor is adopted in getimages by default. Although lower values (f_W → 1) ensure that the size range of the structures of interest is better resolved, in practice, they slow down the procedure without any noticeable gain in accuracy. On the other hand, large factors (f_W ≳ 3) lead to a poorly sampled size range and thus tend to give a less accurate background $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathrm{ℬ}}_{\mathit{\lambda }}\end{array}$ for some structures.

This procedure ensures that the background obtained with small windows for the structures with H_λ ≈ O_λ is preserved after median filtering with much larger windows that are suitable for significantly wider sources or filaments (H_λ ≫ O_λ). The resulting background image $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathrm{ℬ}}_{\mathit{\lambda }}\end{array}$ becomes appropriate for structures within the entire size range O_λ ≤ S_λ ≲ X_λ. The background-subtracted image is readily computed as $\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathrm{𝒮}}_{\mathit{\lambda }}\end{array}\hspace{0.17em}\mathrm{=}\hspace{0.17em}{\mathrm{ℐ}}_{\mathit{\lambda }}\hspace{0.17em}\mathrm{-}\hspace{0.17em}\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathrm{ℬ}}_{\mathit{\lambda }}\end{array}$, which is essentially the original image of the structures containing a small differential contribution ${\mathrm{ℬ}}_{\mathit{\lambda }}\hspace{0.17em}\mathrm{-}\hspace{0.17em}\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathrm{ℬ}}_{\mathit{\lambda }}\end{array}$ induced by the inaccuracies in the estimated background.

To significantly improve $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathrm{ℬ}}_{\mathit{\lambda }}\end{array}$ and reduce contamination of $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathrm{𝒮}}_{\mathit{\lambda }}\end{array}$ by residual background peaks, the getimages method employs iterations. In the iterative formulation, the median-filtering algorithm described by Eqs. (2) and (3) remains the same, with two substitutions: ${\mathrm{ℐ}}_{\mathit{\lambda }}\hspace{0.17em}\mathrm{\leftarrow }\hspace{0.17em}\begin{array}{c}\mathrm{˜}\\ {\mathrm{𝒮}}_{\mathit{\lambda }}^{\hspace{0.17em}\mathit{i}\mathrm{-}\mathrm{1}}\end{array}\mathit{,}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\begin{array}{c}\mathrm{˜}\\ {\mathrm{ℬ}}_{\mathit{\lambda }}\end{array}\hspace{0.17em}\mathrm{\leftarrow }\hspace{0.17em}\stackrel{}{{{\mathrm{ℬ}}_{\mathit{\lambda }}^{\hspace{0.17em}\mathit{i}}}_{\mathrm{˜}}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{(}\mathit{i}\mathrm{=}\mathrm{1}\mathit{,}\mathrm{2}\mathit{,}\mathit{...}\mathit{,M}\mathrm{)}\mathit{,}$(4)where $\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathrm{𝒮}}_{\mathit{\lambda }}^{\hspace{0.17em}\mathrm{0}}\end{array}\hspace{0.17em}\mathrm{=}\hspace{0.17em}{\mathrm{ℐ}}_{\mathit{\lambda }}$ and M is the number of iterations. To construct background-subtracted images exclusively for source detection, it would be sufficient to set M ≈ 5−10. To open the possibility of using the images also for source measurements, more iterations (M ≈ 20−30) may be necessary. The final improved images are computed as $\begin{array}{c}\mathrm{˜}\\ {\mathrm{ℬ}}_{\mathit{\lambda }}\end{array}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{M}}\stackrel{}{{{\mathrm{ℬ}}_{\mathit{\lambda }}^{\hspace{0.17em}\mathit{i}}}_{\mathrm{˜}}}\mathit{,}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\begin{array}{c}\mathrm{˜}\\ {\mathrm{𝒮}}_{\mathit{\lambda }}\end{array}\hspace{0.17em}\mathrm{=}\hspace{0.17em}{\mathrm{ℐ}}_{\mathit{\lambda }}\hspace{0.17em}\mathrm{-}\hspace{0.17em}\begin{array}{c}\mathrm{˜}\\ {\mathrm{ℬ}}_{\mathit{\lambda }}\end{array}\mathit{.}$(5)Experience shows that the iterative scheme effectively reduces the differential contribution ${\mathrm{ℬ}}_{\mathit{\lambda }}\hspace{0.17em}\mathrm{-}\hspace{0.17em}\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathrm{ℬ}}_{\mathit{\lambda }}\end{array}$ of background peaks in $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathrm{𝒮}}_{\mathit{\lambda }}\end{array}$.

Fig. 4 Flattening background-subtracted images (Fig. 2). Shown are the standard deviations ${\mathrm{𝒟}}_{\mathit{\lambda }}\hspace{0.17em}\mathrm{=}\hspace{0.17em}{\mathrm{sd}}_{\mathrm{9}\hspace{0.17em}}\mathrm{\left(}\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathrm{𝒮}}_{\mathit{\lambda }}\end{array}\mathrm{\right)}$ of small-scale fluctuations in $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathrm{𝒮}}_{\mathit{\lambda }}\end{array}$ (upper row), median-filtered scaling (flattening) images $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathrm{ℱ}}_{\mathit{\lambda }}\end{array}$ (middle row), and flattened detection images ${\mathrm{ℐ}}_{\mathit{\lambda }\mathrm{D}}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathrm{𝒮}}_{\mathit{\lambda }}\end{array}\mathit{/}\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathrm{ℱ}}_{\mathit{\lambda }}\end{array}$ (bottom row) at selected Herschel wavelengths. In panel b, thin ring-like structures reflect the off-center temperature peaks in starless cores illuminated by the interstellar radiation field. The maximum sizes Xλ of 25, 100, and 100′′ are the same as those used for background derivation.

Fig. 5 Quality of the flattened images presented in Fig. 4. Shown are the standard deviations sd9 (ℐλD) of small-scale fluctuations in the flattened images ℐλD. The fluctuations outside sources are very uniform, compared to ${\mathrm{𝒟}}_{\mathit{\lambda }}\hspace{0.17em}\mathrm{=}\hspace{0.17em}{\mathrm{sd}}_{\mathrm{9}\hspace{0.17em}}\mathrm{\left(}\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathrm{𝒮}}_{\mathit{\lambda }}\end{array}\mathrm{\right)}$ computed from the background-subtracted images.

The background derivation algorithm is illustrated in Fig. 2. The adopted maximum source sizes X_λ of 25, 100, and 100′′for the images at 70, 160, and 350 μm, respectively, were estimated manually from ℐ_λ according to Eq. (8) in Sect. 3.2. The median filtering procedure described by Eqs. (2)–(5) was applied to ℐ_λ and the background images $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathrm{ℬ}}_{\mathit{\lambda }}\end{array}$ were obtained after M = 30 iterations. With X₇₀ = 25′′, intensity variations of $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathrm{ℬ}}_{\mathrm{70}}\end{array}$ are preserved to smaller spatial scales in comparison to $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathrm{ℬ}}_{\mathrm{160}}\end{array}$ and $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathrm{ℬ}}_{\mathrm{350}}\end{array}$ derived using larger sliding windows to accommodate sources up to X_λ = 100′′. Background-subtracted images $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathrm{𝒮}}_{\mathit{\lambda }}\end{array}$ are much simpler and flatter, but they preserve all structures with sizes H_λ ≲ X_λ that exist in ℐ_λ. Revealing the sources much more clearly than the originals do, they contain a small additional contribution ${\mathrm{ℬ}}_{\mathit{\lambda }}\hspace{0.17em}\mathrm{-}\hspace{0.17em}\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathrm{ℬ}}_{\mathit{\lambda }}\end{array}$ from the median-filtered peaks of the true model background ℬ_λ. The background peaks remaining in $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathrm{𝒮}}_{\mathit{\lambda }}\end{array}$ are narrower at 70 μm than at 160 and 350 μm, according to the difference in X_λ. The background residuals are the reason why, in general, $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathrm{𝒮}}_{\mathit{\lambda }}\end{array}$ can only be used for detection, with the exception of simple and smooth extended backgrounds whose peaks have sizes H_λ ≫ X_λ.

The relative accuracies of $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathrm{ℬ}}_{\mathit{\lambda }}\end{array}$ for all sources are presented in Fig. 3, where only the relevant pixels belonging to the model sources are displayed. The accuracy of derived background depends on the source size and its position. For most of the sources, $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathrm{ℬ}}_{\mathit{\lambda }}\end{array}$ is estimated to within 10%. Negative errors are found in the places where the true fluctuating background has local peaks. Relatively large positive errors (red and white areas in Fig. 3) are observed only inside very extended isolated or overlapping sources, where wide areas of the true fluctuating background ℬ_λ may have deeper hollows. Such an overestimation is impossible to correct, as the true background under real sources is very uncertain (cf. Appendix B in Men’shchikov 2016). When averaged over an entire source, the background accuracy becomes much better. However, the accuracy of $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathrm{ℬ}}_{\mathit{\lambda }}\end{array}$ is irrelevant for detection images that are not meant to be used for measurements. Detection images must only preserve all sources, their positions and sizes, which is clearly the case for the results displayed in Fig. 2.

2.3. Flattening background-subtracted images

Observations demonstrate that background intensities and their fluctuations vary by orders of magnitude across images. They are especially variable in the short-wavelength images that are affected by relatively strong dust temperature deviations induced by the radiation from hot stars. Although background subtraction greatly simplifies ℐ_λ (Sect. 2.2), the removal of an average background $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathrm{ℬ}}_{\mathit{\lambda }}\end{array}$ on spatial scales larger than X_λ does not strongly reduce variable fluctuations across $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathrm{𝒮}}_{\mathit{\lambda }}\end{array}$ on scales ≲ X_λ.

A flattening procedure introduced in getsources (Paper I) attempted to equalize the fluctuation levels in different parts of detection images ℐ_λD by dividing the latter with flattening images obtained from small-scale background fluctuations. However, the original scheme required a complete preliminary source extraction in order to determine background by cutting off the extracted sources. Depending on the extraction quality in the original complex images, the two-step approach was not completely universal and hence not fully satisfactory. In contrast, the background derivation procedure introduced in Sect. 2.2 removes all sources automatically without any need of a prior iterative source extraction. The simple median-filtering scheme of Eqs. (2) and (3) leads to a new straightforward and accurate flattening procedure.

The standard deviations of $\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathrm{𝒮}}_{\mathit{\lambda }}\end{array}\hspace{0.17em}\mathrm{=}\hspace{0.17em}{\mathrm{𝒮}}_{\mathit{\lambda }}\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathrm{\left(}{\mathrm{ℬ}}_{\mathit{\lambda }}\hspace{0.17em}\mathrm{-}\hspace{0.17em}\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathrm{ℬ}}_{\mathit{\lambda }}\end{array}\mathrm{\right)}$ are computed in a small sliding window of 3 × 3 pixels², and the operation is denoted as ${\mathrm{𝒟}}_{\mathit{\lambda }}\hspace{0.17em}\mathrm{=}\hspace{0.17em}{\mathrm{sd}}_{\mathrm{9}\hspace{0.17em}}\mathrm{\left(}\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathrm{𝒮}}_{\mathit{\lambda }}\end{array}\mathrm{\right)}$. For illustration, at selected wavelengths of 70, 160, and 350 μm is presented in Fig. 4 (upper panels). Approximately, the operation can be written as a sum of the individual components: ${\mathrm{𝒟}}_{\mathit{\lambda }}\hspace{0.17em}\mathrm{\approx }\hspace{0.17em}{\mathrm{sd}}_{\mathrm{9}\hspace{0.17em}}\mathrm{\left(}{\mathrm{𝒮}}_{\mathit{\lambda }}\mathrm{\right)}\hspace{0.17em}\mathrm{+}\hspace{0.17em}{\mathrm{sd}}_{\mathrm{9}\hspace{0.17em}}\mathrm{\left(}{\mathrm{ℬ}}_{\mathit{\lambda }}\hspace{0.17em}\mathrm{-}\hspace{0.17em}\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathrm{ℬ}}_{\mathit{\lambda }}\end{array}\mathrm{\right)}$. This formulation is not precise (in pixels where ${\mathrm{𝒮}}_{\mathit{\lambda }}\hspace{0.17em}\mathrm{\approx }\hspace{0.17em}{\mathrm{ℬ}}_{\mathit{\lambda }}\hspace{0.17em}\mathrm{-}\hspace{0.17em}\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathrm{ℬ}}_{\mathit{\lambda }}\end{array}$) and used here only to highlight the fact that contains contributions of the fluctuations induced by both sources and background residuals (in general, also by filaments and noise). As a first step toward flattening, it is necessary to remove from the fluctuations produced by sources, that is, to determine its background ${\mathrm{sd}}_{\mathrm{9}\hspace{0.17em}}\mathrm{\left(}{\mathrm{ℬ}}_{\mathit{\lambda }}\hspace{0.17em}\mathrm{-}\hspace{0.17em}\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathrm{ℬ}}_{\mathit{\lambda }}\end{array}\mathrm{\right)}$.

The obvious similarity of the problems of deriving source-free images from ℐ_λ and makes applying the median-filtering algorithm described in Sect. 2.2 feasible. Using the same set of N sliding windows with radii R_{λ 1},R_{λ 2},...,R_λN, the algorithm median filters the standard deviations and minimizes the resulting set of images: $\mathrm{ℱ̂}\mathit{\lambda }\mathrm{=}\begin{array}{c}\min \end{array}\mathrm{\{}\hspace{0.17em}{\mathrm{mf}}_{{\mathit{R}}_{\mathit{\lambda }\hspace{0.17em}\mathit{j}}}\mathrm{(}{\mathrm{𝒟}}_{\mathit{\lambda }}\mathrm{)}{\hspace{0.17em}}^{\mathrm{\}}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{\left(}\mathit{j}\mathrm{=}\mathrm{1}\mathit{,}\mathrm{2}\mathit{,}\mathit{...}\mathit{,N}\mathrm{\right)}\mathit{.}$(6)To conclude the procedure, the median filtering with the largest window is repeated (twice) to smooth the resulting image: $\begin{array}{c}\mathrm{˜}\\ {\mathrm{ℱ}}_{\mathit{\lambda }}\end{array}\mathrm{=}\begin{array}{c}\min \end{array}\left\{\hspace{0.17em}\mathrm{ℱ̂}\mathit{\lambda }\mathit{,}\hspace{0.17em}{\mathrm{mf}}_{{\mathit{R}}_{\mathit{\lambda }\hspace{0.17em}\mathit{N}}}\mathrm{\left(}\mathrm{ℱ̂}\mathit{\lambda }\mathrm{\right)}\mathit{,}\hspace{0.17em}{\mathrm{mf}}_{{\mathit{R}}_{\mathit{\lambda }\hspace{0.17em}\mathit{N}}}({\mathrm{mf}}_{{\mathit{R}}_{\mathit{\lambda }\hspace{0.17em}\mathit{N}}}\mathrm{\left(}\mathrm{ℱ̂}\mathit{\lambda }\mathrm{\right)})\hspace{0.17em}\right\}\mathit{.}$(7)The above procedure ensures that background fluctuations obtained with small windows for sources or filaments with H_λ ≈ O_λ survive median filtering with much larger windows that are suitable for more extended structures. As a consequence, contributions of all sources with sizes O_λ ≤ S_λ ≲ X_λ are removed from the image $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathrm{ℱ}}_{\mathit{\lambda }}\end{array}$ of background fluctuations.

A flattened detection image is obtained by dividing $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathrm{𝒮}}_{\mathit{\lambda }}\end{array}$ by the flattening (scaling) image: ${\mathrm{ℐ}}_{\mathit{\lambda }\mathrm{D}}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathrm{𝒮}}_{\mathit{\lambda }}\end{array}\mathit{/}\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathrm{ℱ}}_{\mathit{\lambda }}\end{array}$. This procedure effectively equalizes fluctuation levels across the entire image ℐ_λD, while preserving the intensity distributions of sources.

The flattening algorithm is illustrated in Fig. 4. The adopted maximum sizes X_λ of 25, 100, and 100′′for the images at 70, 160, and 350 μm, respectively, are the same as those used in Sect. 2.2 to derive background $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathrm{ℬ}}_{\mathit{\lambda }}\end{array}$. The standard deviations of the background-subtracted $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathrm{𝒮}}_{\mathit{\lambda }}\end{array}$, computed in a 9-pixel sliding window, clearly show the sources and that the residual background fluctuations increase along one diagonal. The amplifying fluctuations are induced by the planar temperature gradient adopted in the synthetic background (Sect. 2.2). Such images are not well suited for a complete and reliable source extraction. It is easy to see that global thresholding methods are bound to either produce spurious sources or leave some faint sources undetected. The intensity gradient of fluctuations is also visible in the flattening images $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathrm{ℱ}}_{\mathit{\lambda }}\end{array}$ obtained using the median-filtering procedure described by Eqs. (6) and (7). By construction, the images ℐ_λD in Fig. 4 must have uniform fluctuation levels and they do appear much flatter than $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathrm{𝒮}}_{\mathit{\lambda }}\end{array}$ (Fig. 2).

Fig. 6 Application of getimages (with Xλ = 12′′ and M = 30) to real-life (Spitzer λ = 8 μm) observations of the Ophiuchus L 1688 star-forming region: a) original image ℐλ, b) estimated background $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathrm{ℬ}}_{\mathit{\lambda }}\end{array}$, c) background-subtracted image $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathrm{𝒮}}_{\mathit{\lambda }}\end{array}$, d) standard deviations , e) scaling image $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathrm{ℱ}}_{\mathit{\lambda }}\end{array}$, f) flattened detection image ℐλD, and g) standard deviations sd9 (ℐλD). Selected results for Xλ = 24′′ are also shown: h) scaling image $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathrm{ℱ}}_{\mathit{\lambda }}\end{array}$ and i) flattened image ℐλD. Intensities (in MJy/sr) are somewhat limited in range and their color scaling is logarithmic, except for panels f and i where it is the square root of intensity, and panel g, where it is linear.

Fig. 7 Application of getsources and getfilaments to the flattened detection images ℐλD shown in Fig. 6i. Different structural components are separated as independent images of filaments (a) and sources (b). For reference, the extraction ellipses for all detected sources (most of them may be background galaxies) are also overlaid on the original observed image ℐλ (c). Intensities (in MJy/sr) are plotted with logarithmic color scaling.

The quality of the flattening algorithm is visualized in Fig. 5 by the image of sd₉ (ℐ_λD) computed from the flat detection images. A comparison with ${\mathrm{𝒟}}_{\mathit{\lambda }}\hspace{0.17em}\mathrm{=}\hspace{0.17em}{\mathrm{sd}}_{\mathrm{9}\hspace{0.17em}}\mathrm{\left(}\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathrm{𝒮}}_{\mathit{\lambda }}\end{array}\mathrm{\right)}$ makes it quite clear that ℐ_λD have remarkably uniform fluctuations (outside the sources). Such flat images are optimal for detecting sources (with H_λ ≲ X_λ) using extraction methods that employ global thresholding.

3. Discussion

The background estimation and flattening method getimages was validated in Sect. 2 using a simulated star-forming region. To demonstrate its performance on very complex observed images, it was applied to a Spitzer image³ of the star-forming region L 1688 in Ophiuchus⁴. The only aim of this section is to clarify various aspects of the new method and any astrophysical analysis based on the image is beyond the scope of this paper.

3.1. Application to L 1688 in Ophiuchus

The Spitzer 8 μm image of L 1688 (1 × 1°with a 2′′resolution, Fig. 6a) displays very complex intensity distribution ℐ_λ, with the background varying on different spatial scales by more than two orders of magnitude. In addition to many compact sources in both faint and bright background areas, it also shows some filamentary structures, the most prominent one running up from the middle of the image. Sources blended with such complex and variable background are very difficult to extract reliably and measure accurately using automated methods. To apply the getimages method, it is necessary to specify a maximum size X_λ for the structures of interest (cf. Sect. 3.2). In this application to L 1688, X_λ = 12″ was adopted and as a variation, a higher value of 24′′was also used in order to demonstrate the effects of this parameter. The results obtained after M = 30 iterations are presented below.

The background $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathrm{ℬ}}_{\mathit{\lambda }}\end{array}$ of the L 1688 image, derived by getimages, is displayed in Fig. 6b. The median-filtering algorithm described by Eqs. (2)–(5) has cleanly erased all sources of various sizes, perhaps with the exception of two large and very bright sources below the image center, which appear to be incompletely removed. The two residual maxima in $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathrm{ℬ}}_{\mathit{\lambda }}\end{array}$ may indeed be parts of the extended power-law sources, or they might also be unrelated background structures. Similarly, the prominent vertical filament seems to have left some residuals in $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathrm{ℬ}}_{\mathit{\lambda }}\end{array}$. Nevertheless, the main part of the structures was adequately truncated and small residuals at that level are unimportant for detection images.

The background-subtracted image $\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathrm{𝒮}}_{\mathit{\lambda }}\end{array}\hspace{0.17em}\mathrm{=}\hspace{0.17em}{\mathrm{ℐ}}_{\mathit{\lambda }}\hspace{0.17em}\mathrm{-}\hspace{0.17em}\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathrm{ℬ}}_{\mathit{\lambda }}\end{array}$ of L 1688 is presented in Fig. 6c. Without the dominating bright and strongly variable background, the image becomes substantially simpler, displaying all sources and structures with widths X_λ ≲ 12′′much more clearly. Obviously, the brightness remains highly fluctuating across $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathrm{𝒮}}_{\mathit{\lambda }}\end{array}$, and subtraction of a large-scale background cannot improve the image to ensure a problem-free source extraction.

The fluctuations are quantified by the image of standard deviations ${\mathrm{𝒟}}_{\mathit{\lambda }}\hspace{0.17em}\mathrm{=}\hspace{0.17em}{\mathrm{sd}}_{\mathrm{9}\hspace{0.17em}}\mathrm{\left(}\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathrm{𝒮}}_{\mathit{\lambda }}\end{array}\mathrm{\right)}$ in L 1688 shown in Fig. 6d. Pixel values spanning several orders of magnitude highlight the difficulties that are encountered when using thresholds in terms of the global standard deviations to separate structures from background and noise. For such images, a solution would be either to measure local fluctuations around each structure or to equalize background and noise fluctuations across the entire image $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathrm{𝒮}}_{\mathit{\lambda }}\end{array}$. The former approach is very problematic, as it would require an accurate prior source extraction, and the latter is the flattening approach adopted by getimages.

To render the fluctuations uniform, the method computes a median-filtered background $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathrm{ℱ}}_{\mathit{\lambda }}\end{array}$ of the standard deviations , according to the algorithm of Eqs. (6) and (7). The scaling (flattening) image $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathrm{ℱ}}_{\mathit{\lambda }}\end{array}$ shown in Fig. 6e represents background and noise fluctuations in L 1688, excluding all structures of widths X_λ ≲ 12′′. Similarly to $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathrm{ℬ}}_{\mathit{\lambda }}\end{array}$, which approximates the local background for all sources and filaments, $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathrm{ℱ}}_{\mathit{\lambda }}\end{array}$ evaluates their local fluctuation levels. To equalize the latter over the entire image, it is sufficient to divide $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathrm{𝒮}}_{\mathit{\lambda }}\end{array}$ by the scaling image $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathrm{ℱ}}_{\mathit{\lambda }}\end{array}$.

The flattened detection image ${\mathrm{ℐ}}_{\mathit{\lambda }\mathrm{D}}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathrm{𝒮}}_{\mathit{\lambda }}\end{array}\mathit{/}\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathrm{ℱ}}_{\mathit{\lambda }}\end{array}$ for L 1688 is presented in Fig. 6f. Compared with the background-subtracted image $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathrm{𝒮}}_{\mathit{\lambda }}\end{array}$ (Fig. 6c), the scaled image of L 1688 appears remarkably flat. This visual impression is further quantified by the image of standard deviations sd₉ (ℐ_λD) in Fig. 6g, and it is confirmed by comparison with the image computed before the flattening (Fig. 6d). Such greatly simplified (background-subtracted and flat) detection images are highly beneficial for extracting sources and filamentary structures.

Possible imperfections of the results produced by getimages are related to an inadequate choice of X_λ and they must become clearly visible in $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathrm{ℬ}}_{\mathit{\lambda }}\end{array}$ and $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathrm{ℱ}}_{\mathit{\lambda }}\end{array}$. Indeed, Figs. 6b and e show some residuals in places of the vertical filament and of the two bright sources, above and below the center of the images, respectively. Such levels of inaccuracies in ℐ_λD are not important for source extraction with getsources (Paper I) because the source measurements are made in the original image ℐ_λ. One might be interested, however, in a more complete reconstruction and accurate profiling of the prominent filament seen in Fig. 6. Filament extraction with getfilaments (Paper II) is made in detection images ℐ_λD, using the scaling images as well to obtain correct intensities. To minimize the residuals in $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathrm{ℬ}}_{\mathit{\lambda }}\end{array}$ and $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathrm{ℱ}}_{\mathit{\lambda }}\end{array}$, it is sufficient to simply rerun getimages with an increased value of X_λ.

An illustration of the effect of a twice higher value X_λ = 24′′on the resulting images for L 1688 is presented in Figs. 6h and i. An inspection of the images shows that the larger maximum size leads to lower residuals in $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathrm{ℱ}}_{\mathit{\lambda }}\end{array}$. The corresponding flattened image ℐ_λD displays structures that are brighter and roughly twice as wide (at zero level) than those in Fig. 6f derived with X_λ = 12′′. With the higher parameter value, the vertical filament is extracted more completely, and its profile can be determined more accurately.

Figure 7 shows an example of the source and filament extraction in the L 1688 field using a flattened detection image produced by getimages (Fig. 6i). The getfilaments method extracts filamentary structures in the image and subtracts them before extracting compact sources with getsources. The three methods carefully separate blended structural components of the original image ℐ_λ (background, filaments, and sources) from each other. The linear observational artifacts seen as orthogonal spikes next to bright sources are automatically identified as “filaments” (Fig. 7a) and therefore removed before the source extraction. Owing to the radically simplified detection image, most of the extracted 1193 sources (Fig. 7b) in the L 1688 field are real. A small number of the remaining spurious sources are located around the unresolved bright peaks that are due to the artifacts in the observed image, which are induced by the complex shape of the observational beam.

3.2. Practical definition of the maximum size X_λ

The maximum structure size is the only free parameter of getimages. Before applying the method to an image ℐ_λ, it is necessary to evaluate X_λ for the sources and filaments of interest in that image. Since human eyes are good at identifying sources even on very complex backgrounds, the maximum size can easily be estimated directly from the observed ℐ_λ or, alternatively, assumed on the basis of additional considerations. The parameter X_λ controls the quality of the derived background, therefore it must be related to the size of a zero-level footprint, not to the half-maximum size H_λ of the structures of interest.

Following Papers I and II, footprints are defined here as areas of connected pixels that make a non-negligible contribution to the total fluxes of sources or filaments. The footprint size Z_λ is measured in ℐ_λ at a radial distance from the peak or crest, where the structure intensity completely merges with background and noise fluctuations (“zero level”). Therefore, Z_λ has the meaning of the major axis of an elliptical source or the largest width of a filamentary structure at zero intensity level after background subtraction. The relationship between H_λ and Z_λ varies for different structures, depending on their intensity distribution, signal-to-noise ratio S/N, and background properties. To define the half-maximum X_λ as a proxy to zero-level Z_λ, it is necessary to eliminate the dependence on unknown intensity profiles.

It makes sense to convert Z_λ into X_λ assuming an equivalent Gaussian intensity profile and determine the getimages parameter as X_λ = η^-1Z_λ, where η is a scaling factor within the range of 2−3. To create ℐ_λD exclusively for detection, it is sufficient to adopt η = 3 (see Fig. 8). To subtract the background more accurately and open the possibility of using ℐ_λD for measurements as well, it is advisable to adopt η = 2, which translates Z_λ into larger X_λ.

Incorporating angular resolution O_λ as a lower limit, getimages defines the parameter as ${\mathit{X}}_{\mathit{\lambda }}\mathrm{=}\begin{array}{c}\max \end{array}\mathrm{\{}\hspace{0.17em}{\mathit{\eta }}^{-1}{\mathit{Z}}_{\mathit{\lambda }}\mathit{,}\hspace{0.17em}{\mathit{O}}_{\mathit{\lambda }}{\hspace{0.17em}}^{\mathrm{\}}}\mathit{.}$(8)This formulation also gives reasonable values of X_λ for faint unresolved sources whose Z_λ values obtained from the observed images are likely to be underestimated. An obvious upper limit for X_λ is a reasonably small fraction (≲ 5%) of the image size to prevent the largest windows from sliding beyond the image edges. In practice, it is recommended to restrict the parameter X_λ to the lowest values possible, as median filtering with unnecessarily large windows affects much greater areas of images and hence might make the results less local.

The above definition is illustrated in Fig. 8, which shows model examples of Gaussian and power-law sources with an FWHM size H = 8′′and peak intensity I_P = 100 (in arbitrary units) on a simple planar background with additional random noise. The background brightness (1 and 10) and noise fluctuations (σ = 1 and 10) were chosen to simulate both bright and faint sources. The bright sources have S/N = 100 (lower curves), whereas the relatively faint sources correspond to S/N = 10 (upper curves).

An estimation of the footprint size would give values Z ≈ 24 and 70′′for the bright Gaussian and power-law sources, respectively (Fig. 8). According to the definition in Eq. (8), the values correspond to X = 8 and 23′′(for η = 3), and the maximum size for the steep Gaussian intensity distribution equals the source size H. For much fainter sources, the low-intensity outskirts of the source intensity distributions dissolve into much stronger background and noise fluctuations (Fig. 8). An evaluation of the zero-level size would therefore underestimate the true source values, giving Z ≈ 18 and 26′′, corresponding to X = 8 and 8.7′′. For the faint sources, Z is underestimated and X = H for the Gaussian source because of the presence of O_λ in Eq. (8).

Fig. 8 Definition of the maximum size X. The two sets of curves represent faint (upper) and bright (lower) sources, with an S/N of 10 and 100. Shown are profiles of Gaussian (blue) and power-law (red) sources with a peak intensity IP = 100 and an FWHM size H = 8′′ on a flat background (green), affected by random Gaussian noise (with σ = 10 and 1). The noise-free profiles are shown by black curves. Dashed horizontal lines (blue and red) visualize zero-level sizes Z that would be estimated for the two sources and the solid horizontal lines (cyan and magenta) show the corresponding X values (assuming η = 3). In the high-S/N case, Z ≈ 24 and 70′′, whereas X = 8 and 23′′. In the low-S/N case, Z ≈ 18 and 26′′, whereas X = 8 and 8.7′′.

From the definition of the maximum size, X_λ>H_λ for power-law sources and hence R_λN = 4 X_λ> 4 H_λ (Sects. 2.2 and 2.3). The model values of truncation factors f_T presented in Table 1 (Sect. 2.1) as a function of R/H are therefore only lower limits for power-law sources. The actual truncation factors and hence the quality of source removal by median filtering in getimages are much better. For example, images of the power-law sources profiled in Fig. 8 show that f_T for R/X = 4 must increase by a factor of 5.6 in comparison with the R/H = 4 values listed in Table 1.

4. Conclusions

This paper described getimages, a new general method of background estimation and image flattening, which solves the two problems in a multiscale median-filtering approach. The method does not need any preliminary source (or filament) extractionor any prior information about the location and sizes of sources. A single free parameter of getimages, the maximum size X_λ of the structures of interest, can easily be evaluated directly from the observed images ℐ_λ (Sect. 3.2). The method produces background-subtracted and flattened images that are radically simplified by the removal of huge intensity variations of large-scale backgrounds (Sect. 2.2) and global equalization of the small-scale fluctuation levels (Sect. 2.3), for all structures with sizes H_λ ≲ X_λ. The resulting flat images ℐ_λD are highly advantageous for detecting sources and other structures, whereas measurements of the physical quantities are usually carried out in the original observed images ℐ_λ. When accurately derived, the background-subtracted images $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathrm{𝒮}}_{\mathit{\lambda }}\end{array}$ can also be used for measurements.

Any source or filament extraction method that clearly distinguishes between the detection and measurement images could greatly benefit from using the flat detection images produced by getimages. Background subtraction and flattening are especially important for the methods that employ global thresholding to detect sources. Specifically, the getsources and getfilaments source and filament extraction methods (Papers I and II) are greatly simplified and improved when using the flat detection images ℐ_λD. Instead of the two (initial and final) extractions of the original approach, a single extraction is now sufficient (see Appendix C for details). Most importantly, the new flattening algorithm of getimages is much more reliable and universal than the original algorithm.

Observational and map-making artifacts are also eliminated or greatly reduced in the flattened detection images. For example, mosaicking of several independently observed images of the same field often produces significantly deviating noise levels in the tiles of a resulting large image. As a bonus, getimages automatically equalizes the levels of noise fluctuations between the different tiles, which makes the composite detection image seamless.

Acknowledgments

This study employed SAOImage DS9 (by William Joye) developed at the Smithsonian Astrophysical Observatory (USA), the CFITSIO library (by William D. Pence) developed at HEASARC NASA (USA), and SWarp (by Emmanuel Bertin) developed at Institut d’Astrophysique de Paris (France). The plot utility and ps12d library used in this work to draw figures directly in the PostScript language were written by the author using the PSPLOT library (by Kevin E. Kohler) developed at Nova Southeastern University Oceanographic Center (USA) and the plotting subroutines from the AZEuS MHD code (by David A. Clarke and the author) developed at Saint Mary’s University (Canada). This work used observations made with the Spitzer Space Telescope, which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under a contract with NASA. The author appreciates collaborative work within the Herschel Gould Belt and HOBYS surveys key projects. Useful comments on a draft made by Pierre Didelon helped improve this paper.

References

Appendix A: Model intensity distributions

Fig. A.1 Intensity profiles of radiative transfer models. The three panels from left to right display results for masses of 0.3, 3, and 30 M⊙. In each panel, the lower and upper sets of curves show background-subtracted profiles of starless cores and protostellar envelopes, respectively, for selected Herschel wavelengths of 160, 250, 350, and 500 μm. For reference, the black dashed curves visualize a Gaussian profile and a power-law slope Iλ ∝ θ2, whereas the dotted vertical line indicates the outer radius of the models.

This section employs radiative transfer modeling to justify the choice of simple intensity shapes in Sect. 2.1. It presents intensity profiles in selected Herschel bands for several radiative transfer models of spherical starless cores and protostellar envelopes with masses 0.3, 3, and 30 M_⊙, calculated and described in detail by Men’shchikov (2016). The starless cores have a flat-topped density structure of an isothermal Bonnor-Ebert sphere (Bonnor 1956). The protostellar envelopes have a power-law radial density distribution ρ(r) ∝ r^-2 and accretion luminosity with the numerical value (in L_⊙) of the model mass. The outer boundary of the models was placed at a radial distance of 10⁴ AU, beyond which an embedding constant-density cloud was assumed. At the adopted distance of 140 pc, the outer model radius corresponds to 71′′. For more details of the radiative transfer models, see Sect. 2 in Men’shchikov (2016).

Figure A.1 displays intensity profiles for both starless cores and protostellar envelopes after subtracting the background emission produced by the embedding cloud. For simplicity, obvious angular resolution effects are ignored, therefore the curves show the true distribution of the model intensities (with a pixel size of 0.67′′). The intensity profiles of the starless cores resemble Gaussians and those of the protostellar envelopes can be approximated by a power law I_λ ∝ θ^-2. Although the profiles of radiative transfer models are wavelength-dependent and more complex than the simple Gaussian and power-law shapes, the latter are reasonable choices to illustrate the source removal by median filtering in Sect. 2.1.

Appendix B: Images of simulated structures

Fig. B.1 Simulated sources and filaments with H = 8′′ and their median-filtered W = 4 counterparts (cf. Sect. 2.1). The two upper rows show the simulated Gaussian structures , , in panels a–d and their median-filtered images in panels e–h. The two bottom rows show the simulated power-law structures , , , in panels i–l and their median-filtered images in panels m–p. The images with an intensity range of [−0.3,100 ] and sizes of 200′′ are profiled in Fig. 1. The images are displayed with logarithmic color scaling.

The structure removal by median filtering described in Sect. 2.1 is further illustrated in Fig. B.1 by images of the simulated sources and filaments with sizes of H = 8′′and median-filtered images obtained with W = 4. The intensity profiles shown in Fig. 1 and truncation factors listed in Table 1 were measured along horizontal lines through the image centers.

Appendix C: Implications for source and filament extraction with getsources and getfilaments

This section describes an updated approach to using getimages in combination with the extraction method presented in Papers I and II. The old strategy involved an image preparation step and two full extractions separated by a flattening step (cf. Fig. 20 in Paper I and Fig. 1 in Paper II). The new approach, outlined in Fig. C.1, requires only a single run of getsources (getfilaments), preceded by the runs of the image preparation script prepareobs (part of getsources) and of the new getimages method. The latter was coded in a Bash script⁵ as a series of calls to the FORTRAN utilities of getsources (cf. Appendix G in Paper I), and it uses a configuration file that is a subset of the getsources configuration file.

With the new approach and a single extraction run, default values of several parameters in the getsources configuration file need to be modified. The maximum sizes X_λ used in the getimages run replace the old default value of 220. For example, if X_λ equals twice the Herschel beam size O_λ, the values are used as the third parameter on the following lines:

                 070  8.4 16.8 0 1 y 0 | ................. 
                 100  9.4 18.8 0 1 y 0 | ................. 
                 160 13.5 27.0 0 1 y 0 | ................. 
                 250 18.2 36.4 0 1 y 0 | ................. 
                 350 24.9 49.8 0 1 y 0 | ................. 
                 500 36.3 72.6 0 1 y 0 | .................
        

When the flattened detection images are produced in a prior run of getimages, the parameter flattening must be set to n:

n <--user.. | flattening dofl ...

The parameters sreliable, stentative, and contranoise must have a plus in front of them, and the default values

of nsigmacutss, nsigmacutfloor must be set to 6:

          +7          ..user..     |   sreliable .........
          +5          ..user..     |   stentative ........ 
             6          ..expert   |   nsigmacutss ....... 
             6          ..expert   |   nsigmacutfloor ....
          +1.3       ..expert   |   contranoise .......

These are the only differences with respect to the default configuration parameters of getsources used in the original strategy described in Paper I.

Fig. C.1 Flowchart of the extraction approach with getsources and getfilaments when using the flattened detection images produced by getimages. The three blocks represent the prepareobs script (blue), the getimages method (green), and the getsources (getfilaments) method (red).

All Tables

All Figures

Fig. 1 Median filtering applied to Gaussian sources (upper row), power-law sources (second row), Gaussian filaments ℱG (third row), and power-law filaments ℱP (bottom row) with sizes H (FWHM) of 8′′ and 16′′ (left) and to the same images after addition of random Gaussian noise with σ = 0.33 (middle) and background ℬ modeled as a large Gaussian with FWHM of 512′′ (right). The radii R of circular sliding windows given by the indices on the curve labels correspond to the relative radii W = R/H = { 1,2,4,8 }. Original intensity profiles of , , ℱG, and ℱP are shown in gray, filtered profiles of the sources and filaments with H = 8 and 16′′ are plotted in red and blue, whereas profiles of and ℬ are colored in cyan and green, respectively. Large sliding windows (W> 2) truncate Gaussian shapes so efficiently (cf. Table 1) that some of the annotated curves become invisible (fall entirely below the lower edge of the plots).

In the text

Fig. 2 Background derivation by median filtering. Shown are the images ℐλ of a simulated star-forming region (upper row), median-filtered background $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathrm{ℬ}}_{\mathit{\lambda }}\end{array}$ (middle row), and background-subtracted images $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathrm{𝒮}}_{\mathit{\lambda }}\end{array}$ (bottom row) at selected Herschel wavelengths (after M = 30 iterations). The maximum source sizes Xλ of 25, 100, and 100′′ for the getimages method were estimated directly from ℐλ (cf. Sect. 3.2). In panel a, small holes are starless cores seen in absorption at 70 μm, and white emission peaks are protostellar sources. Intensities (in MJy/sr) are limited in range, and their color scaling is linear.

In the text

Fig. 3 Quality of the derived background presented in Fig. 2. Shown are the relative accuracies $\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathrm{ℬ}}_{\mathit{\lambda }}\end{array}\mathit{/}{\mathrm{ℬ}}_{\mathit{\lambda }}\hspace{0.17em}\mathrm{-}\hspace{0.17em}\mathrm{1}$ of the background $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathrm{ℬ}}_{\mathit{\lambda }}\end{array}$ obtained by median filtering with respect to the true background ℬλ. Only the relevant pixels within the model sources are shown for values limited by the range [−0.3,0.3 ] with linear color coding. In panel a, the minima visible in the centers of some footprints are caused by the absorption of radiation in the central parts of starless cores; the two deepest minima are at levels of −0.7 and −0.6. In panels b and c, the maxima are at levels of 0.8 and 0.7, respectively.

In the text

Fig. 4 Flattening background-subtracted images (Fig. 2). Shown are the standard deviations ${\mathrm{𝒟}}_{\mathit{\lambda }}\hspace{0.17em}\mathrm{=}\hspace{0.17em}{\mathrm{sd}}_{\mathrm{9}\hspace{0.17em}}\mathrm{\left(}\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathrm{𝒮}}_{\mathit{\lambda }}\end{array}\mathrm{\right)}$ of small-scale fluctuations in $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathrm{𝒮}}_{\mathit{\lambda }}\end{array}$ (upper row), median-filtered scaling (flattening) images $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathrm{ℱ}}_{\mathit{\lambda }}\end{array}$ (middle row), and flattened detection images ${\mathrm{ℐ}}_{\mathit{\lambda }\mathrm{D}}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathrm{𝒮}}_{\mathit{\lambda }}\end{array}\mathit{/}\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathrm{ℱ}}_{\mathit{\lambda }}\end{array}$ (bottom row) at selected Herschel wavelengths. In panel b, thin ring-like structures reflect the off-center temperature peaks in starless cores illuminated by the interstellar radiation field. The maximum sizes Xλ of 25, 100, and 100′′ are the same as those used for background derivation.

In the text

Fig. 5 Quality of the flattened images presented in Fig. 4. Shown are the standard deviations sd9 (ℐλD) of small-scale fluctuations in the flattened images ℐλD. The fluctuations outside sources are very uniform, compared to ${\mathrm{𝒟}}_{\mathit{\lambda }}\hspace{0.17em}\mathrm{=}\hspace{0.17em}{\mathrm{sd}}_{\mathrm{9}\hspace{0.17em}}\mathrm{\left(}\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathrm{𝒮}}_{\mathit{\lambda }}\end{array}\mathrm{\right)}$ computed from the background-subtracted images.

In the text

Fig. 6 Application of getimages (with Xλ = 12′′ and M = 30) to real-life (Spitzer λ = 8 μm) observations of the Ophiuchus L 1688 star-forming region: a) original image ℐλ, b) estimated background $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathrm{ℬ}}_{\mathit{\lambda }}\end{array}$, c) background-subtracted image $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathrm{𝒮}}_{\mathit{\lambda }}\end{array}$, d) standard deviations , e) scaling image $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathrm{ℱ}}_{\mathit{\lambda }}\end{array}$, f) flattened detection image ℐλD, and g) standard deviations sd9 (ℐλD). Selected results for Xλ = 24′′ are also shown: h) scaling image $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathrm{ℱ}}_{\mathit{\lambda }}\end{array}$ and i) flattened image ℐλD. Intensities (in MJy/sr) are somewhat limited in range and their color scaling is logarithmic, except for panels f and i where it is the square root of intensity, and panel g, where it is linear.

In the text

Fig. 7 Application of getsources and getfilaments to the flattened detection images ℐλD shown in Fig. 6i. Different structural components are separated as independent images of filaments (a) and sources (b). For reference, the extraction ellipses for all detected sources (most of them may be background galaxies) are also overlaid on the original observed image ℐλ (c). Intensities (in MJy/sr) are plotted with logarithmic color scaling.

In the text

Fig. 8 Definition of the maximum size X. The two sets of curves represent faint (upper) and bright (lower) sources, with an S/N of 10 and 100. Shown are profiles of Gaussian (blue) and power-law (red) sources with a peak intensity IP = 100 and an FWHM size H = 8′′ on a flat background (green), affected by random Gaussian noise (with σ = 10 and 1). The noise-free profiles are shown by black curves. Dashed horizontal lines (blue and red) visualize zero-level sizes Z that would be estimated for the two sources and the solid horizontal lines (cyan and magenta) show the corresponding X values (assuming η = 3). In the high-S/N case, Z ≈ 24 and 70′′, whereas X = 8 and 23′′. In the low-S/N case, Z ≈ 18 and 26′′, whereas X = 8 and 8.7′′.

In the text

Fig. A.1 Intensity profiles of radiative transfer models. The three panels from left to right display results for masses of 0.3, 3, and 30 M⊙. In each panel, the lower and upper sets of curves show background-subtracted profiles of starless cores and protostellar envelopes, respectively, for selected Herschel wavelengths of 160, 250, 350, and 500 μm. For reference, the black dashed curves visualize a Gaussian profile and a power-law slope Iλ ∝ θ2, whereas the dotted vertical line indicates the outer radius of the models.

In the text

Fig. B.1 Simulated sources and filaments with H = 8′′ and their median-filtered W = 4 counterparts (cf. Sect. 2.1). The two upper rows show the simulated Gaussian structures , , in panels a–d and their median-filtered images in panels e–h. The two bottom rows show the simulated power-law structures , , , in panels i–l and their median-filtered images in panels m–p. The images with an intensity range of [−0.3,100 ] and sizes of 200′′ are profiled in Fig. 1. The images are displayed with logarithmic color scaling.

In the text

	Fig. C.1 Flowchart of the extraction approach with getsources and getfilaments when using the flattened detection images produced by getimages. The three blocks represent the prepareobs script (blue), the getimages method (green), and the getsources (getfilaments) method (red).
In the text