© J.-S. Carrière et al. 2022

1 Introduction

A large variety of methods have already been developed to extract elongated structures from two-dimensional (2D) maps. These approaches may be divided into three main categories: some focus on a purely local analysis of the structures based on local derivatives at each pixel; others adopt a non-local approach, exploring a larger space around each pixel to look for specific scales; the third category of methods propose a global analysis of the whole map, applying a multi-scale decomposition.

Local methods compute either the gradient (first-order derivatives; e.g., Soler et al. 2013; Planck Collaboration XXXV 2016b) or the Hessian matrix (second-order derivatives; e.g., Polychroni et al. 2013; Schisano et al. 2014; Planck Collaboration XXXII 2016a) at each pixel of the considered map. In some cases, the main purpose is to derive the orientations of elongated structures for statistical purposes (e.g., Soler et al. 2013; Planck Collaboration XXXII 2016a; Planck Collaboration XXXV 2016b). In other cases, filament skeletons are extracted by connecting contiguous pixels along the crests of the (intensity or column-density) distribution. For instance, the DisPerSe method, originally developed to recover filament skeletons in cosmic web maps (Sousbie 2011), was successfully applied to Herschel column density maps (Arzoumanian et al. 2011, 2019; Peretto et al. 2012; Palmeirim et al. 2013) and to ¹³CO intensity maps (Panopoulou et al. 2014). A limitation of the local approach is the difficulty of detecting faint structures such as striations¹.

Non-local methods focus on a given scale around each pixel. The template-matching approach (Juvela 2016) allows the user to look for any specific morphology of given dimensions and to build the probability of finding such oriented structures through a kernel convolution on the map. Another efficient approach is the Rolling Hough Transform (RHT, Clark et al. 2014) method, which computes an estimator of the level of linearity of the structures in the neighborhood of a pixel at a given scale, making use of the Hough transform. This method was extensively used in recent studies, in HI, Herschel, and Planck maps (Clark et al. 2015; Clark & Hensley 2019; Malinen et al. 2016; Panopoulou et al. 2016; Alina et al. 2019). A third approach is the filfinder method, which extracts filament skeletons (Koch & Rosolowsky 2015). However, in contrast to DisPerSe, filfinder starts with a spatial filtering at a given scale and covers a dynamic range broad enough to include striations.

Finally, global methods offer a multi-scale and complete analysis of a field. The getfilaments method (Men’shchikov 2013) extracts a filament network with the help of statistical tools and morphological filtering to remove background noise. It also extracts point sources and is able to perform a multi-waveband analysis. This method is very complete, but it requires some fine-tuning and additional tools to extract the filament orientations and scales. It was already applied to Herschel maps (Cox et al. 2016; Rivera-Ingraham et al. 2016, 2017; Arzoumanian et al. 2019). The wavelet-based methods by Robitaille et al. (2019) and by Ossenkopf-Okada & Stepanov (2019) use an anisotropic-wavelet analysis to extract a whole filament network by analyzing fluctuations in the map as functions of spatial scale. This kind of method may circumvent the well-known biases of other typically used methods (Panopoulou et al. 2017) and remains relatively fast. However, wavelet-based methods require additional steps to extract filament orientations. Although these methods are multi-scale, the wavelet analysis implies a logspace scaling, which results in a lower resolution at higher scale.

A full comparison between these different methods would be extremely useful, but to date only a few limited comparisons exist. For instance, Juvela (2016) compared the template-matching and RHT methods applied to simulation data. He found that both methods give equally good results, except in simulations with significant noise and background fluctuations, for which template matching performs better. More recently, Micelotta et al. (2021) compared the gradient and RHT methods, again on simulation data. These authors found both similarities and disparities in the results, and they attributed the disparities to intrinsic differences in the filamentary structures selected by both methods.

Here, we would like to study the relative orientations between filaments and the local magnetic field. To that end, we need a filament-extraction method that can operate in a broad range of Galactic environments, from dense and complex structures to the more diffuse (neutral) medium, and can provide robust and homogeneous filament detection in various fields for a multi-scale statistical analysis. While none of the methods described above are fully satisfactory for our purpose, the closest is probably RHT, which is easy to use and has already given promising results (see, e.g., Clark et al. 2014; Malinen et al. 2016; Panopoulou et al. 2016; Alina et al. 2019). Therefore, we decided to start from RHT and adapt it to match our requirements. This led us to develop a new method, called FilDReaMS (Filament Detection and Reconstruction at Multiple Scales).

In Sect. 2, we present the main features of our new FilDReaMS method, with reference to RHT. In Sect. 3, we provide a detailed description of the FilDReaMS methodology. In Sect. 4, we present the results of several series of simulations designed to validate FilDReaMS. In Sect. 5, we apply FilDReaMS to the Herschel G210 field and compare our results to those previously obtained with RHT. In Sect. 6, we conclude our paper.

2 General overview

2.1 Introduction to FilDReaMS

Let us consider a map of a given quantity I, which, in the astrophysical context, can represent intensity, column density, or temperature. We use this map as an example to present the FilDReaMS method. For simplicity, we refer to I as an intensity, keeping in mind that I actually has a broader meaning.

The purpose of applying FilDReaMS to the map of I is to identify filamentary structures over a range of spatial scales and intensities, from the largest and brightest filaments down to striations. By considering that filaments can be locally approximated by rectangular bars, FilDReaMS is able to extract two of their characteristics: the widths and the orientations of the associated bars.

In the following, the rectangular bar used by FilDReaMS is referred to as the model bar. It is characterized by a width W_b, a length L_b, and an aspect ratio r_b = L_b/W_b. For any filament detected with a model bar of width W_b, W_b is referred to as the bar width of the filament.

The orientation angle of the model bar, ψ_b, is defined with respect to a given north direction (e.g., Galactic north for an astrophysical map) and taken to increase counterclockwise from north. This definition is consistent with the IAU convention for polarization angles. As the model bar is symmetric, ψ_b can be defined over a 180° range, which we choose to be [−90°, +90°].

We adopt the same convention for the orientation angle of a filament, ψ_f, defined in Sect. 3.3. Furthermore, when considering the difference between two angles defined in [−90°, +90°], we require the result to also lie in the range [−90^°, +90^°], possibly by adding or subtracting 180^°. All the parameters related to FilDReaMS are listed in Table 1.

To provide a first application of FilDReaMS, we selected one of the Herschel fields observed in the Galactic cold core (GCC) key-program (Juvela et al. 2010, 2012), the so-called G210 field, which corresponds to the high-Galactic-latitude star-forming region L1642. This region was studied in detail by Malinen et al. (2016), who investigated the relative orientations between the magnetic field (traced with Planck polarization data) and filaments extracted from the G210 Herschel map using the RHT method. G210 is located at Galactic longitude l = 210.90°, Galactic latitude b = −36.55°, and distance d = 140 ± 20pc (Montillaud et al. 2015). Its H₂ column density map has dimensions Δl × Δb = 1.28° × 1.22°, corresponding to 3.1 pc × 3.0pc. The angular size of a pixel, equal to one-third of the beam size, is 12″, corresponding to 0.0081 pc. The H₂ column density, , varies over the range [0.2,7.5] × 10²¹ cm⁻².

