Noise reduction and quantification of fiber orientations in greyscale images

Abstract

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

About the Authors:

Maximilian Witte

Roles Conceptualization, Software, Validation, Visualization, Writing – original draft, Writing – review & editing

Affiliations Center for Free-Electron Laser Science (CFEL), University of Hamburg, Hamburg, Germany, Beiersdorf AG, Hamburg, Germany

ORCID logo http://orcid.org/0000-0002-2923-8907

Sören Jaspers

Roles Supervision

Affiliation: Beiersdorf AG, Hamburg, Germany

Horst Wenck

Roles Supervision

Affiliation: Beiersdorf AG, Hamburg, Germany

Michael Rübhausen

Roles Conceptualization, Supervision, Writing – review & editing

Affiliation: Center for Free-Electron Laser Science (CFEL), University of Hamburg, Hamburg, Germany

Frank Fischer

Roles Conceptualization, Data curation, Supervision, Writing – review & editing

* E-mail: [email protected]

Affiliation: Beiersdorf AG, Hamburg, Germany

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 is given by:(1)where the real and imaginary parts read as(2)(3)

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 :(4)

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

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

(7)

The calculation of is straight forward:(8)

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, , is assumed as a normal-distributed fluctuation of a repeated measurement Ik(x, y) around the measured intensity with . N is the number of measurements (Fig 1).

[Figure omitted. See PDF.]

Fig 1. Schematic outline of the pixel-wise noise simulation.

(A) The photon counting noise results in the standard deviation of each image pixel. In the pixel-wise uncertainty simulation, normal distributions are used to sample the intensity deviation δI(x, y) in each pixel, as exemplary shown in (B) and (C) for pixels (1, 2) and (1, 1), respectively. Thereby, high intensity pixels are associated with a large absolute deviation of potential intensity values (B).

https://doi.org/10.1371/journal.pone.0227534.g001

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 ΔCMC:(9)

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

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 , which is indeed an arbitrary choice but since the focus here is on validation only the assumption is reasonable.

[Figure omitted. See PDF.]

Fig 2. Effect of spatial domain image filters on the power spectrum.

Images in the spatial domain are shown in greyscale (A-C), whereas the respective power spectra are shown in false colors as well as in logarithmic scale for a better visibility (D-H). Power spectrum values are shifted, such that low frequencies are located in the center. (A) shows the original image boat [46]. (B) shows the periodic component of the periodic plus smooth decomposition of Moisan et. al [45]. Apart from small areas near the image boundary, the entire image information of the original image is conserved. However, as shown in (C), the Hamming window reduces the image information gradually towards the image boundary. (D), the discrete Fourier transform of images causes cross-like artifacts (black arrows) to appear in the power spectrum as the image is assumed to be periodical. (E-F), these artifacts disappear in the power spectra if firstly either the periodic decomposition or the Hamming window are applied. Note that the artifacts reappear in (G) and (H) where the absolute difference of the power spectra of the filtered images with respect to the power spectrum of the original image are shown. (G) also reveals that mostly non-directional information is removed from the power spectrum of the original image as only cross-like shapes are pronounced. Other than in (F), where the reduction of image information by the Hamming window also affects the power spectrum, which appears less sharp (white arrows) compared to (D) and (E).

https://doi.org/10.1371/journal.pone.0227534.g002

