Removing visual bias in filament identification: a new goodness-of-fit measure

TL;DR

This paper introduces the use of the mean structural similarity index (MSSIM) as a mathematical goodness-of-fit measure to objectively evaluate and select the best filament skeletons from image analysis algorithms, reducing visual bias and enabling automation.

Contribution

It proposes a novel application of MSSIM for systematic, unbiased comparison of filament skeletons to original images, applicable across different algorithms and facilitating automation.

Findings

MSSIM effectively identifies skeletons most similar to input images.

The method enables systematic parameter optimization for filament identification.

It reduces the need for subjective visual assessment.

Abstract

Different combinations of input parameters to filament identification algorithms, such as Disperse and FilFinder, produce numerous different output skeletons. The skeletons are a one pixel wide representation of the filamentary structure in the original input image. However, these output skeletons may not necessarily be a good representation of that structure. Furthermore, a given skeleton may not be as good a representation as another. Previously there has been no mathematical `goodness-of-fit' measure to compare output skeletons to the input image. Thus far this has been assessed visually, introducing visual bias. We propose the application of the mean structural similarity index (MSSIM) as a mathematical goodness-of-fit measure. We describe the use of the MSSIM to find the output skeletons most mathematically similar to the original input image (the optimum, or `best', skeletons) for a given algorithm, and independently of the algorithm. This measure makes possible systematic parameter studies, aimed at finding the subset of input parameter values returning optimum skeletons. It can also be applied to the output of non-skeleton based filament identification algorithms, such as the Hessian matrix method. The MSSIM removes the need to visually examine thousands of output skeletons, and eliminates the visual bias, subjectivity, and limited reproducibility inherent in that process, representing a major improvement on existing techniques. Importantly, it also allows further automation in the post-processing of output skeletons, which is crucial in this era of `big data'.

Equations

$$SSIM(a,b)$$

$$where$$

$$l(a,b)$$

$$c(a,b)$$

$$s(a,b)$$

$$SSIM(a,b)=\frac{(2{\mu }_a{\mu }_b+C_1)(2{\sigma }_{ab}+C_2)}{({\mu }_a^2+{\mu }_b^2+C_1)({\sigma }_a^2+{\sigma }_b^2+C_2)}$$

Full text

ISM: structure — methods: statistical — methods: data analysis — stars: formation — techniques: image processing

1 Introduction

The advent of the Herschel111Herschel is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA. Space Observatory (Pilbratt et al., 2010) has allowed detailed mapping of the structure of stellar nurseries, molecular clouds. Herschel observations have affirmed the filamentary nature of these clouds and suggested these structures may play an important role in star formation (e.g. Molinari et al., 2010; Arzoumanian et al., 2011; Hill et al., 2011; Palmeirim et al., 2013). This has motivated a boom in the study of these filamentary structures in infra-red dust and radio molecular line data (e.g. Cox et al., 2016; Busquet et al., 2013). At the foundation of the study of these thread-like structures lies our ability to automatically identify them in such data. The construction of a skeleton tracing the filamentary structures is facilitated by filament finding algorithms such as disperse (Sousbie, 2011) and filfinder (Koch & Rosolowsky, 2015). The lack of a goodness-of-fit measure to apply to skeletons output by these algorithms represents a major problem in this field.

The filament skeleton is a one pixel wide representation of the filaments structure, tracing the main filament and its branches. When using the same input image (e.g. a dust column density map) different combinations of input parameters to filament finding algorithms produce different output skeletons. Different algorithms also return different skeletons. This is illustrated in Figure 1 and Figure 2 where we present the input image (Herschel dust column density map for the South-ridge of Vela C of Fissel et al. (2016)), alongside the output skeletons produced by filfinder and disperse with different input parameter combinations.

The filament skeleton forms the foundation for other analyses. These include measurements of length, width and mass per unit length (e.g. Arzoumanian et al., 2011), and relative orientation comparisons with respect to magnetic fields and clumps (Green et al. 2017, in preparation). These analyses are then interpreted in terms of supporting or rejecting star formation hypotheses. It is therefore essential that the filament skeleton is as good a representation of the structures in the input image as possible.

Current works involving filament skeletons can be divided into two groups: 1) those that repeatedly perform an analysis (like width fitting) on multiple skeletons for the same input image, and 2) those that perform the analysis on a single skeleton. In general the former perform a parameter study222By parameter study we refer to the process of repeatedly running the filament finding algorithm with different combinations of input parameters, which sample the sensible values for those inputs., and visually identify the subset of skeletons that most reasonably reproduce the structures within the original input image (e.g. Panopoulou et al., 2014). Analyses are then performed on each skeleton in that selected subset. The latter may also perform a parameter study, and then either a single output skeleton corresponding to one set of input parameters is chosen visually as most closely reproducing the structures of the original image (e.g. Hill et al., 2011) or constraints on a property (such as filament/skeleton length) may be used to whittle the larger pool of output skeletons down to a smaller set, which can then be combined (e.g. Liu et al., 2016).