2.2 Overview of RHT

The RHT method developed by Clark et al. (2014) is based on the Hough transform designed to search for straight lines in images, even when partially filled or dominated by noise. Defining a straight line with two parameters (the orientation θ_Η of its normal and the minimal distance ρ_Η to the origin), the Hough transform allows us to pass from the pixel space (x, y) to this other parameter space (ρ_Η, θ_Η), by quantifying the number of pixels located on the same linear feature. In the case of the RHT method, in order to focus on straight features and to suppress large-scale structures, a binary version of the original image is first built by subtracting the smoothed image with a top-hat and thresholding to zero. The Hough transform is then applied locally on this binary image: for each pixel, a simplified version of the Hough transform can be applied (with ρ_Η = 0) inside a circular area defined around this central pixel in order to pass from pixel space to (θ_Η) space, and to build the distribution of orientations of the linear features around this central pixel. Once applied iteratively on every pixel of the whole image, the last step of the method consists in choosing a common threshold above which the distributions of orientations are stored and used to derive local average orientation or total RHT intensity over the whole image.

This method is extremely powerful to get an estimate of the linearity level inside local regions of images irrespective of the overall brightness of the region. It still suffers from a few limitations. Firstly, the use of the Hough transform, which performs transformation from pixel space to (ρ_Η, θ_Η) space, directly implies that RHT distributions of orientation are dependent of the pixelization of the image, which could lead to subtle bias of the results. Secondly, the choice of the threshold used to cut the distributions of orientations is arbitrary and may impact the analyses from one image to another.

2.3 RHT versus FilDReaMS

FilDReaMS tries to overcome some of the limitations of RHT, as explained in the rest of this section. It also makes it possible to access additional information about the widths (more exactly, the bar widths) of the detected filaments.

The sensitivity of the algorithm to pixelisation is inherent in the basic implementation of the RHT, which uses relative positions of pixel centers to a given pixel (centered on x, y) to build a θ_Η representation. FilDReaMS solves this problem by a fundamental change of philosophy: it starts from a discretization of the θ_Η space to compute the intersection of a rotated bar with any pixel area centered on (x, y).

In FilDReaMS, the arbitrary choice of threshold in RHT for determining linearity significance is alleviated by a comparison to a random distribution obtained locally using ideal template bars. This process takes into account the noise level and the complexity of the region, and provides us a robust assessment of the reliability of the detection.

The determination of a preferred orientation based on the orientation distribution in each pixel is improved in FilDReaMS by performing a search for local maxima, allowing us to derive multiple peaks in the orientation distribution with individual significances, instead of averaging the orientation angle estimated over the whole orientation distribution.

2.4 The main steps of FilDReaMS

Briefly, the main steps of FilDReaMS are as follows: (1) Spatial filtering: a 2D top-hat filtering is applied to the input image, which is then converted to a binary map that contains only structures narrower than a given width; see Sect. 3.1. (2) Building an orientation distribution: the matching between the binary map and a model bar with given width W_b and variable orientation is evaluated at each pixel i in order to build an orientation distribution; see Sect. 3.2. (3) Detection of preferred orientations: local maxima in the orientation distribution at each pixel i are identified with preferred orientations; see Sect. 3.3. (4) Reliability assessment: the significance of each preferred orientation is assessed by comparison with an ideal orientation distribution; see Sect. 3.4. (5) Reconstruction of physical filaments: the true shape and the intensity of physical filaments is reconstructed from the initial image masked by the binary map, and a filament orientation angle at each pixel ί′, , is derived; see Sect. 3.5. (6) Iteration over various bar widths: the procedure is repeated for a range of values of W_b in order to derive the most significant bar width at each pixel i′, , as well as the most prevalent bar widths for the entire map, ; see Sect. 3.6.

While step (1) is fully similar to the RHT processing, steps (2) and (3) have the same objectives as RHT but a totally different implementation, allowing us to address the pixelization sensitivity (see Sect. 3.2) and the multiplicity of the local preferred orientations in the case where linear structures are crossed (see Sect. 3.3). Finally, steps (4) to (6) are entirely new and specific to FilDReaMS.

Fig. 1 Illustration of the FilDReaMS method applied to one pixel i of the H2 column density map of the Herschel G210 field. Top left: initial map A of the entire G210 field. Top middle: smoothed map B obtained by convolving map A with a 2D top-hat kernel. Top right: corresponding binary map C in which the yellow pixels have a value of 1 and the dark pixels a value of 0. Bottom right: subregion Ci of map C centered on pixel i, with the model bar (width Wb, length Lb, and orientation angle ψb) overplotted in green. Bottom middle: measured orientation function, equal to the fraction of the bar area covered by yellow pixels, , as a function of ψb, (blue curve). Also shown are the ideal orientation function, (orange curve; see Sect. 3.4) and the angular window of width ∆ψ (Eq. (1)). Bottom left: subregion Ai of map A centered on pixel i, with the model bar overplotted in green at the orientation angle (ψf)i of a potential filament (corresponding to the peak of ).

3 Detailed methodology

In this section, we describe the successive steps of the FilDReaMS method in more detail, applying them to a map of intensity I. As explained at the beginning of Sect. 2.1, the word “intensity”, used in connection with the symbol I, is to be understood in a broad sense, which includes quantities such as column density and temperature. To illustrate the method, we provide detailed figures that rely on an H₂ column density ( ) map of the Herschel G210 field.

3.1 Spatial filtering

Let us start with a map of intensity I, which we refer to as the initial map A (top left panel of Fig. 1). Our purpose is to identify filaments in this map that can be locally approximated by a rectangular bar of width W_b.

The first step of FilDReaMS is to filter out structures wider than W_b. The spatial filtering is performed with the help of a 2D top-hat kernel of radius R, whose value is adjusted to the value of W_b in the manner explained after the following paragraph.

To begin with, the initial map A is convolved with the 2D top-hat kernel to produce a smoothed map B (top-middle panel of Fig. 1). Roughly speaking, this smoothing removes structures with widths smaller than ~2R. The smoothed map B is then subtracted from the initial map A to produce a map A–B from which structures with widths larger than ~2R are removed. Finally, the map A–B is transformed into a binary map C by setting all the pixels with positive values to 1 (yellow pixels in the right panels of Fig. 1) and all the pixels with negative values to 0 (dark pixels).

