J.-S Carrière

L Montier

K Ferrière

I Ristorcelli

FilDReaMS 1. Presentation of a new method for Filament Detection and Reconstruction at Multiple Scales

Keywords: ISM, clouds -ISM, structures -ISM, magnetic fields -dust -infrared, ISM -submillimiter, ISM -techniques, image processing

Context. Filamentary structures appear to be ubiquitous in the interstellar medium. Being able to detect and characterize them is the first step toward understanding their origin, their evolution, and their role in the Galactic cycle of matter. Aims. We present a new method, called FilDReaMS, to detect and analyze filaments in a given image. This method is meant to be fast, user-friendly, multi-scale, and suited for statistical studies. Methods. The input image is scanned with a rectangular model bar, which makes it possible to uncover structures that can be locally approximated by this bar and to derive their orientations. The bar width can be varied over a broad range of values to probe filaments of different widths. Results. We performed several series of tests to validate the method and to assess its sensitivity to the level of noise, the filament aspect ratios, and the dynamic range of filament intensities. We found that the method exhibits very good performance at recovering the orientation of the filamentary structures, with an accuracy of 0.5 • in nominal conditions, up to 3 • in the worst case scenario with high level of noise. The filaments width is recovered with uncertainties better than 0.5 px (pixels) in most of the cases, which could extend up to 3 px in case of low signal-to-noise ratios. Some attempt to build a correspondence between Plummer-type filament profiles and outcomes of the method is proposed, but remains sensitive to the local environment. Conclusions. Our method is found to be robust and adapted to the identification and the reconstruction of filamentary structures in various environments, from diffuse to dense medium. It allows us to explore the hierarchical scales of these filamentary structures with a high reliability, especially when dealing with their orientation.

Introduction

A large variety of methods have already been developed to extract elongated structures in two-dimensional (2D) maps. Their 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, by 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 et al. 2016b) or the Hessian matrix (second-order derivatives) (e.g., Polychroni et al. 2013;Schisano et al. 2014;Planck Collaboration et al. 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 et al. 2016b,a). 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(Arzoumanian et al. , 2019;;Peretto et al. 2012;Palmeirim et al. 2013) and to 13 CO intensity maps (Panopoulou et al. 2014). A limitation of the local approach is the difficulty of detecting faint structures such as striations. 1Non-local methods focus on a given scale around each pixel. The template-matching approach (Juvela 2016) allows one 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 neighbourhood 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 multiwaveband 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 it remains quite fast. However, it requires additional steps to extract filament orientations. Although this method is 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. They found similarities, but also 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 detections in various fields for a multi-scale statistical analysis. While none of the methods described above is 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.

