Noise reduction and quantification of fiber orientations in greyscale images

Summary

Quantification of the angular orientation distribution of fibrous tissue structures in scientific images benefits from the Fourier image analysis to obtain quantitative information. Measurement uncertainties represent a major challenge and need to be considered by propagating them in order to determine an adaptive anisotropic Fourier filter. Our adaptive filter method (AF) is based on the maximum relative uncertainty δ_cut of the power spectrum as well as a weighted radial sum with weighting factor α. We use a Monte-Carlo simulation to obtain realistic greyscale images that include defined variations in fiber thickness, length, and angular dispersion as well as variations in noise. From this simulation the best agreement between predefined and derived angular orientation distribution is found for evaluation parameters δ_cut = 2.1% and α = 1.5. The resulting cumulative orientation distribution was modeled by a sigmoid function to obtain the mean angle and the fiber dispersion. A comparison to a state-of-the-art band-pass method revealed that the AF method is more suitable for the application on greyscale fiber images, since the error of the fiber dispersion significantly decreased from (33.9 ± 26.5)% to (13.2 ± 12.7)%. Both methods were found to accurately quantify the mean fiber orientation with an error of (1.9 ± 1.5)° and (2.3 ± 2.1)° in case of the AF and the band-pass method, respectively. We demonstrate that the AF method is able to accurately quantify the fiber orientation distribution in in vivo second-harmonic generation images of dermal collagen with a mean fiber orientation error of (6.0 ± 4.0)° and a dispersion error of (9.3 ± 12.1)%.

Keywords:

Imaging techniques – Collagens – Image processing – Grayscale – Aspect ratio – Signal filtering – Fourier analysis – Monte Carlo method

Introduction

The evaluation of the angular distribution of structures in scientific greyscale images is of major importance for various applications like in the analysis of soft tissue fibers e.g. in [1–15], textiles [16, 17], electrospun scaffolds [18–20] or even reinforced concrete [21]. Knowing the angular distribution in fiber reinforced materials gives meaningful insights into their mechanical functionality [22]. For example, in case of biological tissue, the orientation distribution of collagen fibers can be directly inserted into biomechanical material models for finite element simulations [23, 24].

Fiber images from different image modalities such as scanning electron microscopy [25], histology [1], and laser scanning microscopy as e.g. in [2, 3, 7, 8, 10, 15, 22, 26–29] provide diverse image properties in terms of sharpness, contrast and fiber appearance. Thus, a variety of automated image processing techniques were developed including ellipsoidal fitting in the spatial domain [30], fiber tracking [31], the structure tensor method [10] and Fourier domain image processing [7, 9, 12, 32–34].

Looking into the Fourier method in more detail, it can be split into four major steps: image preprocessing, Fourier transform and filtering, calculation of the angular distribution and its quantification.

In the Fourier domain it is of key importance to reduce noise as it was shown that rotationally symmetric band-pass filtering significantly improve the accuracy of the method [6, 33, 34]. In mechanical experiments, (e.g. stretching of fiber reinforced material), fiber properties such as angular distribution and diameter get modified [2, 35]. Accordingly, isotropic as well as anisotropic signal responses in the Fourier domain has to be accounted for, which requires an adaptive anisotropic filtering. Quantification of the angular orientation distribution in terms of mean angle and fiber dispersion is commonly realized by fitting a semi-circular von-Mises distribution to the angular orientation distribution [6, 10, 12, 22, 30, 34, 36, 37].

Here, we investigate the Fourier method by exploiting objective approaches at any of the four mentioned steps. To approach the requirement of adaptive anisotropic filtering we introduce an image filter, which is based on the propagation of measurement uncertainties through the discrete Fourier transform. And finally, in terms of an improved quantification of the angular orientation distribution, we test a sigmoid function to model the cumulative orientation distribution.

To capture the performance of our method, we quantitatively compare it to a band-pass method, introduced by Morrill et. al [34], which was proven to provide an accurate quantification of the angular orientation distribution. Based on realistic Monte-Carlo simulated greyscale fiber images we observe the evolution of accuracy with respect to fiber width, fiber aspect ratio, degree of alignment and image quality.