The adjustment of R to W_b is done iteratively by considering increasing values of R, starting at R = 2 px. For each value of R, the initial map A is convolved with a kernel of radius R (as explained above), and the binary map C is examined in search of regions of nonzero pixels wider than W_b. In practice, this is done by convolving C with a 2D top-hat kernel of diameter (W_b + 1 px); when the normalized convolution reaches a value of 1 at any pixel, we may conclude that this pixel is the center of a disk of diameter (W_b + 1 px) filled with nonzero pixels, which in turn implies that C contains a region of nonzero pixels wider than W_b. This is therefore where the iteration process stops. Thus, the binary map C represents the contrasted structures (yellow in Fig. 1) from the initial map A that are not wider than W_b.

The convolution gives rise to border effects, with a blank band of width R adjacent to the border of the convolved map B. As a result, map B is somewhat smaller than the initial map A. For large values of W_b, the blank band may represent a significant fraction of map A; in that case, large filaments close to the edge of map A may escape detection.

3.2 Orientation distribution

Let us consider a given pixel i in the full binary map C (top right panel of Fig. 1) as well as the surrounding subregion C_i of size L_b, centered on pixel i (bottom right panel). Pixel i must clearly be more distant than L_b /2 from the border of map C. Our purpose is to find structures in C_i that can be locally matched to a model bar of width W_b.

We consider values of the orientation angle of the model bar, ψ_b (defined with the conventions described in Sect. 2.1), spanning the range −90° to +90° in 1° steps. For each value of ψ_b, the model bar is centered on the considered pixel i, and we measure the fraction f_i (ψ_b) of the bar area covered by nonzero pixels (yellow pixels). This gives us the “measured” orientation function, (blue curve in the bottom-middle panel of Fig. 1).

The computation of the area resulting of the intersection of every pixel and the model bar rotated by ψ_b is performed only once at the beginning of the processing for all pixels of the L_b × L_b domain, and obtained through a numerical drizzling approach. In practice, each original pixel is subdivided into N_drizz subpixels, allowing us to compute an estimate of the exact intersection of the bar and the pixel at a desired accuracy. We adopt a value of N_drizz=101, scaling with 5/W_b, so that we keep an accuracy of better than 0.2% on the computation of the intersecting area with the model bar. Once computed, this kernel defined over L_b × L_b can be translated at any pixel i location and multiplied with the full binary map to obtain the “measured” orientation function for this pixel.

3.3 Detection of potential bar-like filaments

The measured orientation function, , makes it possible to detect potential filaments of bar width W_b centered on pixel i and to estimate their orientation angle, (ψ_f)_i. We stress that the measured orientation function is relatively smooth over the range of orientation angles, because it results from a kind of convolution with an extended bar at each orientation angle value, so that the determination of local maxima is not affected by local pixel-scale fluctuations (see bottom-right panel of Fig. 2).

FilDReaMS identifies the local maxima of the measured orientation function, , through an iterative process. It assigns the first peak to the maximum value of over the whole range of orientation angles. For this peak, FilDReaMS then assigns an angular window of width Δψ centered on the peak orientation angle, and corresponding to the expected domain of angular extension of the model bar, defined as twice the angle between the diagonals of the bar:

(1)

where r_b = L_b/W_b is the aspect ratio of the model bar.

Once the first peak is identified with its own angular window, FilDReaMS looks for any other local maximum over the remaining angular domain of the original interval, and assigns to this second peak the corresponding angular orientation and its Δψ angular window. This procedure is repeated until all the local maxima in are identified. In the event that two consecutive peaks, at ψ₁ and ψ₂, are so close that their windows overlap (|ψ₁ – ψ₂| < Δψ), they are considered to actually be part of one and the same potential filament; in that case, only the higher peak is retained together with its window, while the lower peak is ignored. All the peaks remaining after the handling of overlapping windows correspond to potential filaments, whereas the other peaks are considered to be part of the background.

In the case of Fig. 2, FilDReaMS detects two potential filaments with orientation angles ψ₁ and ψ₂ (bottom-right panel), while the common practice with RHT orientation functions is to compute the expectation over a given threshold leading to a single preferred orientation in this specific case (bottom-left panel).

Fig. 2 Illustration of a case where the considered pixel i lies at the intersection of two filaments with different orientations. Top: subregion of the initial map A, highlighting the filaments detected with FilDReaMS (dark blue lines) and when computing expectation over RHT orientation function (red line). Bottom: orientation function obtained with RHT (left) and FilDReaMS (right), showing the filament orientation angles derived in both cases.

3.4 Reliability assessment of potential bar-like filaments

3.4.1 The significance criterion

To test the reality of a potential filament detected at the considered pixel i, we compare the measured orientation function, , to the ideal orientation function, (right panel of Fig. 3), that would be obtained for an ideal filament – similar to a filament with the exact same shape as the model bar – superposed on an empty background, at the orientation angle (ψ_f)i of the potential filament (left panel). Clearly, at ψ_b = (ψ_f)i, the model bar coincides exactly with the ideal filament, such that all the pixels of the model bar have a value of 1 in the binary map of the ideal filament (middle panel). As a result, reaches its peak value, , at ψ_b, = (ψ_f)i (right panel). Quite expectedly, the width of the peak in is ~Δψ (Eq. (1)).

FilDReaMS then compares the measured and ideal orientation functions over an angular window of width Δψ centered on the peak associated with the potential filament. This window is both broad enough to capture the relevant characteristics of the potential filament and narrow enough to avoid contamination by a nearby filament at a slightly different angle. The comparison is performed with the help of the parameter

(2)

where

(3)

Ν_ψ is the number of ψ_b, values within the comparison window Δψ, and σ_f is the intrinsic error in (ψ_b) defined as the median absolute deviation computed over all the angles ψ_b, of all the pixels i of the initial map A. , which is actually similar to the square root of a reduced chi-squared, quantifies how close a potential filament is to a rectangular bar in an empty background.

We assess the significance of the potential filament by comparing to a threshold, (χ_r)_th, derived through Monte-Carlo simulations. As these simulations need to be run only once (for any given W_b) for the entire map A, not for each pixel i separately, we describe them in a separate subsection (Sect. 3.4.2).

We consider that a potential filament detected at pixel i is significant if

(4)

or, equivalently,

(5)

where

(6)

is defined as the significance of the detection. The definition of the threshold (χ_r)_th in Sect. 3.4.2 implies that there is a 5% chance of mistakenly rejecting an ideal filament that was, in reality, significant.

Fig. 3 Illustration of one ideal filament and its ideal orientation function. Left: map of an ideal filament (filament with the exact same shape as the model bar) centered on pixel i and oriented at angle (ψf)i (here (ψf)i < 0). Middle: corresponding binary map, with the model bar overplotted in green. Right: ideal orientation function, from the second step of FilDReaMS applied to the binary map.