The overall mean deviation (Eq 10, N = 105 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 (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.

[Figure omitted. See PDF.]

Fig 3. Relative uncertainty between calculated and simulated uncertainty of the power spectrum.

The relative deviation between the theoretical uncertainty and the Monte-Carlo simulated uncertainty , of the Boat image (Fig 4). Note that the spectrum values are shifted such that low frequencies are located in the origin and that the original image was used here.

https://doi.org/10.1371/journal.pone.0227534.g003

[Figure omitted. See PDF.]

Fig 4. Relative uncertainty and filtering of the power spectrum.

The relative uncertainty (A) of the power spectrum increases with distance from the origin to values above 800%. A filter mask for the power spectrum is achieved from excluding relative uncertainties above δcut. (B-E) show the filtered power spectra for different δcut with the respective inverse Fourier transformations in (F-I). The number of remaining pixels after filtering as well as the error of the uncertainty calculation compared to the Monte-Carlo simulation ΔMC is shown in (J). The mean error, its standard deviation (grey shaded areas) as well as the number of remaining pixels gradually decrease towards a lower threshold δcut. Note that even if 99.5% of the pixels are excluded, the contours of the major image components are still visible (E).

https://doi.org/10.1371/journal.pone.0227534.g004

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:(11)with the dispersion parameters k and defining the width and the center of the distribution. I0(k) is the modified zero order bessel function . A rejection sampling algorithm is used to sample fiber angles. The intensity of each sampled fiber with certain width and aspect ratio is added to the existing image to obtain a greyscale image. The image is then smoothed using a gaussian kernel with standard deviation of two to slightly dissolve sharp edges. After that, intensities are scaled such that the maximum intensity is equal to the maximum intensity of a 16-bit intensity image. Width, aspect ratio and dispersion k of the fibers affect the accuracy of the image processing algorithm [34]. The minimum fiber width amounts for one pixel as the maximum fiber width is confined by the image size and the maximum allowed aspect ratio. As images at 512 × 512 pixels are generated, a maximum aspect ratio of 45 and a maximum fiber width of 10 are chosen to still allow for an effective placement of the fibers within the image boundaries. A minimum aspect ratio of 10 was proposed by Marquez et. al [6] to allow for a reasonable evaluation of orientation. As we are generating greyscale fiber images with overlapping fibers we choose a minimum aspect ratio of 15. Dispersed as well as aligned fiber networks are achieved by sampling k within [0.01, 5].

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 [−δθ, +δθ].(12)

The uncertainty propagates as:(13)

The calculated intensity is normalized such that . 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 . The modified intensity then reads as:(14)with the uncertainty given by:(15)where α defines the impact of NΔθ. A value of α = −1 corresponds to an averaged intensity, whereas α > 0 amplifies the number of remaining pixels within the given angular range. The ideal choice of α strongly depends on the chosen cut-off value δcut for the power spectrum. For example in case of a low cut-off value, the filtered power spectrum might exhibit a high anisotropy, where α > 1 might lead to a more accurate result compared to the unweighted sum. Considering a very high cut-off value, holds and the number of summed pixels should not have any influence. Thus α = −1 might be the value of choice. A somewhat similar approach was used by Polzer et. al [36], which used a factor to empower the entire intensity distribution .

To find the optimal parameter set [θcut, α], a two parameters optimization algorithm is applied on N = 104 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 and the prescribed orientation distribution as objective function, which is minimized:(16)

Termination was enforced as soon as the difference between consecutive iterations of both parameters was smaller than 10−3. 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 103 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 . 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) is quite appealing, since the data are naturally smoothed:(17)

The uncertainty of is given by a quadratic summation:(18)

Since the choice of the starting angle of the summation (Eq 17) is arbitrary, the uncertainty must be independent from θ. Thus, the maximum uncertainty of 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:(19)

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

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:(21)(22)

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 = 33 ⋅ 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 of the FiberFit software were inverted and shifted to match our coordinate system . 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 Δbrel and Δkrel. A small fraction of images were found to exhibit large relative errors Δbrel and Δkrel. 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 and bm 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).

[Figure omitted. See PDF.]

Fig 5. Effect of evaluation parameters δcut and α on the error.

False colours indicate the mean squared difference MSD between the computed angular distributions and their reference on a subset of N = 103 images for a discrete set of 10 × 10 parameter values. Data points represent the optimal parameter set for each image (N = 104) based on the exploited two parameters Nelder-Mead optimization algorithm which used the squared difference between reference and calculated orientation distribution as objective function. Histograms show the frequency distribution of parameter δcut and α. The median values of both distributions(δcut = 2.1%, α = 1.5) were identified as optimal evaluation parameters. The mean number of iterations, until the termination criterion was met, is 23 ± 9. Only 0.2% of the images reached the maximum number of iterations, which was set to 100.

https://doi.org/10.1371/journal.pone.0227534.g005

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 Δkrel, Δbrel 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.

[Figure omitted. See PDF.]

Fig 6. Error of the band-pass method and the AF method vs. the dispersion k of the fiber network.