To validate the applicability of our method on real images, we quantify the mean fiber orientations and dispersions in vivo second-harmonic generation (SHG) images of collagen fibers of human skin [38].

Material and methods

Any calculations were performed using MATLAB [39]. The Image Processing Toolbox [40] as well as the Curve Fitting Toolbox [41] were applied.

Note that in the following the term uncertainty is used for statistical measurement errors, whereas the term error is associated with the deviation of a value to its reference.

Fourier transform and uncertainty propagation

Let I(x, y) be an image with x ∈ [0, X] and y ∈ [0, Y]. The discrete Fourier transform I ( x , y ) → F [ I ( x , y ) ] = I ^ ( u , v ) is given by:

The intensities of the angular distribution are calculated by evaluating the centered power spectrum of I(x, y), which is defined as the squared magnitude of I ^:

Any intensity image exhibits a specific intensity uncertainty ΔI(x, y), which is at least equal to the photon counting uncertainty Δ I ( x , y ) = I ( x , y ) assuming Poissonian statistics [42]. The uncertainty of the real and imaginary part, Δℜ and Δℑ are given by propagating Eqs 2 and 3:

Since Δ ℜ [ I ^ ] and Δ ℑ [ I ^ ] both depend on the image uncertainty ΔI, the covariance Δ I ^ ℜ ℑ has to be taken into account [43, 44]:

The calculation of Δ P ( u , v ) is straight forward:

Noise simulation

To verify the validity of our image transformations and uncertainty calculations, a Monte-Carlo based noise simulation with different test images was carried out.

The uncertainty of each pixel, Δ I ( x , y ) = I ( x , y ), is assumed as a normal-distributed fluctuation of a repeated measurement I_k(x, y) around the measured intensity I ( x , y ) = 1 N ∑ k = 0 N I k ( x , y ) = 1 N ∑ k = 0 N ( I ( x , y ) + δ I k ( x , y ) ) with 1 N ∑ k = 0 N δ I k ( x , y ) = 0. N is the number of measurements (Fig 1).

The propagation of the image uncertainty ΔI(x, y) through an arbitrary image operation I(x, y) → C(I(x, y)), ΔI(x, y) → ΔC(I(x, y)) is compared to the standard deviation of the Monte-Carlo simulated images ΔC_MC:

The result of Eq 9 is then compared to the calculated uncertainty ΔC by computing the relative deviation as:

Preprocessing and filtering

Artifact removal

The computation of the discrete Fourier transform (Eq 1) intrinsically assumes a periodical continuation of the image causing cross-like artifacts to appear in the Fourier domain mainly along the major axis (Fig 2). The magnitude of these artifacts depends on the image intensities near the boundary. Weighting functions in the spatial domain such as the Hamming or Welch window, which gradually reduce the image intensity towards the boundary, are able to remove these artifacts [6, 12, 16, 32, 36, 45]. However, applying window functions is a trade-off between reducing the image content and removing the artifacts in the Fourier domain. To overcome the drawback of weighting functions in the spatial domain, the linear decomposition of the image I(x, y) → per[I(x, y)] into a smooth component s(x, y) and a periodic component p(x, y) = per[I(x, y)] with I(x, y) = p(x, y) + s(x, y) as proposed by Moisan [45] is used here (Fig 2). Since the periodic component per(I) differs only slightly from I (Fig 2), the effect of the periodic mapping on ΔI is believed to be negligible. Verification of the effect of the periodic mapping on the image uncertainty is carried out by applying the Monte-Carlo noise simulation and evaluating Eq 9 to three different test images. Classic image processing greyscale test images including the Lena, Cameraman and Boat image were chosen [46]. The Boat image serves as example in several figures. The image uncertainty ΔI is assumed as Δ I = I, which is indeed an arbitrary choice but since the focus here is on validation only the assumption is reasonable.