3.4.2 Description of the Monte-Carlo simulations

The purpose of the Monte-Carlo simulations (illustrated in Fig. 4) is to derive the threshold, (χ_r)th, below which a potential filament can be considered significant against the background (see Eq. (4)). This threshold depends only on the bar width, W_b, and it applies to the entire map A. At each iteration (top row), a pixel j is drawn at random in map A, and a subregion A_j of size L_b, centered on pixel j, is cut out (leftmost panel of Fig. 4). Pixel j must lie far enough from the border of map A to ensure that A_j is entirely contained within A and that convolution with a 2D top-hat function of radius R will be possible (as explained in Sect. 3.1). A synthetic map is then created (middle-left panel) by superposing an ideal filament with uniform intensity I₀ on subregion A_j, centered on pixel j and oriented at random (over a uniform distribution of angles). For convenience, I₀ is written in terms of the standard deviation of subregion A_j, , and the signal-to-noise ratio of the filament, . (S/N)_fil is a free parameter, the value of which is discussed in Sect. 3.7.

Applying FilDReaMS to the synthetic map (surrounded by a band of width R from map A to allow convolution with a 2D top-hat function of radius R²) leads to a synthetic binary map (middle right panel of Fig. 4) followed by a synthetic orientation function, (blue curve in the right panel of Fig. 4). Adding a contrasted synthetic filament in implies that all the pixels close to the filament are now part of a less contrasted background (with respect to the filament), and together form a thin dark region surrounding the filament in (see Sect. 3.1).

The synthetic orientation function can be compared to the corresponding ideal orientation function, , similar to the orientation function obtained for the same ideal filament superposed on an empty background (orange curve). The associated ; is then derived from Eq. (2) with f^M replaced by f^S. Clearly, quantifies the impact of the background on the ideal filament by comparing the orientation functions f^I and f^S obtained from maps without (e.g., map A_i in Fig. 3) and with (map in Fig. 4) background, respectively.

This Monte-Carlo iteration is repeated 10000 times (each time with a new random region A_j and a new random orientation of the ideal filament). The resulting Monte-Carlo distribution D of is shown in the bottom panel of Fig. 4. The threshold (χ_r)_th is chosen to be the 95th percentile of the distribution of .

Fig. 4 Illustration of our Monte-Carlo simulations for the computation of the threshold (χr)th for a given value of the bar width, Wb. Leftmost: initial map A of the G210 field, with a subregion Aj centered on an arbitrary pixel j. Middle left: synthetic map, , composed of subregion Aj plus an ideal filament centered on pixel j. Middle right: binary map, , from the first step of FilDReaMS applied to , with the model bar overplotted in green. Rightmost: synthetic orientation function, (blue curve), from the second step of FilDReaMS applied to , together with the ideal orientation function, (orange curve). Bottom: Monte-Carlo distribution D of the parameter that compares fS and fI within the angular window of width Δψ (see Eq. (2) with superscript M replaced by S), and threshold (χr)th corresponding to a p-value of 5%.

3.5 Reconstruction of physical filaments

Above, we have used a model bar with given width W_b, and have applied FilDReaMS to a given pixel i of the initial map A. More specifically, we have detected potential bar-like filaments centered at pixel i (Sect. 3.3) and we have retained the significant filaments, similar to filaments with significance S_i > 1 (Sect. 3.4).

To detect all the bar-like filaments of bar width W_b in map A, we repeat the procedure for all the pixels i of map C that are more distant than L_b/2 from the border (see beginning of Sect. 3.2). The significance at all the pixels of C (for W_b = 14 px) is shown in the top panel of Fig. 5, while the significance at the centers of significant filaments (S > 1) is shown in the bottom panel. We can see that most pixels of map C are not significant (S < 1). The most significant filaments have S ≃ 6.

By combining all the significant bar-like filaments, we can then reconstruct the true shape and the intensity of physical filaments, as illustrated in Fig. 6. We first produce a binary map C′ in which the model bars associated with all the significant filaments have their pixels set to 1, while all the other pixels are set to 0. We then create a filament mask by multiplying the binary map C introduced in Sect. 3.1 by the new binary map C′, and we apply this mask to the initial map A, by computing a simple product of both images. The resulting map R (rightmost panel of Fig. 6) reveals the network of physical filaments of bar width W_b.

At this point, a filament orientation angle (denoted by (ψ_f)i at pixel i) is defined only at pixels at which one or more significant bar-like filaments are found (bottom-left panel of Fig. 6). We now assign a filament orientation angle (denoted by , at pixel i′) to all nonzero pixels of map R (rightmost panel of Fig. 6). For each pixel i′, we consider all the significant filaments whose associated model bars pass through i′, and we define , as the orientation angle of the most significant filament, similar to the filament with the highest significance, S_i (Eq. (6)). It is important to realize that for pixels i′ at which (ψ_f)_i′ is defined, , generally differs slightly from , because the most significant bar-like filament passing through i′ is generally not the filament centered on i′, but a filament centered on a neighboring pixel i. It must also be kept in mind that , is a function of the bar width, W_b.

Fig. 5 Maps of the significance, S (Eq. (6)), of the G210 Herschel field for a bar width Wb = 14 px: (top) for all the pixels of the binary map C and (bottom) for the pixels that are the centers of significant filaments (S > 1 ).

3.6 Derivation of the most significant bar widths and most prevalent bar widths

Now that we have laid out the entire procedure for a given bar width, W_b, (Sects. 3.1 to 3.5), we iterate over a range of values of W_b (see Sect. 3.7 for a recommendation on the optimal range). This enables us to assign a “most significant bar width”, , to every nonzero pixel i′ of the different reconstructed maps R(W_b). Namely, for every pixel i′, we consider all the significant filaments whose associated model bars pass through i′, and we define , as the bar width of the most significant filament.

We can then construct a histogram of over all pixels, namely the number of pixels, N_pix, whose most significant bar width is equal to . In practice, to make it easier to compare different maps, we propose to work with a normalized histogram, where N_pix is divided by the total number of pixels in the map, N_map. If the histogram exhibits clear peaks (i.e., local maxima), the corresponding bar widths are denoted by and are referred to as the “most prevalent bar widths”. Examples of histograms (N_pix/N_map) versus are shown in Figs. 9 and 10.

3.7 Instructions for use

We now discuss the optimal values, or ranges of values, of the three important free parameters of FilDReaMS (see Table 2): the width of the model bar, W_b, the aspect ratio of the model bar, r_b = L_b/W_b, and the signal-to-noise ratio of the ideal filament in the Monte-Carlo simulations, (S/N)_fil.

The aspect ratio of the model bar, r_b, sets an approximate lower limit to the aspect ratio of elongated structures that can be detected with FilDReaMS. Here, we consider that an elongated structure can be qualified as a filament if its aspect ratio is at least 3. Accordingly, we recommend that users adopt r_b, = 3. This value was also used by Panopoulou et al. (2014); Arzoumanian et al. (2019).

