Measuring Filament Orientation: A New Quantitative, Local Approach

Elongated filaments of gas and dust are ubiquitous in molecular clouds (e.g., Molinari et al. 2010). These clouds are stellar nurseries and the filaments they host may play an important role in star formation, with the majority of star-forming cores lying along filaments "like beads on a string" (André et al. 2014).

Filaments represent velocity-coherent overdensities of gas and dust, and have aspect ratios greater than at least three (Panopoulou et al. 2014). They can be identified from a 2D astronomical image such as, for instance, a column density map, using skeleton-based filament identification algorithms such as filfinder (Koch & Rosolowsky 2015). For a given input image and set of input parameters, these return a filament skeleton. The skeleton is a one-pixel-wide representation of the filamentary structure in the original image, tracing the main path of the filament and its branches. Clumps and cores are also overdense compared to their surroundings, and are distinguished from filaments by smaller aspect ratios of ∼2 (Tachihara et al. 2000). Clumps are inhomogeneously dense velocity-coherent regions from which a system of stars may form. A core is a dense velocity-coherent region that may form a single star or binary star. Cores are usually found grouped into clumps.

In the study of these filaments, one useful measurement is that of their orientation. Filament orientation is used in the construction of radial profiles used to derive filament width. Filament orientation can also be compared with that of the magnetic field. Magnetic fields are believed to have a dynamically important role in filament formation and stability. In several theories of cloud structure formation, matter is channeled along the field lines, allowing filaments to form through gravitational contraction (Nakamura & Li 2008). In this scenario, dense filaments are aligned perpendicular to the field and less dense filaments are aligned parallel (Li et al. 2008). Goldsmith et al. (2008) and Planck Collaboration et al. (2016) find observational evidence for this scenario.

Filament orientation can be measured from the filament skeleton (Koch & Rosolowsky 2015). The intensity changes at the edge of the skeleton, and this intensity gradient has an associated direction (see Figure 1). Here we propose a new method to derive the filament orientation, exploiting this fact. This is achieved through the use of the Sobel filter, described in Section 2, and some additional post-processing steps discussed in Section 3. This method, which we call the "Sobel-gradient method," returns a quantitative and local map of filament orientation for any filament skeleton, including those with complex interconnected structures.⁶ The map reveals how the orientation changes as the skeleton curves on a local scale. We explore the uncertainties associated with the method in Section 4. Applications for this method are suggested in Section 5.

Zoom In Zoom Out Reset image size

1.1. Motivation

There are two main existing quantitative approaches to measuring filament orientation, the first being a map-based analysis (e.g., the Schisano et al. 2014 Hessian matrix method), and the second being a skeleton-based analysis (e.g., the Koch & Rosolowsky 2015 filfinder algorithm). Prior to their introduction, the predominantly utilized method for measuring filament and field relative orientation was a qualitative, global, visual comparison (e.g., Goldsmith et al. 2008; Busquet et al. 2013; Palmeirim et al. 2013), and this approach is still used in more recent works (e.g., Kusune et al. 2016).

The Schisano et al. (2014) method uses the Hessian matrix to identify filaments and measure their orientation from a 2D astronomical map such as a Herschel⁷ (Pilbratt et al. 2010) dust column density map. However, Schisano et al. (2014, p. 4) state that "for very complex features where filaments are organized in web-like structures, the cross-spine profile fitting often fails to converge." In our related project studying the filamentary structure of the South-ridge (SR) and Center-ridge (CR) of the Vela C Molecular Cloud (C.-E. Green et al. 2017, in preparation), we tried this method to derive the orientation of the filaments shown in Figure 2 panel (a). Indeed, the method failed to converge for this complex interconnecting data set, motivating our search for an alternate method.

Zoom In Zoom Out Reset image size