The overall mean deviation (Eq 10, N = 10⁵ iterations) averaged over all images and the entire image range amounts for Δ_MC = (0.01 ± 2.21)%. The image uncertainty exceeds the Monte-Carlo calculated uncertainty at the boundary by (23 ± 6)%. Thus, the original image uncertainty is slightly underestimated due to the uncertainty of the remaining pixels by −(0.23 ± 0.61)%. This deviation is considered as sufficiently low to reasonably assume Δ(per(I)) = ΔI. Detailed results for every test image are listed in the supplementary (S1 Table).

Power spectrum filtering

The same set of test images is used for validating the calculated uncertainty of the real and imaginary parts (Eqs 5 and 6). The averaged relative deviation (Eq 10) amounts for (0.00 ± 0.22)% for the real and imaginary parts of all images.

Lastly, the Monte-Carlo method was applied on Eq 8 for the same images and same number of iterations. Averaging the relative deviation Δ_MC over the entire image yields a maximum deviation of −(59.32 ± 26.48)% (S2 Table). However, Δ_MC correlates well with the relative uncertainty of the power spectrum Δ P / P (Fig 3), which ranges from minor values up to relative uncertainties above 800%. Restricting the area of evaluation to pixels which exhibit a maximum relative uncertainty of 100% increases the accuracy of the uncertainty calculation to −(8.09 ± 5.46)% at maximum (S2 Table, Fig 4). Excluding values with a high relative uncertainty naturally filters the power spectrum, while leaving the back-transformed image relatively unaffected.

In order to find the optimal cut-off value δ_cut, images with known ground truth have to be used. Since fibrous images are in the scope of this work, Monte-Carlo generated greyscale fiber images serve as test images for determining the optimal cut-off value.

Monte-Carlo image generation

Besides the Monte-Carlo simulation for noise simulation, another Monte-Carlo routine was implemented for the generation of images containing fibers with a known angular distribution.

A single fiber is defined by its width, length, orientation and location. While several publications dealt with the generation of binary fiber images [12, 34] we focus on the generation of random greyscale intensity images including noise knowing that the contrast of the fibers influence their contribution within the power spectrum [36]. Since our goal here is to find the optimal evaluation parameters for a large variety of realistic images similar to the SHG images of dermal collagen that are used later on, the usage of binary images is not reasonable. In addition, our approach aims on evaluating images without any preprocessing, thus the test images should cover different image qualities.

The orientation of the fibers is sampled from a semi-circular π-periodic von-Mises distribution:

To enforce different image qualities, we introduce a noise factor. The noise factor ranges from a minimum of 0, which corresponds to a completely denoised image up to a maximum of 1 which corresponds to random speckle noise which can value up to half of the maximum intensity of the image.

Angular distribution generation

Quantifying the angular orientation of the fibers by means of mean angle and fiber dispersion requires the calculation of the total intensity of each angle of the filtered power spectrum. This is realized by using a radial summation. The total intensity of a certain angle is given by the sum over all pixels of the power spectrum, that are included within the angular range [−δθ, +δθ]. The contribution of each pixel is weighted by the percentage of area, which is included in [−δθ, +δθ].

The uncertainty propagates as:

The calculated intensity is normalized such that ∑ θ I ( θ ) = 1. An issue that is faced here are high intensity pixels close to the center of the power spectrum, which do not provide any directional information. The intensity of these pixels is several magnitudes higher than the intensities of interest, which causes artifacts in the angular distribution. A sensitivity analysis of a set of test images showed that zeroing pixel with a distance smaller than 3 pixels from the center is sufficient to remove these artifacts (S1 Fig).

Additionally, a modified intensity is defined, which exploits the anisotropy of the introduced filter. Let N_Δθ be the number of non-zero pixels within the angular range Δθ of the filtered power spectrum P ′. The modified intensity then reads as:

To find the optimal parameter set [θ_cut, α], a two parameters optimization algorithm is applied on N = 10⁴ Monte-Carlo generated images. For this purpose we use the MATLAB-implemented Nelder-Mead simplex method (fminsearch) [47]. To account for different types of images, fiber properties, namely width, aspect ratio and dispersion as well as the noise factor are randomly sampled from the respective interval specified before. We use the mean squared difference (MSD) between the calculated orientation distribution I MC ( θ ) and the prescribed orientation distribution I ( θ ) as objective function, which is minimized:

Termination was enforced as soon as the difference between consecutive iterations of both parameters was smaller than 10⁻³. MATLAB uses the same termination criterion for both parameters. Therefore, we scaled α by a factor of 0.1 to match the magnitude of δ_cut. A subset of 10³ images was evaluated prior optimization on a 10 × 10 grid with α ∈ [−1, 8] and with δ_cut ∈ [1, 100]% to ensure that the calculated minimum is global and to get an estimate for the starting values of each parameter. We estimated that α = 2 and δ_cut = 2% provide a suitable set of starting values.

Angular distribution quantification

Conventional approach

The common approach for quantifying the mean orientation and dispersion is to fit a semi-circular von-Mises function (Eq 11) to the distribution I ( θ ). This approach provides accurate results in case of aligned fibers but looses accuracy in case of isotropic distributions [34]. Additionally, the use of an arbitrary averaging range to smooth the angular orientation distribution prevents a meaningful and objective interpretation of the data in terms of e.g. the number of dominant fiber directions.

New fitting approach

To find an objective evaluation procedure which reliably can deal with unprocessed data, fitting the cumulative distribution function (CDF) C ( θ ) is quite appealing, since the data are naturally smoothed:

The uncertainty of Δ C is given by a quadratic summation:

Since the choice of the starting angle of the summation (Eq 17) is arbitrary, the uncertainty Δ C must be independent from θ. Thus, the maximum uncertainty of C ( θ ) is set as uncertainty for all angles. A peak in the angular distribution corresponds to a step in the cumulative distribution, which approaches a first order polynomial for isotropic distributions. To model the cumulative orientation distribution a sigmoid function is chosen:

The advantage of fitting a von-Mises distribution is its semi-circularity, namely P_VM(θ) = P_VM(θ + 180°). Since the cumulative distribution function is monotonically increasing, the corresponding condition reads as C ( θ ) = C ( θ - 180 ° ) + 1 = C ( θ + 180 ° ) - 1. In order to meet this condition, Eq 19 is modified by adding neighboured sigmoid functions accounting for the added offsets:

Validation

Comparison to band-pass filtering

In the following we refer to our method as AF method (adaptive filtering method) for convenience. To capture the performance of the implemented procedure we compare the AF method to the state of the band-pass method, which has been proven to provide more accurate results than other methods [34]. In order to observe the influence of fiber and image properties on the accuracy of both methods, we created a dataset, where each property is altered separately using the introduced Monte-Carlo method for generating greyscale images. The following sets of values were used for creating the image dataset:

In order to enable a reasonable statistical evaluation, 20 images with random mean fiber orientation θ ¯ are created for each category, which adds up to a total of N = 3³ ⋅ 21 ⋅ 20 = 11, 340 images. Each image is evaluated using the AF method using the optimal evaluation parameters δ_cut = 2.1% and α = 1.5 and using the band-pass method of Morrill et. al by applying their provided software FiberFit. We follow their recommendations of the upper and lower cutoff parameter t = 2 and t = 32 yielding a lower cutoff frequency of 8 and an upper cutoff frequency of 128. The resulting evaluation parameters, mean orientation θ ¯ and dispersion coefficients b and k are compared to the reference parameters. To account for Monte-Carlo sampling errors, reference mean orientation and distribution parameter k are obtained from fitting the frequency distribution of sampled fiber angles. The reference dispersion parameter b is obtained from fitting the cumulative frequency distribution of sampled fiber angles using the modified sigmoid function (Eq 20). Mean orientations θ ¯ F F of the FiberFit software were inverted and shifted to match our coordinate system θ ¯ F F ′ = 180 ° - θ ¯ F F. Additionally, angles exceeding the interval [0°, 180°] were shifted by 180°. Comparisons to the reference parameters are performed by considering the absolute error of the mean orientation angle angle Δ θ ¯ and the absolute relative error of the dispersion parameters Δb_rel and Δk_rel. A small fraction of images were found to exhibit large relative errors Δb_rel and Δk_rel. Those are classified as outliers if they are exceeding three times the interquartile range of the respective distribution.

Statistics