The signal-to-noise ratio of the ideal filament in the Monte-Carlo simulations, (S/N)_fil, sets an approximate lower limit to the intensity I (in a broad sense) of filaments that have a high probability of being detected with FilDReaMS. If (S/N)_fil is too low, some of the small-scale noise structures might be mistaken for physical filaments. On the other hand, if (S/N)_fil is too high, some physical filaments might escape detection. As a compromise, we recommend adopting (S /N)_fil = 3.

The width of the model bar, W_b, is only constrained by the pixel size at the low end and by the size of the initial map A at the high end. For the lower limit, a good choice is typically (W_b)_min = 5 px; this choice is particularly relevant in the case of the Herschel fields, where the beam size equals three times the pixel size. For the upper limit, we take (W_b)_max = (L_b)_max/r_b, with (L_b)_max equal to one-third the size of map A; in the Herschel G210 field, (W_b)max = 27 px.

4 Validation

FilDReaMS is in principle able to detect filaments of different sizes and orientations in a given map of intensity I. Before applying FilDReaMS to real scientific data, we need to validate the method and estimate its reliability.

In the following, we explore the impact of the filament profile, the noise level, and the aspect ratio on the bar width and orientation estimates. We finally investigate the case with superposition of filaments with variable intensities, in combination with the above parameters.

We first describe the set of simulations used to performed these analyses in Sect. 4.1. We then test the ability of FilDReaMS to recover the widths of the input filaments in Sect. 4.2, and their orientations in Sect. 4.3.

4.1 Set of simulations

We create several series of simulated maps composed of synthetic filaments embedded in realistic environments, including noise and incoherent structures. To study the impact of the filament radial profile, we consider in these simulations two types of synthetic filaments that by default (i.e., unless explicitly stated otherwise) have the following characteristics:

• Ideal filaments: they have the shape of the rectangular model bar, with width W_b, length L_b, and aspect ratio r_b = L_b/W_b, and they have uniform intensity, I₀.

• Plummer-type filaments: their intensity can be described by a 2D Plummer-type profile with half-width R_flat and aspect ratio rb:

(7)

where x and y are the coordinates across and along the long axis of the filament, p is the Plummer power-law index, and I₀ is the central intensity. We adopt p = 2.2 (median value obtained by Arzoumanian et al. 2019, for a sample of 599 filaments including G300).

The realistic environment is simulated as a combination of white and power-law (Brownian) noises, with root mean square (rms) intensity σ_W and σ_B, respectively. The white noise is meant to represent instrumental noise, which was estimated to be ≤7% of the intensity signal for Herschel SPIRE data (Juvela et al. 2012), while the Brownian noise (with a power-law index of −2) is meant to represent astrophysical fluctuations (e.g., Miville-Deschênes et al. 2003). The two types of noise have different realizations in each map. In the following, we consider two different configurations denoted as “default” and “high” levels of noise. In the first case, we chose σ_W = 0.05 I₀ and σ_Β = 0.3 I₀, while in the second case we chose σ_W = 0.1 I₀ and σ_Β = I₀.

We perform four different sets of simulations:

• SimSet 1: r_b = 3 and default noise level.

• SimSet 2: r_b = 3 and high noise level.

• SimSet 3: r_b = 10 and default noise level.

• SimSet 4: Overlapping filaments with variable aspect ratios and intensities, with default noise level.

For the first three sets of simulations, the configuation consists of the co-addition of a map of synthetic filaments with realizations of realistic noise. Each original map is made of 91 nonoverlapping synthetic filaments of common width, W_b, and aspect ratio, r_b, but different orientations spanning linearly from 0^° to +90^° in 1^° steps. A series of 16 maps is thus created for each value of the width, W_b, which ranges from 5 to 20 px (in steps of 1 px), on which 2000 realizations of realistic noise are co-added.

The last set of simulations, SimSet 4, consists of a single series of 100 maps containing 64 synthetic filaments: 48 filaments with W_b = 9px and 16 filaments with W_b = 15 px. Both groups of filaments cover approximately the same total surface area. Each filament is randomly assigned one of eight possible values of the filament aspect ratio, r_b, in the linear range [3,10] and (independently) one of eight possible values of the filament intensity, I₀, in the logarithmic range [0.2, 5]. Each value of r_b and each value of I₀ is assigned to exactly eight filaments. The filaments are located at random positions, have random orientations, and can possibly overlap. For each map realization, we co-add a noise realization using the default configuration.

All sets of simulations are finally repeated for ideal and Plummer-type filaments, and their simulation parameters are displayed in Table. 3. In the case of SimSet 1, we also produce two other sets of simulations using values of the Plummer-type filaments with power-law index, p = 1.2 and 4. The value p = 4 corresponds to an isolated filament in isothermal, nonmagnetic hydrostatic balance (Ostriker 1964); this value was also used by Suri et al. (2019) to fit filament profiles. The value p = 1.2 is the lower limit explored by Arzoumanian et al. (2019) in their study. The larger and smaller values yield more and less contrasted filaments, respectively.

Fig. 6 Illustration of the method used to reconstruct physical filaments of bar width (Wb = 14 px). Top left: initial map A of the G210 field. Top middle: binary map C from the first step of FilDReaMS applied to map A (same pixels as in the top panel of Fig. 5). Bottom left: pixels at which at least one significant bar-like filament of bar width Wb is detected (same pixels as in the bottom panel of Fig. 5). Bottom middle: binary map C′ in which the model bars associated with all the significant bar-like filaments are set to 1 and the background is set to 0. Right: filament mask obtained by multiplying C by C′. Rightmost: reconstructed physical filaments obtained by applying the filament mask to the initial map A.

4.2 Filament bar widths

In the following, we construct the histogram (N_pix/N_map) versus of every map (as explained in Sect. 3.6), and we average over the 2000 noise realizations for the first three sets of simulations (SimSet 1, 2, and 3) and over the 100 noise realizations for the fourth set of simulations (SimSet 4). For this analysis, we can identify a most prevalent bar width , averaged over the multiple noise realizations, and its uncertainty (or precision), which is equivalent to the standard deviation of the histogram. The most prevalent bar width is then compared with the input width of the simulation.

4.2.1 Impact of the noise level

For this study, we focus on the first two sets of simulations, SimSet 1 and 2.

In the case of ideal filaments, we find for each input value of W_b that the average histogram is dominated by a well-defined peak at , which means that FilDReaMS exhibits a very good accuracy by easily recovering the right width of the input filament, even in the presence of high noise. However, the peak does not stand out as clearly at the high noise level (SimSet 2) as at the default noise level (SimSet 1). The precision of the width estimate is very stable and is only moved from 0.01 px to 0.4 px with increasing noise level (see Table 4). This is illustrated for W_b = 12 px in Fig. 7, where the histogram at the high noise level (middle panel, blue) can be compared to the histogram at the default noise level (top panel, blue).