These approaches all involve an assessment of similarity between the output skeletons and the structures in the original input image. The previous lack of a mathematical goodness-of-fit measure forced astronomers to make this assessment visually, and to hence potentially introduce visual bias. We propose that the mean structural similarity index (MSSIM) can be used as a goodness-of-fit measure to find the mathematically most similar output skeleton(s) to the input image, solving this problem. The MSSIM will be discussed in Section 2. This measure makes it possible to examine the relative similarities of the output skeletons (with respect to the input image) as a function of the input parameters, in order to find the range of optimum input values in parameter studies, as will be demonstrated in Section 3.

2 Mathematical similarity with the MSSIM

2.1 Current filament finding algorithms

filfinder is a new algorithm for filament identification (Koch & Rosolowsky, 2015). It has four independent input parameters that influence the output: the “skeleton threshold”, the “branch threshold”, the “global threshold”, and the “flattening threshold”. These parameters respectively set the minimum length of the skeletons in pixels, minimum length of the branches in pixels, a noise threshold as a percent, and a threshold for the arctan flattening step as a percent (which removes the effect of compact, bright regions like cores from the filament mask). The output skeleton differs significantly when different input parameters are used, as illustrated in Figure 1.

disperse is another algorithm used for filament identification (Sousbie, 2011). With this algorithm the user sets an input image and two main input parameters called the “persistence” and “robustness” thresholds. disperse then outputs a skeleton. The persistence threshold is used to filter noise. Components of the topology of the data are represented by pairs of critical points, one negative and one positive, and these are called persistence pairs. The absolute difference of the pair is its persistence. Noise only creates or destroys topological components of persistence lower than its local value. When setting the persistence threshold the user is removing topological components with persistence lower than that threshold and is thus filtering out noise. The robustness parameter is a measure of the contrast between the filaments and the background. The robustness parameter can be set as a ratio called the robustness ratio. Again, the output skeleton differs significantly when different persistence and robustness ratio thresholds are set as shown in Figure 2, as would be the case for other algorithms.

Commonly a parameter study is performed with these algorithms to explore the effect of different input parameter combinations. These parameter studies can produce thousands of output skeletons. These may not all be good representations of the structures in the input image (e.g. the skeleton in panel ix) of Figure 1). Additionally some of these skeletons may be mathematically better representations than others by some predefined metric. Visually, however, some may appear to be very similar (e.g. the skeletons in panels iii) and iv) of Figure 1), and many may seem to be reasonable representations.

2.2 Problems with visual selection of skeletons

Identifying individual skeletons that are visually reasonable is tedious and time consuming. It is also difficult to justify why skeleton A was chosen over skeleton B when both are similar. Not only does it introduce visual bias, but the subjective nature of visual selection of skeletons means that it is often difficult to reproduce results upon repetition. The application of the MSSIM as a mathematical goodness-of-fit measure, addresses all of these problems.

With the imminent influx of enormous amounts of data that astronomers face in the near future, automation of data processing is essential. The MSSIM further facilitates the automation of filament identification by eliminating the need for manual, visual selection (or rejection) in the post-processing of the skeletons. The application of the MSSIM is particularly advantageous for large datasets e.g. ALMA studies of filamentary regions.

2.3 Solution: the MSSIM

The similarity between two images (a ‘perfect’ reference image and a test image) with the same dimensions can be measured using the MSSIM (Wang et al., 2004). The MSSIM is a value between 1 (perfect match) and -1 (no match), providing a very simple and easy to interpret similarity measurement (Wang, private communication 2017). It is designed to be a universal index independent of the images being compared and the individual observers (Wang et al., 2002). It is symmetric so the same MSSIM would be obtained if the reference and test image were switched in the comparison.

The local structural similarity index (SSIM) measures the similarity between two images, $a$ and $b$, in terms of “luminance” (mean intensity, $l(a,b)$), “contrast” (standard deviation of intensity, $c(a,b)$) and “structure”333Structure is also defined as the spatial correlation and/or connectedness of pixels in the image that represent objects in the scene. In the context of the MSSIM it can also be described as any change between the two images that is not a simple shift in luminance or contrast. (pattern stored in the image, $s(a,b)$), on a pixel-by-pixel basis. For each pixel the multiplicative combination of these three terms results in the local SSIM (Wang et al., 2004):