The Koch & Rosolowsky (2015) filament identification algorithm filfinder has an inbuilt skeleton-based filament orientation calculator that uses a "line-based" approach, the Rolling-Hough Transform (RHT). In this type of approach, test lines with different angles are fit to groups of pixels along the filament skeleton, finding the best-fit line and thus the associated angle of the skeleton segments. Koch & Rosolowsky (2015) define filament orientation to be the weighted directional mean of the distribution of angles for the skeleton returned by the RHT. filfinder thus returns to the user a single orientation value over long filament segments (e.g., ∼40 pixels, which corresponds to 1.2 pc in our example Vela C SR data in Figure 2 panel (a)), whereas our more complex filaments curve on a smaller scale of ∼5 pixels (0.15 pc). This definition of filament orientation is therefore not compatible with our goal of a quantitative, local, "position-by-position" filament orientation.

This motivated our search for an alternate, fully automated filament orientation measurement method that:

• 1.  
  returns a quantitative, local, "position-by-position" measurement of filament orientation,
• 2.  
  can be applied to complex interconnecting filaments (such as the SR shown in Figure 2 panel (a)), as well as simpler, more linear filaments.

As previous map-based approaches such as Schisano et al.'s (2014) Hessian matrix method have not worked for these more complicated "looped" (in the 2D image) filaments we focused our search on a skeleton-based approach that could provide a measurement on a smaller scale than the RHT method built into filfinder. This led us to develop the new Sobel-gradient method that we propose here, exploiting the image intensity gradient to arrive at a map of filament orientation.

Filament skeletons are generally output by filament identification algorithms in the Flexible Image Transport System (FITS⁸ ) file format where the pixels "on" the skeleton have a value of one, and those "off" the skeleton have a value of zero. These can be trivially converted to a grayscale image matrix, where the "on" skeleton pixels are white, with a value of 255 and the "off" skeleton pixels are black, with a value of zero. We use the python scipy ndimage implementation of the Sobel filter where it is necessary to use this convention. White could also be represented as a value of one, black as a value of zero, and gray shades as decimal values in between, if a different implementation were used. As the skeleton is a binary image, the image intensity gradient only exists at the edges of the filament skeleton where the intensity changes. This intensity gradient has a magnitude and a direction. An example of the image intensity gradient vector for an idealized case is shown in Figure 1. With some minor adjustments, the skeleton orientation can be derived from the intensity gradient direction.

2.1. Intensity Gradient Direction

To calculate the direction of the intensity gradient, we need the first x and y derivatives (G_x and G_y, respectively) of the skeleton image matrix, I. The direction, Θ, of the gradient is calculated as

$$\begin{eqnarray}&&{\rm{\Theta }}={\tan }^{-1}({G}_{y}/{G}_{x}).\end{eqnarray} \tag{ 1 }$$

The Sobel filter is commonly used in computer vision to estimate these derivatives (Gao et al. 2010). It is already built into matlab and the python scipy ndimage library, so this method can be quickly and easily implemented. The Sobel filter itself is computationally inexpensive. Its speed is a major advantage because orientation measurements are generally repeated for the multiple different skeletons produced by different combinations of input parameters to the filament identification algorithm. In some ways this approach is similar to the Histogram of Relative Orientations (HRO) method of Soler et al. (2013), which also uses Gaussian derivatives (of which the Sobel filter is one of the simplest types) to measure the orientation of molecular cloud structure. Our approach differs in that we aim to find the orientation only of strictly defined and identified filaments, by measuring the orientation of the one-pixel-wide filament skeleton. This is in contrast to the HRO method, which makes no structure definitions, and involves finding the orientation of all structures of all scales within a column density map.

2.1.1. The Sobel Filter

The Sobel filter is a discrete differential operator consisting of two 3 × 3 matrices of coefficients. When convolved with the image matrix, I, two new image matrices are created, representing estimates of G_x and G_y as follows (where ∗ represents the convolution operation; Gao et al. 2010):