We show the case of the Plummer-type filaments for the same sets of simulations in the top and middle panels of Fig. 7 (orange histograms), with R_flat = 12 px. In the case of default noise level (top panel), the histogram of is clearly shifted compared to the input value and more spread than in the case of the ideal filament. This shift is expected and will depend on the power-law of the Plummer profile, as illustrated in Sect. 4.2.3. The larger spread can be explained here by the larger impact of the noise on the tails of the Plummer profile. In addition, the central intensity I₀ is only reached along the crest of the Plummer-type filament. The integral of the transverse profile of a Plummer-type filament over a bar width W_b = 2 R_flat represents 75% of the integration over the corresponding ideal filament (when computed with the default Plummer power-law index equal to 2.2). Therefore, there will be a slightly higher degradation from the noise as one moves away from the crest. Despite the shift observed in the case of the Plummer-type filament with default noise, we still obtain a precision of 0.4 px. However, the histograms (N_pix/N_map) versus of Plummer-type filaments may exhibit one or two spurious (lower) peaks.

The impact of increasing the noise level is much more pronounced than in the ideal filament case (see middle panel), leading to much larger uncertainties, 2.7 px, stable over the whole range of input width (see Table 4). The increase in noise level increases the frequency and the strength of the spurious peaks. The observed shift between the widths 2 R_flat and can be empirically predicted, as we discuss in Sect.4.2.3.

Fig. 7 Average histograms (Npix/Nmap) versus obtained for ideal (blue) and Plummer-type (orange, with power-law index p = 2.2) filaments with bar width Wb = 12 px and Plummer half-width Rflat = 12 px, respectively, in three different sets of simulations, SimSet 1 (top), Sim-Set 2 (middle) and SimSet 3 (bottom).

4.2.2 Impact of the filament aspect ratio

For this analysis, we use the two sets of simulations SimSet 1 and 3, between which the aspect ratio has been changed from 3 to 10. The peak becomes increasingly dominant as r increases, leading to a gain in precision by a factor of between 2 and 5, for Plummer and ideal cases, respectively (see Table 4). This is in accordance with our expectation that more elongated filaments are easier to recognize. In the case of Plummer-type filaments, the shift between 2 R_flat and appears to be unchanged. This is illustrated in Fig. 7, where the histograms obtained for r_b = 10 (bottom panel) can be compared to the histograms obtained for r_b = 3 (top panel).

4.2.3 Plummer-to-ideal width relation

The sets of simulations with Plummer-type filaments allow us to derive the empirical relation between 2 R_flat and for different noise levels and also for different values of the Plummer index, p. The curves linking the input 2 R_flat to the output are plotted in the top panel of Fig. 8 for the default and high noise levels, SimSet 1 (blue triangles) and SimSet 2 (dark blue triangles), respectively. Also plotted are the line = 2 R_flat (dashed line), the linear fits to the two curves versus 2 R_flat (solid lines in the corresponding colors), as well as the standard deviations computed on each average histogram ( Npix/ Nmap) versus (error bars on each symbol in the corresponding colors). For reference, the parameters a and b of the linear fits (2 R_flat) + b are given in Table 5.

The error in the fit for SimSet 1 can be roughly estimated at ≲1px over the range R_flat = [5,20] px, from a comparison with linear fit constrained to pass through the origin. This error adds to the statistical uncertainty ≃0.5 px associated with the dispersion of the measured points about the fits and to the standard deviations (see Table 4). In the case of high noise level (SimSet 2), a slightly decreases with increasing noise, whereas b takes on larger values. This, combined with the high statistical uncertainty and standard deviations (error bars), indicates that FilDReaMS is not well suited to recovering the width of input Plummer-type filaments in cases where there are high levels of noise.

For completeness, we also look at the impact of the Plummer index on the curves versus 2 R_flat. Plotted in the bottom panel of Fig. 8 are the curves obtained for p = 1 . 2 (orange squares), p = 2.2 (blue triangles), and p = 4 (green circles) with their respective standard deviations (error bars in the corresponding colors). The linear fits to the three curves versus 2 R_flat are again shown in solid lines, and their parameters given in Table 5. The parameter a increases with decreasing p. The increase in slope can be understood by noting that, for a given R_flat, a smaller p yields a more spread-out Plummer-type filament (see Eq. (7)), which in turn results in a detection at larger . As for p = 2.2, the error in the fits for the other two values of p can be roughly estimated at ≲1 px over the range R_flat = [5,20] px.

Finally, the curves versus 2 R_flat obtained for increasing values of r_b are identical, with the standard deviations (error bars in Fig. 8) decreasing down to half those obtained for r_b = 3. From this, we conclude that the empirical relation between 2 R_flat and derived in SimSet 3 remains valid independently of r_b.

Fig. 8 Average relation between the width 2 Rflat of an input Plummer-type filament and the most prevalent bar width derived by FilDReaMS, for two different noise levels (top) and for three different values of the Plummer index, p (bottom). The solid lines are linear fits to the relations versus 2 Rflat of the same colors, and the dashed line is the straight line .

Fig. 9 Visualization of the set of simulations SimSet 4 for ideal filaments, showing how well FilDReaMS recovers the widths Wb of overlapping filaments with different widths, aspect ratios, and intensities. Top left: one of the 100 test maps, containing 64 ideal filaments, with width Wb = 9 px or 15 px, aspect ratio in the range [3, 10], and intensity in the range [0.2,5]. Top right: corresponding map of all the filaments detected and reconstructed by FilDReaMS. Bottom: number of pixels, Npix, whose most significant bar width is , normalized to the number of pixels in the map, Nmap, as a function of .

4.2.4 Overlapping filaments of different properties

To study the reliability of the width estimate in a more realistic case, we focus now on the fourth set of simulations, SimSet 4, for ideal and Plummer-type filaments. In the case of ideal filaments, FilDReaMS is able to detect and reconstruct all of them (top right panel of Fig. 9). The histogram (N_pix/N_map) versus , averaged over the 100 maps, has two pronounced peaks at px and (bottom panel of Fig. 9). Hence, FilDReaMS recovers again the widths of the input ideal filaments. The histogram also contains non-negligible values away from the peaks because the superposition of two input filaments can either lead to a structure that will be identified as a thicker filament or conversely cause the broader filament to effectively hide part of the narrower filament.