[TABLE]

where $\mu_{a}$, $\mu_{b}$ are the local means, $\sigma_{a}$, $\sigma_{b}$ are the standard deviations and $\sigma_{ab}$ is the cross-variance for the two images. The constants $C_{1}$, $C_{2}$ and $C_{3}$, are small values to avoid the situation where the denominators would otherwise be equal to zero. In Equation 1, $\alpha=\beta=\gamma=1$ and $C_{3}=C_{2}/2$, simplifying the local SSIM to (Wang et al., 2004):

[TABLE]

The local SSIM values for each pixel over the images are then averaged to give the mean SSIM (MSSIM). The MSSIM is “mostly insensitive” (Brooks et al., 2008) to differences between images due to changes in luminance and contrast. The luminance term is present to take care of the fact that local differences between images are less obvious in brighter regions. The contrast term accounts for local differences between images being less obvious when there is ‘texture’444Texture in this context can be described as rapid changes in contrast between spatially proximate pixels. The contrast term accounts for the fact that it is hard to see distortion in a highly textured image because of rapid changes in contrast, e.g. in an image of gravel it would be hard to see a distortion because of the changes in contrast between lighter rocks and their darker shadows. in the image. According to Brooks et al. (2008) the MSSIM is instead sensitive to the changes between images that “break down natural spatial correlation” (i.e. to changes in structure), so can be used effectively to compare the structures in the original image input to the filament finding algorithm to those in the output skeleton images.

The MSSIM was originally used to measure image quality in television transmission, comparing the originally transmitted (reference) image to the received (degraded test) image. In this context, the output skeletons can be viewed as a very degraded version of the input image. The MSSIM can be used to compare the skeleton (degraded test) image to the input (reference) image to measure the similarity or ‘quality’ of the skeleton with respect to the input image. The MSSIM has been used in medicine for a similar purpose, to turn endoscopic videos into still images by finding video frames that are most different to each other, reducing a lengthy video into a number of still images that form the basis of a diagnosis (Low et al., 2011). In this letter we propose for the first time that the MSSIM be utilised in the post-processing of identified filament skeletons to find the skeletons most similar to the input image and thus apply it as a goodness-of-fit measure.

2.4 Using the MSSIM

The individual details of projects in filament identification differ. What they have in common is that many output skeletons are produced by different input parameters (where one combination produces one skeleton, and a different combination results in a different skeleton) and that these skeletons are usually saved as image arrays, generally as Flexible Image Transport System (FITS555http://fits.gsfc.nasa.gov/fits_primer.html) files. The MSSIM can then be calculated between the original greyscale image and the black and white output skeletons using inbuilt functions within matlab or the python scikit-image library (van der Walt et al., 2014).

The MSSIMs will not be large (approaching unity) as a binary black/white skeletal image is being compared to a greyscale image including more diffuse structures. In our examples in Figure 1 and  Figure 2 the resultant MSSIM values for the different skeletons are very close. The differences in the MSSIM values are real, representing genuine differences in the relationship between the skeleton image and input image. We are mainly interested in a method of comparison between the images that is quantitative and repeatable so even these small differences in MSSIM are sufficient. This does not negatively affect the ability of the MSSIM to discriminate between skeletons to select a quantitatively-defined ‘best’ skeleton(s). The MSSIM effectively standardises the selection of the ‘best’ skeleton(s) by rigorous mathematical and reproducible means.

The MSSIM is a method to compare the similarity of a skeleton to its corresponding input image. We suggest it is best utilised to find the optimum skeletons for a given input image and filament finding algorithm. It could also be used to find the best skeletal representations of the input image independently of the algorithm. However we urge caution in that approach because filament identification algorithms operate on different mathematical bases, implicitly defining filaments in different ways. They also have different input parameters and their skeletons may require different amounts of post-processing before they meet the users definition of a filament (e.g. in the removal of skeleton sections that are too short to be a filament by the users definition) and before the outputs from different algorithms are directly comparable. Additionally the MSSIM does not indicate that a skeleton and the input parameters which created it are physically sensible, nor does it indicate anything about the ‘goodness’ or suitability of an algorithm, or input image (e.g. in terms of noise properties), for the purpose of filament identification. It only quantifies you how good a representation of the input image the corresponding skeleton is.

We have proposed that the MSSIM method can be used for skeleton-based algorithms, but it could also be applied for the same purpose to the outputs of non-skeleton based filament identification algorithms. These algorithms, such as the Schisano et al. (2014) Hessian matrix approach, generally produce a ‘region of interest’ (ROI) mask, in FITS file format, that defines the extent (in terms of length and width) of the detected filaments. The ROI mask is later applied to the original image to extract only the filamentary regions. Like filament skeletons, the ROI mask is a representation of the filaments in the image. To measure how good a representation of the filaments an ROI mask is, it could be compared back to the original input image using the MSSIM. Alternatively the mask could be skeletonised, and then compared. In this way the ROI mask(s) best representing the input image could be found from the pool output by the algorithm.

3 Parameter optimisation with the MSSIM

The MSSIM makes it possible to perform systematic parameter studies exploring the impact of input parameters on skeletons. It can be used to locate the optimum range of input parameters for filament finding. An example of input parameter optimisation using the MSSIM is presented in Figure 3 for both disperse and filfinder. We use the same input image as previously.

For disperse we first perform a coarsely incremented study in Figure 3 panel i), exploring persistence in the range of 0.1-4.1$\times$1022 cm*-2* with an increment of 0.2$\times$1022 cm*-2*. The robustness ratio is varied from 0 (proxy for unset) to 10 with an increment of 0.5. That figure shows that the highest MSSIMs are located in the persistence range of 0.1-2.1$\times$1022 cm*-2*, and a robustness ratio range of 3-8. We then zoom in on this range in panel ii) with a more finely incremented study, with increments of 0.1$\times$1022 cm*-2*, and 0.25 respectively. This two-stage approach of using a coarsely incremented study over a large parameter range then narrowing that range and performing a fine parameter study saves computing time compared to performing a finely incremented study over the whole parameter range. From the resulting pool of skeletons we can locate the one with the highest MSSIM, which was produced by a persistence of 0.1$\times$1022 cm*-2* and a robustness ratio of 4.5.