A paired t-test is used to calculate the level of significance (p-value) between the error of both evaluation methods for the same subset of images in terms of width, aspect ratio, dispersion and noise. If a value is classified as outlier, the corresponding value of the pair is excluded for p-value calculation. An unpaired t-test is used for p-value calculations among groups with different noise factor. Significance levels of 0.05, 0.01 and 0.001 are marked as (*, **, ***), respectively.

Application on experimental data

Finally, we checked the applicability of the AF method on in vivo experimental data. We use SHG-images of dermal collagen that were recorded using a CE-certified multi-photon microscope (DermaInspect) for in vivo applications, which was developed in collaboration with Jenlab GmbH (Jena, Germany). The SHG signal was captured using an excitation wavelength of 820 nm together with a specific band-pass filter (410 ± 10 nm, AQ 410/20m-2P, Chroma Technology Corp., Bellows Falls, VT). A scan time of 13 s and image dimensions of 512 × 512 pixels with a 220 × 220 μm field of view were chosen as image acquisition parameters. For a detailed description of the microscope refer to [48–50]. In total, ten SHG images of dermal collagen were taken from ten different volunteers at a depth of 30 μm under the basal membrane located at the cheek.

This study was conducted according to the recommendations of the current version of the Declaration of Helsinki and the Guideline of the International Conference on Harmonization Good Clinical Practice, (ICH GCP). In addition, this study was approved and cleared by the institutional ethics review board (Beiersdorf AG, Hamburg, Germany). Written informed consent was obtained from each volunteer.

To get the reference angular distribution we manually trace the collagen fibers. Statistical variance is achieved by tracing each image three times in random order. The angular orientation distribution is achieved from adding up the length of each fiber being oriented along a certain angle in 1° steps. A smooth distribution is obtained by filtering the data using a Gaussian kernel with a standard deviation of 1°. Reference parameters θ ¯ m and b_m are obtained by fitting the modified sigmoid function (Eq 20).

Results

Angular distribution generation

Fig 5 shows the result of the optimization procedure of the evaluation parameters δ_cut and α. Apart from a minor fraction of outliers, data points of minimal difference spread within δ_cut < 4% and 0 < α < 3. Since the frequency distributions p(δ_cut) and p(α) are not normally distributed, the median value of both parameters is considered as estimation of the overall minimum (δ_cut = 2.1%, α = 1.5). The mean squared difference minimizes for α in [1, 3] and for δ_cut < 10% indicating that a global minimum is considered. Subsequent calculations using the AF method are performed with δ_cut = 2.1% and α = 1.5. Median values of the mean squared difference as a function of fiber dispersion k, image quality (NF) and fiber geometry (width, aspect ratio) are provided in the supplement (S2 Fig).

Validation

Comparison to band-pass filtering

The overall error of the mean fiber orientation, Δ θ ¯ amounts for (2.2 ± 1.8)° and (1.8 ± 1.4)° (p < 0.001) for the band-pass and the AF method, respectively (Table 1). Note that for calculating the overall error Δ θ ¯ highly dispersed fiber networks (k < 1) were excluded from statistical analysis since they do not feature a mean orientation. The overall mean error of the dispersion parameters Δk_rel, Δb_rel amounts for (33.9 ± 26.5)% and (13.2 ± 12.7)% (p < 0.001) using the band-pass and the AF method, respectively. It was found that for both methods the error of the mean fiber orientation Δ θ ¯ mainly depends on the degree of the alignment of the fiber network (Fig 6A) and is rather independent from the fiber geometry or image quality (not shown). The error of both methods reduces towards an increasing degree of fiber alignment.

The error of the dispersion parameters Δb_rel and Δk_rel is found to strongly depend on the degree of alignment, fiber geometry and image noise (Figs 6B and 7). With increasing the degree of alignment, the fiber dispersion error of the AF method, Δb_rel, continuously decreases from a maximum error of 10.9% at k = 0.01 to a plateau of ∼5% for k > 2.2. The fiber dispersion error of the band-pass method, Δk_rel, exhibits a large error at k = 0.01 Δk_rel = 101.4%, which is first reduced to a minimum error of Δk_rel ∼ 12% around k ∼ 1 but then is increased again towards aligned fiber networks to an error of 27.7% at k = 5 (Fig 6B). Except for k = 1 and k = 1.25 the error of both methods shows a significant (p < 0.05), for most values of k even highly significant (p < 0.001) difference.