In the case of Plummer-type filaments, FilDReaMS is also able to detect and reconstruct all of them (top right panel of Fig. 10). The histogram (N_pix/N_map) versus , averaged over the 100 test maps, has three peaks at , and 19 px that present a larger spread compared to the ideal case (bottom panel of Fig. 10). The two peaks at and 11 px are associated with the same input filament width (R_flat = 9px). This can be explained by the log-scale distribution of the central intensity of the filaments, I₀, which leads to two different noise regimes, with one-third of the filaments in the high noise level, and two-thirds in the default noise level. Following top panel of Fig. 8, a 2 R_flat value of 18 px roughly corresponds to bar widths of W_b = 9 and 11 px in the high and default noise levels, respectively, which is consistent with the two observed peaks. The same kind of behavior is not observed in this simulation for the population of larger filaments, because they are more affected by the impact of overlapping filaments, leading to a spread distribution masking this effect.

The corresponding values of 2 R_flat can be derived from the linear fits , with the values of α and b given in Table 5. We use the case p = 2.2 and the high noise level for the first peak at 9 px (as explained above), and the case p = 2.2 with default noise level for the others. This leads to 2 R_flat = 16.3 px, 17.3 px, and 28.5 px, in good agreement with the input R_flat (9 px and 15 px). This agreement supports the general validity of the empirical relation between 2 R_flat and obtained in Sect. 4.2.3. As in the case of ideal filaments, the fact that the histogram contains more than the two expected peaks except the second peak at 9 px is not an indication that FilDReaMS performs poorly, but rather a consequence of the subjective way of identifying filaments.

Fig. 10 Same as Fig. 9, with the input ideal filaments of width Wb replaced by Plummer-type filaments of width 2 Rflat.

4.3 Filament orientations

Again, to study the impact of noise level and aspect ratio on the reliability of the filament orientation estimates, we use the first three sets of simulations SimSet 1, 2, and 3 for ideal and Plummer-type filaments. Hence, for each series, the derived orientation angles, , are compared to the input orientations, ψ_f over the 2000 noise realizations of each set of simulations. This allows us to estimate both the bias and the standard deviation σ_ψ of .

We first observe that the orientations of the reconstructed filaments do not show any deviation on average from the input orientation angles, and do not exhibit any specific dependency with the input value of the orientation angle. We also observe from these simulations that σ_ψ decreases with the filament width, meaning that the results obtained with W_b and R_flat = 5 px can be considered as upper limits. We consider these conservatives values as our final uncertainty on the filament orientation, displayed in Table 4.

We can see that the precision is very good, that is, ≤0.2° for ideal filaments and ≤1.8° for Plummer-type filaments for the set of simulations SimSet 1 (default noise level and short aspect ratio). These uncertainties increase with increasing noise (Sim-Set 2), up to ≤0.4° (ideal) and ≤3.1° (Plummer). The impact of noise level on the filament orientation estimate is much less important than that observed for the width, with a degradation of only a factor of ≤2 for both kinds of synthetic filament profiles.

The uncertainties decrease with increasing aspect ratio (Sim-Set 3), down to ≤0.1° (ideal) and ≤1.2° (Plummer), which is again in accordance with our expectation that more elongated filaments are easier to recognize.

Fig. 11 Application of FilDReaMS to the 250 μm intensity map of the Herschel G210 field. Top row: 250 μm intensity map of G210. Second row: reconstructed filaments with bar widths Wb = 6px (left), Wb = 10 px (middle), and Wb = 17 px (right). Third row: map of all the reconstructed filaments with bar width in the range [5,27] px. Bottom row: filaments detected with RHT using a 2D top-hat radius R = 11 px (Malinen et al. 2016). The color code of the reconstructed filaments indicates their column density .

5 First astrophysical application and comparison with previous work

As a first astrophysical illustration, we apply FilDReaMS to the Herschel 250 μm intensity map of the G210 field (top panel of Fig. 11). The left, middle, and right panels in the second row of Fig. 11 show the networks of reconstructed filaments with bar widths W_b = 6 px, 12 px, and 17 px. The corresponding radius of the 2D top-hat kernel used to filter out the large scales in the initial map (to obtain map B; see top middle panel of Fig. 1) is R = 11 px, 24 px, and 45 px, respectively. We immediately see that the majority of the small filaments are substructures of larger ones, while only a small fraction of them cover areas not already included in larger structures, stressing the hierarchical organization of these filamentary structures in the interstellar medium.

The map of all the reconstructed filaments with bar widths between W_b = 5 px and 27 px is displayed in the third row of Fig. 11. This map bears some obvious resemblance to the initial map (top row), while a direct inspection by eye indicates that FilDReaMS recovers most of the filamentary structures independently of the intensity.

The striations in the northern-central part of the region are better reconstructed (i.e., over longer path lengths) when integrating all bar widths compared to when using any single value of W_b. These striations are in fact composed of filaments of different bar widths, mostly parallel to each other, with the thinner corresponding to the crests of the larger.

We now compare the above results to those of Malinen et al. (2016), who applied the RHT method to the Herschel 250 μm intensity map of G210 (top panel of Fig. 11), using a spatial filtering with R = 11 px. The filaments detected with RHT are displayed in the bottom row of Fig. 11, where they can be compared to the filaments detected with FilDReaMS (left panel in the second row). It appears that both methods yield similar filamentary networks, but with a few noticeable differences: RHT detects more filaments in the most diffuse part of the region, while in the densest part FilDReaMS detects more strands and structures branching out from the main central filament. The long and thin striations that seem to escape detection with FilDReaMS for R = 11 px could be missed because of the short aspect ratio of the FilDReaMS model bar. However, they are recovered when considering a short range of values around R = 11 px.

6 Concluding remarks

In this paper, we present a new method to detect filaments of different widths in a map of intensity (in a broad sense) I. We call our method Filament Detection and Reconstruction at Multiple Scales, or FilDReaMS. Briefly, FilDReaMS uses a rectangular model bar of width W_b and aspect ratio r_b, where W is meant to cover a broad range of values and r_b is a free parameter (typically set to r_b = 3). For any given value of W_b, FilDReaMS (1) detects potential filaments that can be locally approximated by the model bar, (2) retains the potential filaments with significance S > 1 (Eq. (5)), (3) reconstructs the true shape and the intensity of physical filaments from the initial map of I together with the associated binary map (see Fig. 6), and (4) assigns a filament orientation angle, , to each pixel i of a reconstructed filament of width W_b. After repeating the procedure for all the values of W_b, a most significant bar width, , is derived for each pixel i of all the reconstructed filaments, and the most prevalent bar widths of the map, , are inferred from the peaks of the histogram of . The can then be cautiously converted to the widths of the often-used Plummer-type profiles, 2 R_flat (see Fig. 8 and Table 5).

Thus FilDReaMS makes it possible to detect filaments of a given bar width, to identify the most prevalent bar widths (and corresponding Plummer widths) in a given map, and to derive the local orientation angles of the detected filaments. The main assets of FilDReaMS are the ability to detect filaments over a broad range of widths; the small number of free parameters (only three; see Table 2); the speed of execution: typically, for a given map, running FilDReaMS takes about 20–30s for each value of W_b and roughly 10–20 min to cover the entire range of W_b; and its user-friendliness, which makes it particularly suited for statistical studies.

