Issue 
A&A
Volume 668, December 2022



Article Number  A41  
Number of page(s)  15  
Section  Interstellar and circumstellar matter  
DOI  https://doi.org/10.1051/00046361/202243506  
Published online  02 December 2022 
FilDReaMS
I. Presentation of a new method for Filament Detection and Reconstruction at Multiple Scales
IRAP, Université de Toulouse, CNRS,
9 avenue du Colonel Roche,
BP 44346,
31028
Toulouse Cedex 4, France
email: jeansebastienpaulcarriere@gmail.com
Received:
9
March
2022
Accepted:
13
September
2022
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, userfriendly, multiscale, 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, and up to 3° in the worstcase scenario with high levels of noise. The width of the filament is recovered with uncertainties of better than 0.5 px (pixels) in most cases, which could extend up to 3px in the case of low signaltonoise ratios. Some attempt to build a correspondence between Plummertype filament profiles and the outcomes of the method is proposed, but remains sensitive to the local environment.
Conclusions. We find our FilDReaMS to be robust and adapted to the identification and 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.
Key words: ISM: clouds / ISM: structure / ISM: magnetic fields / infrared: ISM / submillimeter: ISM / techniques: image processing
© J.S. Carrière et al. 2022
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the SubscribetoOpen model. Subscribe to A&A to support open access publication.
1 Introduction
A large variety of methods have already been developed to extract elongated structures from twodimensional (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 nonlocal 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 multiscale decomposition.
Local methods compute either the gradient (firstorder derivatives; e.g., Soler et al. 2013; Planck Collaboration XXXV 2016b) or the Hessian matrix (secondorder 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 columndensity) 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 ^{13}CO intensity maps (Panopoulou et al. 2014). A limitation of the local approach is the difficulty of detecting faint structures such as striations^{1}.
Nonlocal methods focus on a given scale around each pixel. The templatematching 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 multiscale 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 finetuning and additional tools to extract the filament orientations and scales. It was already applied to Herschel maps (Cox et al. 2016; RiveraIngraham et al. 2016, 2017; Arzoumanian et al. 2019). The waveletbased methods by Robitaille et al. (2019) and by OssenkopfOkada & Stepanov (2019) use an anisotropicwavelet 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 wellknown biases of other typically used methods (Panopoulou et al. 2017) and remains relatively fast. However, waveletbased methods require additional steps to extract filament orientations. Although these methods are multiscale, 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 templatematching 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 filamentextraction 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 multiscale 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) keyprogram (Juvela et al. 2010, 2012), the socalled G210 field, which corresponds to the highGalacticlatitude starforming 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_{2} 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 onethird of the beam size, is 12″, corresponding to 0.0081 pc. The H_{2} column density, , varies over the range [0.2,7.5] × 10^{21} cm^{−2}.
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 largescale structures, a binary version of the original image is first built by subtracting the smoothed image with a tophat 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.
List of all the symbols used in the paper.
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 tophat 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 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 tophat 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 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, , 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 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 ). 
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_{2} 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 tophat 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 tophat kernel to produce a smoothed map B (topmiddle 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 tophat 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 bottommiddle 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 barlike 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 pixelscale fluctuations (see bottomright 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:
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 ψ_{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} (bottomright 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 (bottomleft 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 barlike 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
where
Ν_{ψ} 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 chisquared, 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 MonteCarlo 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
or, equivalently,
where
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 MonteCarlo simulations
The purpose of the MonteCarlo 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 tophat function of radius R will be possible (as explained in Sect. 3.1). A synthetic map is then created (middleleft panel) by superposing an ideal filament with uniform intensity I_{0} on subregion A_{j}, centered on pixel j and oriented at random (over a uniform distribution of angles). For convenience, I_{0} is written in terms of the standard deviation of subregion A_{j}, , and the signaltonoise 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 tophat function of radius R^{2}) 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 MonteCarlo 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 MonteCarlo 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 MonteCarlo 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 subregion A_{j} centered on an arbitrary pixel j. Middle left: synthetic map, , composed of subregion A_{j} 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: MonteCarlo distribution D of the parameter 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 pvalue 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 barlike 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 barlike 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 barlike 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 barlike filaments are found (bottomleft 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 barlike 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 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 ). 
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 signaltonoise ratio of the ideal filament in the MonteCarlo 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 signaltonoise ratio of the ideal filament in the MonteCarlo 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 smallscale 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 onethird 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_{0}.
Plummertype filaments: their intensity can be described by a 2D Plummertype profile with halfwidth R_{flat} and aspect ratio rb:
where x and y are the coordinates across and along the long axis of the filament, p is the Plummer powerlaw 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 powerlaw (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 powerlaw index of −2) is meant to represent astrophysical fluctuations (e.g., MivilleDeschê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_{0} and σ_{Β} = 0.3 I_{0}, while in the second case we chose σ_{W} = 0.1 I_{0} and σ_{Β} = 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, the configuation consists of the coaddition 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 coadded.
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_{0}, in the logarithmic range [0.2, 5]. Each value of r_{b} and each value of I_{0} 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 coadd a noise realization using the default configuration.
All sets of simulations are finally repeated for ideal and Plummertype 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 Plummertype filaments with powerlaw 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 (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 barlike 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 barlike 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. 
Free parameters of FilDReaMS, with their definitions (second column) and the corresponding characteristics of the filaments to be detected by FilDReaMS (third column).
Description of the parameters of the various sets of simulations performed to test the robustness of the width and orientation estimates recovery by FilDReaMS.
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 welldefined 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 Plummertype 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 powerlaw 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 Plummertype filament. The integral of the transverse profile of a Plummertype 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 powerlaw 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 Plummertype filament with default noise, we still obtain a precision of 0.4 px. However, the histograms (N_{pix}/N_{map}) versus of Plummertype 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 (N_{pix}/N_{map}) versus obtained for ideal (blue) and Plummertype (orange, with powerlaw index p = 2.2) filaments with bar width W_{b} = 12 px and Plummer halfwidth R_{flat} = 12 px, respectively, in three different sets of simulations, SimSet 1 (top), SimSet 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 Plummertype 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 Plummertoideal width relation
The sets of simulations with Plummertype 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 Plummertype 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 spreadout Plummertype 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 R_{flat} of an input Plummertype 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 R_{flat} of the same colors, and the dashed line is the straight line . 
Uncertainties of FilDReaMS obtained in three different cases (left) in the estimation of width (center) and orientation (right).
Fig. 9 Visualization 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 , normalized to the number of pixels in the map, N_{map}, 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 Plummertype 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 nonnegligible 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 Plummertype 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 logscale distribution of the central intensity of the filaments, I_{0}, which leads to two different noise regimes, with onethird of the filaments in the high noise level, and twothirds 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 W_{b} replaced by Plummertype filaments of width 2 R_{flat}. 
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 Plummertype 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 Plummertype 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 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 (SimSet 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 W_{b} = 6px (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 tophat 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 tophat 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 northerncentral 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 oftenused Plummertype 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 userfriendliness, which makes it particularly suited for statistical studies.
FilDReaMS opens wide perspectives of application to astrophysical 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, UCLMSSL, 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
 Alina, D., Ristorcelli, I., Montier, L., et al. 2019, MNRAS, 485, 2825 [Google Scholar]
 Arzoumanian, D., André, P., Didelon, P., et al. 2011, A&A, 529, L6 [Google Scholar]
 Arzoumanian, D., André, P., Könyves, V., et al. 2019, A&A, 621, A42 [Google Scholar]
 Carrière, J.S., Ferrière, K., Ristorcelli, I., & Montier, L. 2022, A&A, 668, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Clark, S. E., & Hensley, B. S. 2019, ApJ, 887, 136 [Google Scholar]
 Clark, S. E., Peek, J. E. G., & Putman, M. E. 2014, ApJ, 789, 82 [Google Scholar]
 Clark, S. E., Hill, J. C., Peek, J. E. G., Putman, M. E., & Babler, B. L. 2015, Phys. Rev. Lett., 115, 241302 [Google Scholar]
 Cox, N. L. J., Arzoumanian, D., André, P., et al. 2016, A&A, 590, A110 [Google Scholar]
 Goldsmith, P. F., Heyer, M., Narayanan, G., et al. 2008, ApJ, 680, 428 [Google Scholar]
 Juvela, M. 2016, A&A, 593, A58 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Juvela, M., Ristorcelli, I., Montier, L. A., et al. 2010, A&A, 518, A93 [Google Scholar]
 Juvela, M., Ristorcelli, I., Pagani, L., et al. 2012, A&A, 541, A12 [Google Scholar]
 Koch, E. W., & Rosolowsky, E. W. 2015, MNRAS, 452, 3435 [Google Scholar]
 Malinen, J., Montier, L., Montillaud, J., et al. 2016, MNRAS, 460, 1934 [Google Scholar]
 Men’shchikov, A. 2013, A&A, 560, A63 [Google Scholar]
 Micelotta, E. R., Juvela, M., Padoan, P., et al. 2021, A&A, 647, A121 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 MivilleDeschênes, M. A., Levrier, F., & Falgarone, E. 2003, ApJ, 593, 831 [Google Scholar]
 Montillaud, J., Juvela, M., RiveraIngraham, A., et al. 2015, A&A, 584, A92 [Google Scholar]
 Narayanan, G., Heyer, M. H., Brunt, C., et al. 2008, ApJS, 177, 341 [NASA ADS] [CrossRef] [Google Scholar]
 OssenkopfOkada, V., & Stepanov, R. 2019, A&A, 621, A5 [Google Scholar]
 Ostriker, J. 1964, ApJ, 140, 1056 [Google Scholar]
 Palmeirim, P., André, P., Kirk, J., et al. 2013, A&A, 550, A38 [Google Scholar]
 Panopoulou, G. V., Tassis, K., Goldsmith, P. F., & Heyer, M. H. 2014, MNRAS, 444, 2507 [Google Scholar]
 Panopoulou, G. V., Psaradaki, I., & Tassis, K. 2016, MNRAS, 462, 1517 [Google Scholar]
 Panopoulou, G. V., Psaradaki, I., Skalidis, R., Tassis, K., & Andrews, J. J. 2017, MNRAS, 466, 2529 [Google Scholar]
 Peretto, N., André, P., Könyves, V., et al. 2012, A&A, 541, A63 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXXII. 2016a, A&A, 586, A135 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXXV. 2016b, A&A, 586, A138 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Polychroni, D., Schisano, E., Elia, D., et al. 2013, ApJ, 777, L33 [Google Scholar]
 RiveraIngraham, A., Ristorcelli, I., Juvela, M., et al. 2016, A&A, 591, A90 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 RiveraIngraham, A., Ristorcelli, I., Juvela, M., et al. 2017, A&A, 601, A94 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Robitaille, J. F., Motte, F., Schneider, N., Elia, D., & Bontemps, S. 2019, A&A, 628, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schisano, E., Rygl, K. L. J., Molinari, S., et al. 2014, ApJ, 791, 27 [Google Scholar]
 Soler, J. D., Hennebelle, P., Martin, P. G., et al. 2013, ApJ, 774, 128 [Google Scholar]
 Sousbie, T. 2011, MNRAS, 414, 350 [NASA ADS] [CrossRef] [Google Scholar]
 Suri, S., SánchezMonge, Á., Schilke, P., et al. 2019, A&A, 623, A142 [Google Scholar]
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 magneticfieldaligned structures detected in the diffuse ^{12}CO emission from the Taurus molecular cloud (Goldsmith et al. 2008; Narayanan et al. 2008).
All Tables
Free parameters of FilDReaMS, with their definitions (second column) and the corresponding characteristics of the filaments to be detected by FilDReaMS (third column).
Description of the parameters of the various sets of simulations performed to test the robustness of the width and orientation estimates recovery by FilDReaMS.
Uncertainties of FilDReaMS obtained in three different cases (left) in the estimation of width (center) and orientation (right).
All Figures
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 tophat 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 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, , 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 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 ). 

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 MonteCarlo 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 subregion A_{j} centered on an arbitrary pixel j. Middle left: synthetic map, , composed of subregion A_{j} 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: MonteCarlo distribution D of the parameter 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 pvalue of 5%. 

In the text 
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 ). 

In the text 
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 barlike 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 barlike 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 (N_{pix}/N_{map}) versus obtained for ideal (blue) and Plummertype (orange, with powerlaw index p = 2.2) filaments with bar width W_{b} = 12 px and Plummer halfwidth R_{flat} = 12 px, respectively, in three different sets of simulations, SimSet 1 (top), SimSet 2 (middle) and SimSet 3 (bottom). 

In the text 
Fig. 8 Average relation between the width 2 R_{flat} of an input Plummertype 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 R_{flat} 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 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 , normalized to the number of pixels in the map, N_{map}, as a function of . 

In the text 
Fig. 10 Same as Fig. 9, with the input ideal filaments of width W_{b} replaced by Plummertype filaments of width 2 R_{flat}. 

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 W_{b} = 6px (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 tophat radius R = 11 px (Malinen et al. 2016). The color code of the reconstructed filaments indicates their column density . 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.