Fibers with a width of 5, aspect ratio of 30 and an image noise factor of 1 were considered here. (A) shows the error of the mean fiber orientation . Inlets show exemplary Monte-Carlo generated greyscale fiber images from respective distributions. (B) shows the relative error of the fiber dispersion parameters Δkrel, Δbrel. Sample greyscale images, that were simulated by the implemented Monte-Carlo method are shown for k = 0.01 and k = 5. (*p < 0.05, **p < 0.01, ***p < 0.001).

https://doi.org/10.1371/journal.pone.0227534.g006

[Figure omitted. See PDF.]

Table 1. Overall error of mean orientation and dispersion determined by the band-pass method and by the AF method.

Note that f indicates the percentage of data points that were classified as outlier.

https://doi.org/10.1371/journal.pone.0227534.t001

The error of the dispersion parameters Δbrel and Δkrel 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, Δbrel, 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, Δkrel, exhibits a large error at k = 0.01 Δkrel = 101.4%, which is first reduced to a minimum error of Δkrel ∼ 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.

[Figure omitted. See PDF.]

Fig 7. Error of the dispersion parameter of the band-pass method and the AF method for different fiber geometries, dispersions and noise factors.

(A) shows the error of the dispersion parameter for fibers with a width of 1 pixel for aspect ratios 15 and 45 as well as noise factors 0 and 1. (B) shows the same as in (A) but for fibers with a width of 10 pixels. Each error is given for a dispersed (k = 0.01) and an aligned fiber network (k = 5). Inlets show exemplary Monte-Carlo generated greyscale fiber images from respective distributions. (*p < 0.05, **p < 0.01, ***p < 0.001).

https://doi.org/10.1371/journal.pone.0227534.g007

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 Δkrel for every combination of aspect ratio, dispersion and noise. In case Δbrel, the increase in error can be noted for aligned networks only. A raised noise factor significantly (p < 0.001) increases Δkrel and Δbrel 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 Δbrel 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 , whereas the mean relative error of the fiber dispersion is Δbrel = (9.3 ± 12.1)%. The mean coefficient of determination was R2 = 0.99 ± 0.01.

[Figure omitted. See PDF.]

Fig 8. Evaluation example of an in vivo SHG-image of dermal collagen.

Left column ((B),(D),(F)): Evaluation steps of the implemented image processing algorithm. Right column ((C),(E),(G)): Manual fiber tracing of fibers, which serves as ground truth. (A) shows the original SHG-image as measured with the multi-photon microscope. (B), the filtered power spectrum for δcut = 2.1%. (C), manually traced fibers. (D) and (E) show the resulting angular orientation distributions. For (D), α = 1.5 was used. (F) and (G) show the respective cumulative orientation distributions and the fit parameters obtained from sigmoid fitting.

https://doi.org/10.1371/journal.pone.0227534.g008

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 R2 = 0.99 and R2 = 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 R2 = 0.98 and R2 = 0.8.

[Figure omitted. See PDF.]

Fig 9. Calculated angular distribution parameters vs. reference distribution parameters for the experimental data.

(A) Mean fiber orientation vs. mean fiber orientation achieved from manual segmentation. (B) Fiber dispersion b vs. fiber dispersion bm achieved from manual segmentation. R2 values are given for the fitted curve (solid, orange) as well as for the ideal curve (dashed, black) with a slope of one. Errorbars represent 95% confidence bounds of fitted parameters.