$$\begin{eqnarray}{G}_{x}=\left[\begin{array}{ccc}-1 & 0 & +1\\ -2 & 0 & +2\\ -1 & 0 & +1\end{array}\right]\ast I,\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}{G}_{y}=\left[\begin{array}{ccc}+1 & +2 & +1\\ 0 & 0 & 0\\ -1 & -2 & -1\end{array}\right]\ast I.\end{eqnarray} \tag{ 3 }$$

Angular measurements have a reference point and a direction of increase. For the image gradient returned by the Sobel filter, these are the horizontal and the anticlockwise directions. Therefore, in that convention the gradient angle needs to be rotated by 90° to give the angle of the edge. However in the astronomical convention the reference point is north (often vertically upward in astronomical images), with an anticlockwise direction of increase. We are operating in the domain of [90, −90], so we therefore only need to perform a simple sign reversal to arrive at an estimate of the orientation of the skeleton edge⁹ when working in the astronomical convention.

Throughout this work we will use, for the purpose of illustration, filament skeletons identified with filfinder from the Fissel et al. (2016) Herschel dust column density images of the SR and CR of Vela C. For the SR these are shown in Figure 2 panels (a) and (b), respectively. In the SR and CR data one pixel corresponds to 0.03 pc. The skeletons were selected as belonging to the group of optimum skeletons, most similar to, and therefore the best representation of, the original column density image, using the mean structural similarity index as a goodness-of-fit measure, as described in Green et al. (2017). Together the selected skeletons from the SR and CR comprise a representative data set, containing interconnected "loops" and curvature on a small scale, along with some more linear segments. They were selected because they contained the largest number of "difficult" features for the algorithm to tackle. The SR skeleton selected was produced by filfinder input parameters of: a skeleton threshold (skeleton length cutoff) of 10 pixels (0.3 pc, corresponding to an aspect ratio of 3 (Panopoulou et al. 2014), given an assumed width of 0.1 pc (Arzoumanian et al. 2011)), a branch threshold (branch length cutoff) of 3 pixels (0.09 pc), a global threshold (noise threshold) of 69%, and a flattening threshold (threshold for arctan flattening, which removes the impact of compact sources like clumps in the masking step) of 60%. The CR skeleton was produced by a skeleton threshold of 10 pixels (0.3 pc), a branch threshold of 5 pixels (0.15 pc), a global threshold of 74%, and a flattening threshold of 96%.

Before deriving the filament skeleton orientation we first automatically remove junction points¹⁰ from the skeleton. Junction points belong to all of the intersecting filaments involved and therefore have an undefined orientation. Removing these breaks up the skeleton into many components, as shown in Figure 2 panel (c), which in computer vision are called connected components. The image gradient is defined at the pixels immediately surrounding the skeleton, therefore the gradients may overlap and overwrite each other at the new endpoints created by deleting the junctions since they are so close together. To avoid this issue we automatically locate and label the connected components¹¹ and repeat the Sobel-gradient method described in the following on each component separately, collecting the final orientation maps of each component into a "master map," the final filament skeleton orientation map.

To derive the skeleton orientation of each component we calculate the x and y image derivatives using the Sobel filter (Equations (2) and (3), respectively), and then calculate the skeleton image gradient using Equation (1). The sign of the skeleton image gradient is then reversed for consistency with astronomical conventions. This process is illustrated in Figure 2 panels (d)–(f).

The sign-reversed gradient direction map in Figure 2 panel (f) gives the orientation of the edge of the skeleton. To move from this to a map of the orientation of the skeleton's path we perform some simple post-processing steps, which are demonstrated in Figure 3. These correct minor, partially cosmetic issues that arise as a direct consequence of the nature of the Sobel filter image gradient approach. In these steps we (1) correct the branch ends, (2) infill the center pixels, (3) smooth the map, and (4) select the orientation values at the positions along the original skeleton to save into the master map. These steps are illustrated in Figure 3.

Zoom In Zoom Out Reset image size