Other than the error of the mean fiber orientation, the error of the dispersion parameter strongly depends on the fiber geometry, image noise and choice of evaluation method (Fig 7). An at least significant decrease (p < 0.05) in error was found for every group if the image was evaluated with the AF method.

A decreased fiber width strongly increases Δk_rel for every combination of aspect ratio, dispersion and noise. In case Δb_rel, the increase in error can be noted for aligned networks only. A raised noise factor significantly (p < 0.001) increases Δk_rel and Δb_rel for aligned networks of fibers with a width of 1. Additionally, a significant (p < 0.05) decrease and a highly significant increase (p < 0.001) in Δb_rel was found for thick, short fibers (width = 10, AR = 15).

Application on experimental data

Our implemented algorithm was applied on an exemplary set of ten in vivo recorded SHG images of dermal collagen. Parameters of the cumulative angular distribution θ ¯ and b were calculated using the AF method and compared to reference parameters θ ¯ and b achieved from manual fiber tracing (Fig 8, S1 Dataset). The absolute mean error between calculation and manual segmentation for the mean orientation amounts for Δ θ ¯ = ( 6 . 0 ± 4 . 0 ) °, whereas the mean relative error of the fiber dispersion is Δb_rel = (9.3 ± 12.1)%. The mean coefficient of determination was R² = 0.99 ± 0.01.

Value pairs of calculated and reference parameter were fitted using a first order polynomial (Fig 9). High Pearson correlations were found for both parameters, namely R² = 0.99 and R² = 0.90 for the mean orientation angle and dispersion, respectively. Calculated Pearson correlations with respect to an ideal calculation (slope equal to one) amount for R² = 0.98 and R² = 0.8.

Discussion

We report on a robust method to quantify the angular distribution of fibers in noisy greyscale fiber image. The whole image processing procedure covers: the application of the periodic decomposition to remove cross-like artifacts in the Fourier domain, Fourier filtering by only permitting values below a certain relative uncertainty δ_cut, computation of the angular distribution by weighting the number of pixels with N^α and quantification of the mean angle and dispersion by fitting a modified sigmoid function to the cumulative orientation distribution.

In comparison to conventionally used window functions like in [12, 16, 32, 36], the periodic decomposition has the advantage of conserving almost the entire image information while completely removing artifacts in the Fourier domain, as shown in Fig 2. Therefore, we omit a quantitative analysis. The Monte-Carlo noise simulation revealed that the effect on the uncertainty can be neglected since the periodic mapping only has significant effects on the image boundary only.

Filtering the power spectrum by excluding values above a predefined relative uncertainty allows for the definition of an adaptive filter. Optimal evaluation parameters were identified by applying a two parameter Nelder-Mead optimization, which succeeded for 99.8% of the images. A maximum cut-off error of δ_cut = 2.1% and a weighting factor of α = 1.5 was calculated, while the non-locality of the optimum was ensured. Polzer et. al introduced a similar weighting factor which raises the entire intensity distribution I(θ) [36]. However, the optimum value of their weighting factor seems to suffer from large fluctuations, whereas the optimal value of α is stable throughout the degree of alignment of fiber networks and throughout different fiber properties and image noise (S2 Fig).

Previously used filtering methods like the bandpass method [6, 34, 36] work with rotationally symmetric filters, which disregard the anisotropy of the signal. The derivation of optimal band-pass range is based on the fiber diameter [6]. With the AF method we apply a relative uncertainty criterion without any assumptions on the underlying fiber properties. The optimal value of δ_cut is found to be completely independent from the chosen fiber geometries, fiber disperion and image noise (S2A, S2C, S2E and S2G Fig). As a consequence, the fiber dispersion can be quantified with a significantly lower error for completely dispersed and highly aligned fiber networks in comparison to the band-pass method (Fig 6B). Especially the quantification of highly dispersed fiber networks (k < 1) is much more reliable using the AF method, whereas the dispersion calculated by the band-pass method suffers from large errors (Δk_rel > 100%) (Figs 6B and 7). This is in accordance with the observations of Morrill et. al [34], who measured an error of ∼ 30% for k ∼ 0.2 using binary Monte-Carlo images.