https://doi.org/10.1371/journal.pone.0227534.g009

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 (Δkrel > 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, Δkrel 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 . 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 R2 = 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

[Figure omitted. See PDF.]

S1 Table. Periodic decomposition: Relative deviation of the error propagation simulated by Monte-Carlo and the basic image error ΔI for N = 105 iterations.

https://doi.org/10.1371/journal.pone.0227534.s001

(PDF)

S2 Table. Power spectrum: Relative deviation between calculated and Monte-Carlo simulated uncertainty for N = 105 iterations.

https://doi.org/10.1371/journal.pone.0227534.s002

(PDF)

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

(A) Squared difference between calculated angular orientation distribution and sampled fiber angle distribution on N = 104 images (Monte-Carlo generated). The same dataset was also used for evaluation parameter optimization. Unfiltered power spectra are used to isolate the effect of the cut-off radius. (B-D) Evaluation example for cut-off radii (0,3,6). The artifact dramatically reduces from a 0 pixel cut-off to a radial cut-off of 3 pixels. Cut-off radii above 3 pixels barely influence the angular distribution. In order to save angular information, we choose a cut-off radius of 3 pixels.

https://doi.org/10.1371/journal.pone.0227534.s003

(EPS)

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

Plotted are median values and errorbars represent the interquartile range (25, 75). (A,C,E,G) represent plots of δcut vs. dispersion k, noise factor NF, aspect ratio (AR) and width. (B,D,F,H) are plots of the weighting factor α vs. dispersion k, noise factor NF, aspect ratio (AR) and width.

https://doi.org/10.1371/journal.pone.0227534.s004

(EPS)

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

The von-Mises fit was applied to the angular orientation distribution as calculated by the AF method, whereas the sigmoid fit was applied to the respective cumulative orientation distribution.

https://doi.org/10.1371/journal.pone.0227534.s005

(EPS)

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

https://doi.org/10.1371/journal.pone.0227534.s006

(CSV)

Citation: Witte M, Jaspers S, Wenck H, Rübhausen M, Fischer F (2020) Noise reduction and quantification of fiber orientations in greyscale images. PLoS ONE 15(1): e0227534. https://doi.org/10.1371/journal.pone.0227534

References

1. Ní Annaidh A, Bruyère K, Destrade M, Gilchrist MD, Otténio M. Characterization of the anisotropic mechanical properties of excised human skin. J Mech Behav Biomed Mater. 2012;5(1):139–148. pmid:22100088

2. Bancelin S, Lynch B, Bonod-Bidaud C, Ducourthial G, Psilodimitrakopoulos S, Dokládal P, et al. Ex vivo multiscale quantitation of skin biomechanics in wild-type and genetically-modified mice using multiphoton microscopy. Sci Rep. 2015;5(December):1–14.

3. Chen X, Nadiarynkh O, Plotnikov S, Campagnola PJ. Second harmonic generation microscopy for quantitative analysis of collagen fibrillar structure. Nat Protoc. 2012;7(4):654–669. pmid:22402635

4. Frisch KE, Duenwald-Kuehl SE, Kobayashi H, Chamberlain CS, Lakes RS, Vanderby R. Quantification of collagen organization using fractal dimensions and Fourier transforms. Acta Histochem. 2012;114(2):140–144. pmid:21529898

5. Levillain A, Orhant M, Turquier F, Hoc T. Contribution of collagen and elastin fibers to the mechanical behavior of an abdominal connective tissue. J Mech Behav Biomed Mater. 2016;61:308–317. pmid:27100469

6. Marquez JP. Fourier analysis and automated measurement of cell and fiber angular orientation distributions. Int J Solids Struct. 2006;43(21):6413–6423.

7. Mega Y, Robitaille M, Zareian R, McLean J, Ruberti J, DiMarzio C. Quantification of lamellar orientation in corneal collagen using second harmonic generation images. Opt Lett. 2012;37(16):3312–4. pmid:23381241

8. Mercatelli R, Ratto F, Rossi F, Tatini F, Menabuoni L, Malandrini A, et al. Three-dimensional mapping of the orientation of collagen corneal lamellae in healthy and keratoconic human corneas using SHG microscopy. J Biophotonics. 2017;10(1):75–83. pmid:27472438

9. Nesbitt S, Scott W, Macione J, Kotha S. Collagen Fibrils in Skin Orient in the Direction of Applied Uniaxial Load in Proportion to Stress while Exhibiting Differential Strains around Hair Follicles. Materials (Basel). 2015;8(4):1841–1857.

10. Rezakhaniha R, Agianniotis A, Schrauwen JTC, Griffa A, Sage D, Bouten CVC, et al. Experimental investigation of collagen waviness and orientation in the arterial adventitia using confocal laser scanning microscopy. Biomech Model Mechanobiol. 2012;11(3-4):461–473. pmid:21744269

11. Ribeiro JF, dos Anjos EHM, Mello MLS, de Campos Vidal B. Skin Collagen Fiber Molecular Order: A Pattern of Distributional Fiber Orientation as Assessed by Optical Anisotropy and Image Analysis. PLoS One. 2013;8(1):5–7.

12. Schriefl AJ, Reinisch AJ, Sankaran S, Pierce DM, Holzapfel GA. Quantitative assessment of collagen fibre orientations from two-dimensional images of soft biological tissues. J R Soc Interface. 2012;9(76):3081–3093. pmid:22764133

13. Stender CJ, Rust E, Martin PT, Neumann EE, Brown RJ, Lujan TJ. Modeling the effect of collagen fibril alignment on ligament mechanical behavior. Biomech Model Mechanobiol. 2018;17(2):543–557. pmid:29177933

14. Van Zuijlen PPM, Ruurda JJB, Van Veen HA, Van Marle J, Van Trier AJM, Groenevelt F, et al. Collagen morphology in human skin and scar tissue: No adaptations in response to mechanical loading at joints. Burns. 2003. pmid:12880721

15. Wu S, Li H, Yang H, Zhang X, Li Z, Xu S. Quantitative analysis on collagen morphology in aging skin based on multiphoton microscopy. J Biomed Opt. 2011.

16. Pourdeyhimi B, Dent R, Davis H. Measuring Fiber Orientation in Nonwovens Part III: Fourier Transform. Text Res J. 1997;67(2):143–151.

17. Ghassemieh E, Acar M, Versteeg H. Microstructural analysis of non-woven fabrics using scanning electron microscopy and image processing. Part 1: Development and verification of the methods. Proc Inst Mech Eng Part L J Mater Des Appl. 2002;216(3):199–207.

18. Ayres CE, Jha BS, Meredith H, Bowman JR, Bowlin GL, Henderson SC, et al. Measuring fiber alignment in electrospun scaffolds: A user’s guide to the 2D fast Fourier transform approach. J Biomater Sci Polym Ed. 2008;19(5):603–621. pmid:18419940

19. D’Amore A, Stella JA, Wagner WR, Sacks MS. Characterization of the complete fiber network topology of planar fibrous tissues and scaffolds. Biomaterials. 2010;31(20):5345–54. pmid:20398930

20. Liu C, Zhu C, Li J, Zhou P, Chen M, Yang H, et al. The effect of the fibre orientation of electrospun scaffolds on the matrix production of rabbit annulus fibrosus-derived stem cells. Bone Res. 2015;3(January).

21. Redon C, Chermant L, Chermant JL, Coster M. Assessment of fibre orientation in reinforced concrete using Fourier image transform. J Microsc. 1998;191(3):258–265. pmid:9767490

22. Stender CJ, Rust E, Martin PT, Neumann EE, Brown RJ, Lujan TJ. Modeling the effect of collagen fibril alignment on ligament mechanical behavior. Biomech Model Mechanobiol. 2018;17(2):543–557. pmid:29177933

23. Gasser TC, Ogden RW, Holzapfel GA. Hyperelastic modelling of arterial layers with distributed collagen fibre orientations. J R Soc Interface. 2006;3(6):15–35. pmid:16849214

24. Fan R, Sacks MS. Simulation of planar soft tissues using a structural constitutive model: Finite element implementation and validation. J Biomech. 2014. pmid:24746842

25. Wu H, Fan J, Chu CC, Wu J. Electrospinning of small diameter 3-D nanofibrous tubular scaffolds with controllable nanofiber orientations for vascular grafts. J Mater Sci Mater Med. 2010;21(12):3207–3215. pmid:20890639

26. Frahs SM, Oxford JT, Neumann EE, Brown RJ, Keller-Peck CR, Pu X, et al. Extracellular Matrix Expression and Production in Fibroblast-Collagen Gels: Towards an In Vitro Model for Ligament Wound Healing. Ann Biomed Eng. 2018;46(11):1882–1895. pmid:29873012

27. Lau TY, Ambekar R, Toussaint KC. Quantification of collagen fiber organization using three-dimensional Fourier transform-second-harmonic generation imaging. Opt Express. 2012;20(19):21821. pmid:23037302

28. Lutz V, Sattler M, Gallinat S, Wenck H, Poertner R, Fischer F. Impact of collagen crosslinking on the second harmonic generation signal and the fluorescence lifetime of collagen autofluorescence. Ski Res Technol. 2012;18(2):168–179.

29. Lutz V, Sattler M, Gallinat S, Wenck H, Poertner R, Fischer F. Characterization of fibrillar collagen types using multi-dimensional multiphoton laser scanning microscopy. Int J Cosmet Sci. 2012;34(2):209–215. pmid:22235828

30. Annaidh AN, Karine Bruyère, Destrade M, Gilchrist MD, Maurini C, Otténio M, et al. Automated estimation of collagen fibre dispersion in the dermis and its contribution to the anisotropic behaviour of skin. Ann Biomed Eng. 2012;40(8):1666–1678.

31. Bredfeldt JS, Liu Y, Pehlke CA, Conklin MW, Szulczewski JM, Inman DR, et al. Computational segmentation of collagen fibers from second-harmonic generation images of breast cancer. J Biomed Opt. 2014;19(1):016007.

32. Kim A, Lakshman N, Petroll WM. Quantitative assessment of local collagen matrix remodeling in 3-D Culture: The role of Rho kinase. Exp Cell Res. 2006.

33. Sander EA, Barocas VH. Comparison of 2D fiber network orientation measurement methods. J Biomed Mater Res—Part A. 2009;88(2):322–331.

34. Morrill EE, Tulepbergenov AN, Stender CJ, Lamichhane R, Brown RJ, Lujan TJ. A validated software application to measure fiber organization in soft tissue. Biomech Model Mechanobiol. 2016;15(6):1467–1478. pmid:26946162

35. Yang W, Sherman VR, Gludovatz B, Schaible E, Stewart P, Ritchie RO, et al. On the tear resistance of skin. Nat Commun. 2015;6:6649. pmid:25812485

36. Polzer S, Gasser TC, Forsell C, Druckmüllerova H, Tichy M, Staffa R, et al. Automatic identification and validation of planar collagen organization in the aorta wall with application to abdominal aortic aneurysm. Microsc Microanal. 2013;19(6):1395–1404. pmid:24016340

37. Schriefl AJ, Wolinski H, Regitnig P, Kohlwein SD, Holzapfel GA. An automated approach for three-dimensional quantification of fibrillar structures in optically cleared soft biological tissues. J R Soc Interface. 2012.

38. Puschmann S, Rahn CD, Wenck H, Gallinat S, Fischer F. In vivo quantification of human dermal skin aging using SHG and autofluorescence. Multimodal Biomed Imaging VII. 2012;8216:821608.

39. MATLAB. version 9.4.0.813654 (R2018a). The MathWorks Inc.; 2018.

40. MATLAB. Image Processing Toolbox (version 10.2). The MathWorks Inc.; 2018.

41. MATLAB. Curve Fitting Toolbox (version 3.5.7). The MathWorks Inc.; 2018.

42. Paul H, Jex I, Paul H, Jex I. Photon statistics. Introd to Quantum Opt. 2010; p. 127–154.

43. Becker RI, Morrison N. The errors in FFT estimation. IEEE Trans Signal Process. 2002;44(8):133–135.

44. Withayachumnankul W, Fischer BM, Lin H, Abbott D. Uncertainty in terahertz time-domain spectroscopy measurement. J Opt Soc Am B. 2008;25(6):1059.

45. Moisan L. Periodic plus smooth image decomposition. J Math Imaging Vis. 2011;39(2):161–179.

46. University of Wisconsin. Public-Domain Test Images for Homeworks and Projects;. Available from: https://homepages.cae.wisc.edu/~ece533/images/.

47. Lagarias C J, Reeds A J, Wright H M, Wright E P. Convergence Properties of the Nelder-Mead Simplex Method in Low Dimensions. SIAM J Optim. 1998;9(1):112–147.

48. Koenig K, Riemann I. High-resolution multiphoton tomography of human skin with subcellular spatial resolution and picosecond time resolution. J Biomed Opt. 2003;8(3):432.

49. Rahn CD, Meine H, Gallinat S, Wenck H, Wittern KP, Fischer F. Fully automated data acquisition and fast data interpretation in a customized multimodal multiphoton microscope. In: Three-Dimensional Multidimens. Microsc. Image Acquis. Process. XVII; 2010.

50. Schwarz M, Riemann I, Stracke F, Huck V, Gorzelanny C, Schneider SW, et al. A comparative study of different instrumental concepts for spectrally and lifetime-resolved multiphoton intravital tomography (5D-IVT) in dermatological applications. In: Imaging, Manip. Anal. Biomol. Cells, Tissues VIII; 2010.