First, at the ends of the connected components there are pixels with angle values that are roughly orthogonal to the rest of the component. This is because of the additional exposed pixel edges at the component ends. We automatically detect and delete¹² the handful of affected pixels at component ends in the map. We then give them the value of the circular vector average of the closest¹³ unaffected values in the branch, resolving the issue as shown in Figure 3 panel (b). Wherever angles undergo averaging we use circular vector averaging. This ensures angles are averaged correctly, accounting for the fact that they are a circular quantity that wraps back around such that, e.g., 90° and −90° are equal and represent the horizontal in the astronomical definition if north is vertical.

Second, the intensity gradient only exists at the edge of the skeleton, thus leaving a partially blank center.¹⁴ We automatically infill the pixels along the positions of the original skeleton with the circular vector average of their neighboring pixels within a 5 × 5 pixel window (also including the center pixel in the average if it is not blank). The result is shown in Figure 3 panel (c).

Finally, we smooth the map. A 5 × 5 pixel window is passed over the image and we consider only pixel positions lying on the original skeleton that we center in the window. The circular vector average of the pixels inside the window is calculated and is saved to the corresponding pixel position of the central pixel in a new map. This results in a smoothed orientation map for that connected component. The smoothed map for each connected component is saved into the master map. After repeating the process for each component, we arrive at a map of the filament skeleton's orientation as illustrated in Figure 3 panel (d).

These three post-processing steps take us from the map of the gradient direction to a quantitative skeleton orientation map that reveals how the orientation changes as the skeleton curves on a local scale. For the first time, we present the filament orientations for the SR and CR of Vela C on this small scale in Figure 4.

Zoom In Zoom Out Reset image size

We validate the Sobel-gradient method against the known analytic case of the circle. We generate circles with radii of 7–500 pixels and calculate the theoretical orientation of their tangents at each point around them. We then apply the Sobel-gradient method to them and compare the orientations at each point around the circle. Circles drawn digitally with radii smaller than 7 pixels are essentially squares with the middle pixel along each side pushed out by one position. We therefore only consider radii larger than this.

For each radius we calculate the difference between the theoretical tangent orientation and the Sobel-gradient orientation at each point around the circle, and then find the maximum, mean, and standard deviation of those differences. These are plotted against their corresponding radii in Figure 5. The average of the maximum differences at each radius was 78, the average of the mean was 21, and the average of the standard deviations was 16. In calculating these values we include only the circles with radii of 57 pixels or greater, as these have 360 pixels (and therefore 360 unique angles) around them. Circles digitally generated with radii less than this are still very effected by the "squaring effect" present at small radii. This means that their theoretical orientation deviates greatly from the orientation that is actually drawn. They are thus not a reliable or accurate point of validation.

Zoom In Zoom Out Reset image size

The maximum differences are dominated by digitization errors, which we estimate to be up to ∼5°. The average of the standard deviations of the differences between the theoretical and Sobel-gradient orientations for each radii thus provides a more appropriate estimate of the uncertainty associated with the Sobel-gradient method. Consequently, we estimate the uncertainty of the Sobel-gradient method to be ∼2° based on this circle test analysis.

Circles have no start or endpoints, have no branches, and do not wiggle back and forth, changing their direction as real filaments do. To further gauge the uncertainty associated with the Sobel-gradient method in a realistic scenario we compare the Sobel-gradient orientation maps of two filamentary regions, the Vela C SR and CR, to those measured manually.¹⁵

A difference map was calculated for each region between the Sobel-gradient and manually measured orientation maps, and these are shown in Figure 6. The histograms of the difference maps for both regions are shown in Figure 7. The majority of orientations differed by less than one degree. The maximum difference for the SR was 71, the mean difference was 19, and the standard deviation was 18, while that for the CR were 72, 12, and 13, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

When measuring orientation manually, the nonlinear skeletal path was essentially decomposed into small linear sections. The Sobel-gradient method does not decompose the path in this way; it instead has a smooth transition along the filament that better reflects the curvature of the filament's path. The human-defined section does not always align perfectly with the corresponding section in the Sobel map; sometimes they are shifted off each other by 1–2 pixels. The larger orientation differences in the difference map mostly occur at locations where these shifts exist. This indicates that the larger difference values in the distribution in Figure 7 are likely caused by this affect. This issue is unavoidable in constraining the uncertainty on the Sobel-gradient method in a realistic scenario—there is currently no other method on this scale besides manual measurement to provide a comparison.