If we use the AF method in conjunction with a von-Mises fit of the angular orientation distribution, large errors of the dispersion coefficient can be noted towards dispersed fiber networks (S3 Fig). Contrary to the band-pass method, Δk_rel rapidly decreases towards aligned networks to a level of ∼ 10%. This effect is solely related to the adaptive filter and the weighted radial summation of the AF method, whereas the sigmoid fit ensures a reliable dispersion parameter estimation for dispersed networks.

Fibers with a width of 1 somehow represent an exception in comparison to the images containing fibers with width > 1 pixel. Regarding the result of the optimization of α with respect to different image and fiber properties (S2 Fig), it can be seen that α is quite stable throughout the dispersion k, noise factor and different fiber geometries. α fluctuates within the interval [1, 2] except for fibers with a width of 1 pixel, where the median value α amounts for α_width=1 = 2.4. The application of the gauss-filter to dissolve sharp edges might even reduce the true fiber dimension down to a size which is near the resolution limit of power spectrum based methods resulting in large deviations of the dispersion for highly aligned fiber networks.

Using the Monte-Carlo approach to create artificial fibrous images for validation purposes is a common tool. However, it is rather difficult to draw a comparison to the accuracy of other methods found in the literature since mostly a low quantity of binary images was investigated [12, 34]. For example, Morrill et. al measured an overall error of (0.71 ± 0.43)° for the mean orientation and (7.4 ± 3.0)% for the fiber dispersion using binary fiber images, whereas we measured errors of (2.3 ± 2.1)° and (33.9 ± 26.5)% [34]. The use of greyscale images will generally result in a lower accuracy compared to the evaluation of binary images, since the superimposition of fibers might generate intensity deviations in the power spectrum.

The AF method was applied to in vivo SHG images of dermal collagen. The multi-photon images that were used provide a sufficiently high image intensity I(x, y) (16-bit), which comes along with a sufficiently low relative intensity error Δ I / I = 1 / I. Therefore, δ_cut = 2.1% provides a reasonable filtering value, which results in an accurate quantification of the angular orientation distribution in terms of mean fiber orientation and fiber dispersion.

Conclusion

The proposed adaptive filter method modifies common Fourier methods at different stages, namely artifact removal in the Fourier domain, filtering of the power spectrum and quantification of the angular signal.

The adaptive filter conserves the anisotropy of the angular signal in the Fourier domain, which ensures a stable error for disordered as well as highly aligned fiber networks. Using the cumulative distribution naturally averages the data, which spares out any averaging of the angular distribution. The high mean goodness of the fit R² = 0.99 which was measured for both, Monte-Carlo images and SHG-images, indicates that the modified sigmoid function provides a perfect model of the cumulative distribution function.

The adaptive filtering method was found to be a reliable and accurate tool for quantifying the angular orientation distribution in fibrous SHG images of dermal collagen, even for images suffering from a low image quality. Aside from its benefits concerning accuracy and reliability, the AF method considers measurement uncertainties, which are of key importance in any kind of scientific experiment.

Supporting information

S1 Table [pdf]
Periodic decomposition: Relative deviation of the error propagation simulated by Monte-Carlo and the basic image error Δ for = 10 iterations.

S2 Table [pdf]
Power spectrum: Relative deviation between calculated and Monte-Carlo simulated uncertainty for = 10 iterations.

S1 Fig [a]
Sensitivity analysis of the central cut-off radius.

S2 Fig [ar]
Optimal evaluation parameters vs. fiber properties and image noise.

S3 Fig [eps]
Error of the dispersion parameter and of the von-Mises fit and the sigmoid fit of the AF method for fibrous Monte-Carlo images (width = 5, AR = 30).

S1 Dataset [csv]
The .csv file provides the reference parameter as well as the parameter calculated by the AF method for each image.