A similar parameter study can be performed for filfinder. The four parameters could be optimised in a four dimensional study, and the MSSIM used to identify the optimum skeletons from the pool produced. Alternatively the flattening and global thresholds could be optimised with the MSSIM in a two dimensional study, holding the skeleton and branch thresholds constant at values appropriate for the dataset (these are minimum cutoffs so could be reasonably set to minimum values deemed suitable). We show an example of a two stage, two dimensional filfinder parameter study in  Figure 3. Panel iii) of that figure shows the coarsely incremented parameter study, exploring both the global and flattening thresholds from 0-100%, with an increment of 5%. For the purpose of this illustration the skeleton and branch thresholds are held constant at 10 pixels (0.3 pc, corresponding to the filament definition of an aspect ratio of 3 (Panopoulou et al., 2014), given an assumed width of 0.1 pc (Arzoumanian et al., 2011)) and 6 pixels (0.18 pc) respectively. That panel shows that the highest MSSIMs occur in the global threshold range of 60-80% so we zoom in on this range in the finely incremented parameter study shown in panel iv) with an increment of 1%. The skeleton with the highest MSSIM in that study was produced by a global threshold of 70% and a flattening threshold of 60%. Regardless of the filament identification algorithm used, the MSSIM can be utilised to optimise the input parameters and identify the optimum skeleton(s) for subsequent analyses.

4 Summary

Previously there has been no ‘goodness-of-fit’ measure to apply to skeletons output by filament finding algorithms. Different input parameter combinations can produce thousands of output skeletons. Different algorithms also produce different skeletons. Until now the only way to select skeletons for further analysis has been visually. This comes with inherent problems, including visual bias, large time-cost and limited reproducibility. We propose the application of an existing measure used in television transmission and in medical imaging, the mean structural similarity index (MSSIM), as a ‘goodness-of-fit’ measure for the output skeletons of filament identification algorithms. The MSSIM allows the best output skeletons for a given input image to be found for a given algorithm, or independently of the algorithm. This measure makes possible the optimisation of parameter space searches (parameter studies). It enables the user to examine the similarity of the output skeletons to the input image and immediately zoom in on the region of parameter space in which best fit solutions are located. The MSSIM for the first time presents a method by which best fit filament skeletons may be found in a completely automated and reproducible way for large volumes of data. Such a measure is essential in this era of ‘big data’.

Acknowledgements

We sincerely thank the referee for their helpful comments that improved this letter. Herschel is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA. 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. LMF 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 astropy (http://www.astropy.org) (Astropy Collaboration et al., 2013), aplpy (http://aplpy.github.com), Karma visualisation tools (Gooch, 1996), scikit-image (van der Walt et al., 2014) and NASA’s Astrophysics Data System.