The maximum of the difference maps of 7° is therefore a poor measure of the actual uncertainty of the Sobel-gradient method. It is more appropriate to use the standard deviation of the combined difference distribution of the SR and CR of ∼2° to estimate this uncertainty. This value is in agreement with that obtained from the circle test analysis. We therefore conservatively estimate the uncertainty of the Sobel-gradient method to be ∼2°. This is acceptable considering that, for the Vela C BLASTPol data of Fissel et al. (2016), the average uncertainty of the magnetic field maps is ∼2°, reaching up to ∼16° in places.

The Sobel-gradient method is slightly more accurate than manual measurement, but its strength is that it is significantly faster. To measure the orientation of one skeleton (such as those presented here) on this small scale manually takes most of a working day. When the additional time to then input the manual orientation measurement into a FITS file is taken into account, the manual process takes about one to two working days per skeleton. In studies of filaments, the filament and field relative orientation and radial profile measurements that involve filament orientation are often repeated for hundreds of skeletons (all corresponding to the same input image, but to different combinations of input parameters to the filament identification algorithm). Consequently, a fast and accurate orientation measurement method is essential. The Sobel-gradient method allows automation of the filament orientation measurement and is therefore crucial in this era of "big data."

The Sobel-gradient method described is a technique to derive the orientation of filaments from their skeletons. This measurement has a number of astrophysical applications. One of the most significant is its use for calculating the relative orientation between magnetic fields and filaments, which provides clues on the role of magnetic fields in the formation and stability of filaments. We perform this calculation and present the results for the filaments of Vela C in C.-E. Green et al. (2017, in preparation). There are a number of other potential applications of the method, including investigating relations between filament orientation and filament column density, mass, spatial width, or molecular line width.

We have described a fully automated method to derive the orientation of a filament skeleton from the direction of the image intensity gradient that is suitable for complex, "looping" filamentary structures. We call this the "Sobel-gradient method." It allows a local measurement of filament orientation that reflects the often rapid changes in orientation as a filament curves. This means that the filament orientation calculated from the intensity gradient can be directly compared to a map of the magnetic field, giving a quantitative, local measure of relative orientation as opposed to the qualitative, global and "by-eye" technique that is the current predominantly adopted method. It also has a number of other applications for investigating relationships involving filament orientation, such as that between filament orientation and column density. We have found this method to have a high degree of accuracy, with an uncertainty of ∼2°. This computer vision technique provides the significant advantage that it can be easily automated, saving a significant amount of time compared to manual measurement, which is imperative in this era of "big data." It also has broader applications and can be applied to any image containing lines or edges to find their orientation.

The authors would like to thank the referee Erik Rosolowsky and the anonymous statistical editor for their helpful comments, which improved this work. The authors are also grateful to Eric Koch for helpful discussions surrounding the filfinder filament skeletons of this work. Herschel is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA. C.-E.G. gratefully acknowledges the support of the Layne Beachley Foundation. L.M.F. and G.N. acknowledge support from NASA grant NNX13AE50G and from the Center for Interdisciplinary Exploration and Research in Astrophysics. L.M.F. was supported in part by an NSERC Postdoctoral Fellowship. L.M.F. is a Jansky Fellow of the National Radio Astronomy Observatory (NRAO). NRAO is a facility of the National Science Foundation (NSF) operated under cooperative agreement by Associated Universities, Inc. This research made use of NASA's Astrophysics Data System.

Software: This project utilized astropy (http://www.astropy.org) (Astropy Collaboration et al. 2013), aplpy (http://aplpy.github.com), filfinder (Koch & Rosolowsky 2015), Karma visualization tools (Gooch 1996), and scipy (http://www.scipy.org/).