FilDReaMS opens wide perspectives of application to astro-physical data, and can be used to study the processes of stellar formation inside interstellar filamentary structures. As an initial test case, we investigate the interplay between the orientation of the filaments and the Galactic magnetic field in a companion paper (Carrière et al. 2022), generalizing the analysis of Malinen et al. (2016) to four other fields. A broader statistical analysis over 116 Herschel fields is also in preparation.

Acknowledgements

We extend our deepest thanks to our referee, Gina Panopoulou, for her careful reading of our paper and for her many constructive comments and suggestions. We also acknowledge useful discussions with Dana Alina, Susan Clark, Mika Juvela, and Julien Montillaud. Herschel SPIRE has been developed by a consortium of institutes led by Cardiff University (UK) and including University Lethbridge (Canada); NAOC (China); CEA, LAM (France); IFSI, University Padua (Italy); IAC (Spain); Stockholm Observatory (Sweden); Imperial College London, RAL, UCL-MSSL, UKATC, University Sussex (UK); Caltech, JPL, NHSC, University Colorado (USA). This development has been supported by national funding agencies: CSA (Canada); NAOC (China); CEA, CNES, CNRS (France); ASI (Italy); MCINN (Spain); SNSB (Sweden); STFC (UK); and NASA (USA).

References

All Tables

All Figures

Fig. 1 Illustration of the FilDReaMS method applied to one pixel i of the H2 column density map of the Herschel G210 field. Top left: initial map A of the entire G210 field. Top middle: smoothed map B obtained by convolving map A with a 2D top-hat kernel. Top right: corresponding binary map C in which the yellow pixels have a value of 1 and the dark pixels a value of 0. Bottom right: subregion Ci of map C centered on pixel i, with the model bar (width Wb, length Lb, and orientation angle ψb) overplotted in green. Bottom middle: measured orientation function, equal to the fraction of the bar area covered by yellow pixels, , as a function of ψb, (blue curve). Also shown are the ideal orientation function, (orange curve; see Sect. 3.4) and the angular window of width ∆ψ (Eq. (1)). Bottom left: subregion Ai of map A centered on pixel i, with the model bar overplotted in green at the orientation angle (ψf)i of a potential filament (corresponding to the peak of ).

In the text

Fig. 2 Illustration of a case where the considered pixel i lies at the intersection of two filaments with different orientations. Top: subregion of the initial map A, highlighting the filaments detected with FilDReaMS (dark blue lines) and when computing expectation over RHT orientation function (red line). Bottom: orientation function obtained with RHT (left) and FilDReaMS (right), showing the filament orientation angles derived in both cases.

In the text

	Fig. 3 Illustration of one ideal filament and its ideal orientation function. Left: map of an ideal filament (filament with the exact same shape as the model bar) centered on pixel i and oriented at angle (ψf)i (here (ψf)i < 0). Middle: corresponding binary map, with the model bar overplotted in green. Right: ideal orientation function, from the second step of FilDReaMS applied to the binary map.
In the text

Fig. 4 Illustration of our Monte-Carlo simulations for the computation of the threshold (χr)th for a given value of the bar width, Wb. Leftmost: initial map A of the G210 field, with a subregion Aj centered on an arbitrary pixel j. Middle left: synthetic map, , composed of subregion Aj plus an ideal filament centered on pixel j. Middle right: binary map, , from the first step of FilDReaMS applied to , with the model bar overplotted in green. Rightmost: synthetic orientation function, (blue curve), from the second step of FilDReaMS applied to , together with the ideal orientation function, (orange curve). Bottom: Monte-Carlo distribution D of the parameter that compares fS and fI within the angular window of width Δψ (see Eq. (2) with superscript M replaced by S), and threshold (χr)th corresponding to a p-value of 5%.

In the text

	Fig. 5 Maps of the significance, S (Eq. (6)), of the G210 Herschel field for a bar width Wb = 14 px: (top) for all the pixels of the binary map C and (bottom) for the pixels that are the centers of significant filaments (S > 1 ).
In the text

Fig. 6 Illustration of the method used to reconstruct physical filaments of bar width (Wb = 14 px). Top left: initial map A of the G210 field. Top middle: binary map C from the first step of FilDReaMS applied to map A (same pixels as in the top panel of Fig. 5). Bottom left: pixels at which at least one significant bar-like filament of bar width Wb is detected (same pixels as in the bottom panel of Fig. 5). Bottom middle: binary map C′ in which the model bars associated with all the significant bar-like filaments are set to 1 and the background is set to 0. Right: filament mask obtained by multiplying C by C′. Rightmost: reconstructed physical filaments obtained by applying the filament mask to the initial map A.

In the text

	Fig. 7 Average histograms (Npix/Nmap) versus obtained for ideal (blue) and Plummer-type (orange, with power-law index p = 2.2) filaments with bar width Wb = 12 px and Plummer half-width Rflat = 12 px, respectively, in three different sets of simulations, SimSet 1 (top), Sim-Set 2 (middle) and SimSet 3 (bottom).
In the text

	Fig. 8 Average relation between the width 2 Rflat of an input Plummer-type filament and the most prevalent bar width derived by FilDReaMS, for two different noise levels (top) and for three different values of the Plummer index, p (bottom). The solid lines are linear fits to the relations versus 2 Rflat of the same colors, and the dashed line is the straight line .
In the text

Fig. 9 Visualization of the set of simulations SimSet 4 for ideal filaments, showing how well FilDReaMS recovers the widths Wb of overlapping filaments with different widths, aspect ratios, and intensities. Top left: one of the 100 test maps, containing 64 ideal filaments, with width Wb = 9 px or 15 px, aspect ratio in the range [3, 10], and intensity in the range [0.2,5]. Top right: corresponding map of all the filaments detected and reconstructed by FilDReaMS. Bottom: number of pixels, Npix, whose most significant bar width is , normalized to the number of pixels in the map, Nmap, as a function of .

In the text

	Fig. 10 Same as Fig. 9, with the input ideal filaments of width Wb replaced by Plummer-type filaments of width 2 Rflat.
In the text

Fig. 11 Application of FilDReaMS to the 250 μm intensity map of the Herschel G210 field. Top row: 250 μm intensity map of G210. Second row: reconstructed filaments with bar widths Wb = 6px (left), Wb = 10 px (middle), and Wb = 17 px (right). Third row: map of all the reconstructed filaments with bar width in the range [5,27] px. Bottom row: filaments detected with RHT using a 2D top-hat radius R = 11 px (Malinen et al. 2016). The color code of the reconstructed filaments indicates their column density .

In the text