General overview

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. This map will be our reference to present the FilDReaMS method. For simplicity, we will refer to I as being 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 polarisation angles. Since 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 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(Juvela et al. , 2012)), the so-called G210 field, which corresponds to the high Galactic latitude star-forming region L1642. This region was in particular studied in detail by Malinen et al. (2016), who investigated the relative orientations between the magnetic field (traced with Planck polarisation 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 ± 20 pc (Montillaud et al. 2015). Its H 2 column density map has dimensions ∆l × ∆b = 1.28 • × 1.22 • , corresponding to 3.1 pc × 3.0 pc. The angular size of a pixel, equal to one-third of the beam size, is 12", corresponding to 0.0081 pc. The H 2 column density, N H 2 , varies over the range [0.2, 7.5] × 10 21 cm -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 θ H of its normal and the minimal distance ρ H to the origin), the Hough transform allows us to pass from the pixel space (x, y) to this other parameter space (ρ H , θ H ), 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 ρ H =0) inside a circular area defined around this central pixel in order to pass from pixel space to (θ H ) 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 over wich 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 Article number, page 2 of 15 J.-S. Carrière et al.: FilDReaMS. 1. A new method to detect filaments and determine their orientations Most prevalent bar width for the entire map 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 (ρ H , θ H ) space, directly implies that RHT distributions of orientation are dependent of the pixelisation 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.

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 θ H representation. FilDReaMS solves this problem by a funda-mental change of philosophy: it starts from a discretisation of the θ H 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 detections.

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.

The main steps of FilDReaMS

In brief, the main steps of FilDReaMS are the following: 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 i ′ , (ψ ⋆ f ) i ′ , 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 ′ , (W ⋆ b ) i ′ , as well as the most prevalent bar widths for the entire map, W ⋆peak b , see Sect. 3.6.

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

Detailed methodology

In this section, we describe in more detail the successive steps of the FilDReaMS method applied 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 2 column density (N H 2 ) map of the Herschel G210 field.

Spatial filtering

Let us start with a map of intensity I, which we will refer to as the initial map A (top left panel of Fig. 1). Our purpose is to identify in this map filaments 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 next 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 adjustement 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 non-zero 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 non-zero pixels, which in turn implies that C contains a region of non-zero pixels wider than W b . Hence 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.

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 sub-region C i of size L b , centered on pixel i (bottom right panel). Obviously, pixel i must 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 non-zero pixels (yellow pixels). This gives us the "measured" orientation function, f M i (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 sub-divided into N drizz sub-pixels, 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 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.

Detection of potential bar-like filaments

The measured orientation function, f M i , 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 angle, since 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, f M i , through an iterative process. It assigns the first peak to the maximum value of f M i over the whole range of orientation angle. It then assigns for this peak, 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 bar's diagonals:

∆ψ = 4 arctan 1 r b , (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 f M i are identified. In the event that two consecutive peaks, at ψ 1 and ψ 2 , are so close that their windows overlap (|ψ 1 -ψ 2 | < ∆ψ), 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 ψ 1 and ψ 2 (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).

Reliability assessment of potential bar-like filaments

The significance criterion

To test the reality of a potential filament detected at the considered pixel i, we compare the measured orientation function, f M i , to the ideal orientation function, f I i (right panel of Fig. 3), that would be obtained for an ideal filament -similar to a filament having 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, f I i reaches its peak value,

f I i = 1, at ψ b = (ψ f ) i (right panel).
Quite expectedly, the width of the peak in f I i 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

(χ M r ) i = χ 2 N ψ = 1 N ψ (ψ f ) i +∆ψ/2 ψ b =(ψ f ) i -∆ψ/2 f I i (ψ b ) -f M i (ψ b ) 2 σ 2 f , (2)
where

σ f = median f M i (ψ b ) -median f M i (ψ b ) , (3)
N ψ is the number of ψ b values within the comparison window ∆ψ, and σ f is the intrinsic error in f M i (ψ b ) defined as the median absolute deviation computed over all the angles ψ b of all the pixels i of the initial map A. (χ M r ) i , 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 (χ M r ) i to a threshold, (χ r ) th , derived through Monte-Carlo simulations. Since 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

(χ M r ) i < (χ r ) th , (4)
or, equivalently,

S i > 1 , (5)
where

S i = (χ r ) th (χ M r ) i (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 actually significant.

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 sub-region 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 A S j is then created (middle left panel) by superposing on the sub-region A j an ideal filament with uniform intensity I 0 , centered on pixel j and oriented at random (over a uniform distribution in angle). For convenience, I 0 is written in terms of the standard deviation of sub-region A j , σ A j , and the signal-to-noise ratio of the filament, (S/N) fil : I 0 = σ A j (S/N) fil . (S/N) fil is a free parameter, whose value is discussed in Sect. 3.7.

Applying FilDReaMS to the synthetic map A S j (surrounded by a band of width R from map A to allow convolution with a 2D top-hat function of radius R2 ) leads to a synthetic binary map C S j (middle right panel of Fig. 4) followed by a synthetic orientation function, f S j (blue curve in the right panel of Fig. 4). Adding a contrasted synthetic filament in A S j 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 the surrounding filament in C S j (see Sect. 3.1). The synthetic orientation function f S j can be compared to the corresponding ideal orientation function, f I j , similar to the orientation function obtained for the same ideal filament superposed on an empty background (orange curve). The associated (χ S r ) j is then derived from Eq. 2 with f M replaced by f S . Clearly, (χ S r ) j quan-tifies the impact of the background on the ideal filament by comparing the orientation functions f I and f S obtained from maps without (map A i in Fig. 3, for example) and with (map A S j in Fig. 4) background, respectively. This Monte-Carlo iteration is repeated ten thousand 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 χ S r is shown in the bottom panel of Fig. 4. The threshold (χ r ) th is chosen to be the 95-th percentile of the distriof χ S r .

Reconstruction of physical filaments

So far, we have used a model bar with given width W b , and we 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) are 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 (ψ ⋆ f ) i ′ at pixel i ′ ) to all non-zero 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 (ψ ⋆ f ) i ′ 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, (ψ ⋆ f ) i ′ generally differs slightly from (ψ f ) i ′ , 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. Also keep in mind that (ψ ⋆ f ) i ′ is a function of the bar width, W b .

Derivation of 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", (W ⋆ b ) i ′ , to every non-zero pixel i ′ of the different reconstructed maps R(W b ). Namely, for every pixel i ′ , we consider all the significant

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 Table 2. Free parameters of FilDReaMS, with their definitions (second column) and the corresponding characteristics of the filaments to be detected by FilDReaMS (third column). Here, we consider that an elongated structure can be qualified as a filament if its aspect ratio is at least 3. Accordingly, we recommend to adopt r 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 to adopt (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.

Validation

FilDReaMS is in principle able to detect filaments with 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 in Sect. 4.1 the set of simulations used to performed these analyses. 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.

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, which, 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 0 .

-Plummer-type filaments: Their intensity can be described by a 2D Plummer-type profile with half-width R flat and aspect ratio r b :

I P (x, y) = I 0        1 + x R flat 2 + y r b R flat 2        -(p-1)/2 , ( 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 0 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 r.m.s. 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) 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, non-magnetic 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.

Filament bar widths

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

Impact of the noise level

For this study we focus on the first two sets of simulations, Sim-Set 1 & 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 W ⋆peak b = W b , 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 (Sim-Set 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. 5). This is il-Table 3. Description of the parameters of the various sets of simulations performed to test the robustness of the width and orientation estimates recovery by FilDReaMS.

SimSet ID

Simulation parameters width orientation aspect ratio central intensity overlap noise level

[px] r b I 0 (σ W , σ B )% × I 0 SimSet 1 [5, 20] [0 • , 90 • ] a 3 1 no (5, 30) SimSet 2 [5, 20] [0 • , 90 • ] a 3 1 no (10, 100) SimSet 3 [5, 20] [0 • , 90 • ] a 10 1 no (5, 30) SimSet 4 9; 15 [0 • , 90 • ] b [3, 10] c [0.2, 5] d yes (5, 30)
Notes. The 'bar width' (in pixel), 'orientation', 'aspect ratio' (r b ), and 'central intensity' (I 0 ) provide the ranges of those parameters used to build the simulated bars, which are specified to present 'overlap' or not ; the 'noise level' is described with two components, a level of white noise (σ W ) and Brownian noise (σ B ). a All values are considered within the interval with a step of 1 • . b Randomly picked in the considered range.

c Randomly assigned to one out 8 possible values in the interval, following linear steps of 1 px.

d Randomly assigned to one out 8 possible values in the interval, distributed in logarithmic scale. lustrated 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 W ⋆ b 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 0 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 = 2R 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 W ⋆ b 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 . 5). The increase of noise level increases the frequency and the strength of the spurious peaks. The observed shift between the widths 2R flat and W ⋆,peak b can be empirically predicted, as we will discuss in Sect.4.2.3.

Impact of the filament aspect ratio

For this analysis we use the two sets of simulations, SimSet 1 & 3, between which the aspect ratio has been changed from 3 to 10. The peak becomes increasingly dominant as r b increases, leading to a gain of a factor of 2 to 5 of the precision, for Plummer and ideal cases respectively (see Table. 5). 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 2R flat and W ⋆,peak b 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.

The error in the fit for SimSet 1 can be roughly estimated at ≲ 1 px 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 . 5). 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 recover the width of input Plummer-type filaments in the case of high noise level.

For completeness, we also look at the impact of the Plummer index on the curves W ⋆peak b versus 2R 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 W ⋆peak b versus 2R flat are again shown in solid lines, and their parameters given in Table 4. 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 W ⋆ b . As for

Overlapping filaments of different properties

To study reliability of 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 W ⋆ b , averaged over the 100 maps, has two pronounced peaks at W = 15 px (bottom panel of Fig. 9). Hence, FilDReaMS recovers again the widths of the input ideal filaments. The reason why the histogram also contains non-negligible values away from the peaks is 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 W ⋆ b , averaged over the 100 test maps, has three peaks at W ⋆peak b = 9, 11 and 19 px which present a larger spread compared to the ideal case (bottom panel of Fig. 10). The two peaks at W ⋆peak b = 9 and 11 px are associated with the same input filament width (R flat = 9 px). This can be explained by the log-scale distribution of the central intensity of the filaments, I 0 , which leads to two different noise regimes, with 1/3 of the filaments in the high noise level, and 2/3 in the default noise level. Following top panel of Fig. 8, a 2R 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 behaviour is not observed in this simulation for the population of larger filaments, since they are more affected by the impact of overlapping filaments, leading to a spread distribution masking this effect.

The corresponding values of 2R flat can be derived from the linear fits W 4. We use the case p = 2.2 and the high noise level for the first peak at 9 rmpx (as explained above), and the case p = 2.2 with default noise level for the others. This leads to 2R 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 2R flat and W ⋆peak b 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.

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, ψ ⋆ f , are compared to the input orientations, ψ f over the 2000 noise realisations of each set of simulation. This allows to estimate both the bias and the standard deviation σ ψ of ψ ⋆ f . We first observe that the orientations of the reconstructed filaments don't show any deviation in average from the input orientation angles, and don't exhibit any specific dependency with the input value of the orientation angle. We also observe from these simulations that σ ψ decreases with the filament width, so 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 5.

We can see that the precision is very good, ≤ 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 (SimSet 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 the one observed for the width, with a degradation of only a factor of ≤ 2 for both kind 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.

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 sub-structures of larger ones, while only a small fraction of them covers areas not already included in larger structures, stressing the hierarchical organisation 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 by-eye inspection 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 than for 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, with however 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.

Concluding remarks

In this paper, we presented a new method to detect filaments of different widths in a map of intensity (in a broad sense) I.

A&A proofs: manuscript no. output Table 5. Uncertainties of FilDReaMS obtained in three different cases (left) in the estimation of width (center) and orientation (right). We called our method FilDReaMS, for Filament Detection and Reconstruction at Multiple Scales. In brief, FilDReaMS uses a rectangular model bar of width W b and aspect ratio r b , where W b 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, (ψ ⋆ f ) i , 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, (W ⋆ b ) i , is derived for each pixel i of all the reconstructed filaments, and the most prevalent bar widths of the map, W can then be cautiously converted to the widths of the often-used Plummer-type profiles, 2R flat (see Fig. 8 and Table 4).

SimSet ID

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 -30 sec for each value of W b and roughly 10 -20 min to cover the entire range of W b ; the user-friendliness, which makes it particularly suited for statistical studies.

FilDReaMS opens wide perspectives of application to astrophysical data, in order to study the processes of the stellar formation inside the interstellar filamentary structures. To start with we investigate the interplay between the orientation of the filaments and the Galactic magnetic field in a companion paper (Carrière et al. 2022), generalising the analysis of Malinen et al. (2016) to four other fields. A broader statistical analysis over 116 Herschel fields is also in preparation.

A

  Initial map of intensity I I Intensity in a broad sense (intensity, column density, or temperature) B Smoothed map resulting from the convolution of A with a 2D top-hat kernel of radius R R Radius of the 2D top-hat kernel C Binary map derived from A-B i Index of the considered pixel of map C for the detection of bar-like filaments ∆ψ Width of the angular window over which f M is compared to f I χ M r Measure of the normalized difference between f M and f I (Eq. 2) I 0 Central intensity of synthetic filaments j Index of the considered pixel of map A in one Monte-Carlo iteration σ A j Standard deviation of sub-region A j of map A, centered on pixel j (S/N) fil = I 0 /σ A j Signal-to-noise ratio of the ideal filament in the Monte-Carlo simulations the most significant filament σ W Standard deviation of white noise in the sets of simulations σ B Standard deviation of Brownian noise in the sets of simulations N pix Number of pixels whose most significant bar width is W ⋆ b N map Total number of pixels in the map W ⋆peak b

  Fig. 1. Illustration of the FilDReaMS method applied to one pixel i of the H 2 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: Sub-region C i of map C centered on pixel i, with the model bar (width W b , length L b , and orientation angle ψ b ) overplotted in green. Bottom Middle: Measured orientation function, equal to the fraction of the bar area covered by yellow pixels, f M i , as a function of ψ b (blue curve). Also shown are the ideal orientation function, f I i (orange curve; see Sect. 3.4) and the angular window of width ∆ψ (Eq. 1). Bottom Left: Sub-region A i 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 f M i ).

Fig. 2 .

 2 Fig. 2. Illustration of a case where the considered pixel i lies at the intersection of two filaments with different orientations. Top: Sub-region 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.

Fig. 3 .

 3 Fig. 3. Illustration of one ideal filament and its ideal orientation function. Left: Map of an ideal filament (filament having 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, f I i , from the second step of FilDReaMS applied to the binary map.

Fig. 4 .

 4 Fig. 4. Illustration of our Monte-Carlo simulations for the computation of the threshold (χ r ) th for a given value of the bar width, W b . Leftmost: Initial map A of the G210 field, with a sub-region A j centered on an arbitrary pixel j. Middle left: Synthetic map, A S j , composed of sub-region A j plus an ideal filament centered on pixel j. Middle right: Binary map, C S j , from the first step of FilDReaMS applied to A S j , with the model bar overplotted in green. Rightmost: Synthetic orientation function, f S j (blue curve), from the second step of FilDReaMS applied to C S j , together with the ideal orientation function, f I j (orange curve). Bottom: Monte-Carlo distribution D of the parameter (χ S r ) j that compares f S and f I 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 %.

Fig. 5 .

 5 Fig. 5. Maps of the significance, S (Eq. 6), of the G210 Herschel field for a bar width W b = 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).

Fig. 6 .

 6 Fig. 6. Illustration of the method used to reconstruct physical filaments of bar width (W b = 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 W b 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.

  FilDReaMS model bar Width of the detected filaments r b = L b /W b Aspect ratio of the FilDReaMS model bar Minimum elongation of the detected filaments (S/N) fil Signal-to-noise ratio of the ideal filament in the Monte-Carlo simulations Minimum contrast of the detected filaments detected with FilDReaMS.

  , 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). Both types of noise have different realisations 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 0 and σ B = 0.3 I 0 , while in the second case we chose σ W = 0.1 I 0 and σ B = I 0 .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, it consists of the coaddition of a map of synthetic filaments with realisations of realistic noise. Each original map is made of 91 non-overlapping 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 are thus created for each value of the width, W b ranging from 5 to 20 px (in steps of 1 px), on which 2000 realisations of realistic noise are co-added.For the last set of simulations, SimSet 4, it consists of a single series of 100 maps containing 64 synthetic filaments: 48 filaments with W b = 9 px 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 out of 8 possible values of the filament aspect ratio, r b , in the linear range[3, 10] and (independently) one out of 8 possible values of the filament intensity, I 0 , in the logarithmic range [0.2, 5]. Each value of r b and each value of I 0 are assigned to exactly eight filaments. The filaments are located at random positions, they have random orientations, and they can possibly overlap. For each map realisation, we co-add a noise realisation using the default configuration.

Fig. 7 .

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

  4.2.3. Plummer-to-ideal width relationThe sets of simulations with Plummer-type filaments allow us to derive the empirical relation between 2R flat and W ⋆peak b for different noise levels and also for different values of the Plummer index, p. The curves linking the input 2R flat to the output W ⋆peak b 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 W ⋆peak b = 2R flat (dashed line), the linear fits to the two curves W ⋆peak b versus 2R flat (solid lines in the corresponding colors), as well as the standard deviations computed on each average histogram (N pix /N map ) versus W ⋆ b (error bars on each symbol in the corresponding colors). For future reference, the parameters a and b of the linear fits W ⋆peak b = a (2R flat ) + b are given in Table

Table 4 .Fig. 8 .

 48 Fig. 8. Average relation between the width 2R flat of an input Plummertype filament and the most prevalent bar width W ⋆peak b 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 W ⋆peak b versus 2R flat of the same colors, and the dashed line is the straight line W ⋆peak b

Fig. 9 .

 9 Fig. 9. Visualisation of the set of simulations SimSet 4 for ideal filaments, showing how well FilDReaMS recovers the widths W b of overlapping filaments with different widths, aspect ratios, and intensities. Top left: One of the 100 test maps, containing 64 ideal filaments, with width W b = 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, N pix , whose most significant bar width is W ⋆ b , normalized to the number of pixels in the map, N map , as a function of W ⋆ b .

Fig. 10 .

 10 Fig. 10. Same as Fig. 9, with the input ideal filaments of width W b replaced by Plummer-type filaments of width 2R flat .

  ⋆peak b = a (2R flat ) + b, with the values of a and b given in Table

  from the peaks of the histogram of W ⋆ b . The W ⋆peak b

Fig. 11 .

 11 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 W b = 6 px (left), W b = 10 px (middle), and W b = 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 N H 2 .

  • , +90 • ], possibly by adding or subtracting 180 • .

Table 1 .

 1 List of all the symbols used in the paper.

Here, we use the term striations to refer to the faint and periodic structures seen in the Herschel maps. These are similar to the periodic magnetic-field-aligned structures detected in the diffuse 12 CO emission from the Taurus molecular cloud(Goldsmith et al.

2008;Narayanan et al. 2008).

For each value of the bar width, W b , the value of R is calculated once and for all for the entire map A (see Sect.

3.1); it is not recalculated at each Monte-Carlo iteration.

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