Fourier transform-based method for quantifying the three-dimensional orientation distribution of fibrous units

Several materials and tissues are characterized by a microstructure composed of fibrous units embedded in a ground matrix. In this paper, a novel three-dimensional (3D) Fourier transform-based method for quantifying the distribution of fiber orientations is presented. The method allows for an accurate identification of individual fiber families, their in-plane and out-of-plane dispersion, and showed fast computation times. We validated the method using artificially generated 3D images, in terms of fiber dispersion by considering the error between the standard deviation of the reconstructed and the prescribed distributions of the artificial fibers. In addition, we considered the measured mean orientation angles of the fibers and validated the robustness using a measure of fiber density. Finally, the method is employed to reconstruct a full 3D view of the distribution of collagen fiber orientations based on in vitro second harmonic generation microscopy of collagen fibers in human and mouse skin. The dispersion parameters of the reconstructed fiber network can be used to inform mechanical models of soft fiber-reinforced materials and biological tissues that account for non-symmetrical fiber dispersion.

Biphasic solids, characterized by a microstructure of fibrous units embedded in a ground substance, are common in both engineered materials and biological tissues.For example, fiber-reinforced materials have emerged as versatile and indispensable components in a wide range of structural engineering applications, ranging from aerospace and automotive to construction, including fiber-reinforced polymers 1-4 , fiber-reinforced concrete 5,6 and non-woven fabrics 7 .Their importance goes beyond traditional applications and finds increasing application in the medical sciences and for smart materials [8][9][10] .
In addition, researchers have thoroughly investigated the mechanical behavior of biological tissues with a fibrous microstructure, including cornea, cartilage, skin, and blood vessels.Collagen is of particular interest, it is also the most abundant protein in the human body and plays a role in providing structural integrity and load-bearing functions of tissues.Understanding the distribution and organization of collagen fibers has significant implications for the preliminary diagnosis of pathological situations.For example, alterations in collagen organization have been associated with diseases such as cancer [11][12][13] , abdominal aortic aneurysms [14][15][16] , mitral valve disease 17 , keratoconus [18][19][20][21] , and genetic disorders such as Marfan syndrome 22 .Quantification of the underlying microstructure of soft tissues is also of paramount importance for constitutive models developed for fibrous biological tissues.Accurate models rely on microstructural parameters, such as the mean fiber orientation and dispersion, to accurately reproduce the nonlinear anisotropic response of tissues under complex loading scenarios [23][24][25][26] .
Precise imaging techniques play a crucial role in assessing the architecture of collagen fibers.Several methods, including histological staining 27,28 , polarized light microscopy 29,30 , optical coherence microscopy 31,32 , and second harmonic generation (SHG) microscopy [33][34][35] are used to visualize collagen fibers.Among these techniques, SHG microscopy has gained prominence due to its ability to provide three-dimensional, high-resolution images of collagen fibers that can be detected at depth without the need for staining.However, the quantitative threedimensional (3D) assessment of the fiber orientation distribution (FOD) remains a challenge.

www.nature.com/scientificreports/
The importance of automated techniques to measure fiber distribution cannot be overlooked.These methods have the potential to revolutionize the analysis of collagen organization by reducing human error, increasing efficiency and enabling the measurement of a greater number of features.Current approaches to quantitatively measure fiber distribution include texture analysis 36,37 , fiber segmentation and tracking algorithms (e.g., the CT-FIRE algorithm) 38,39 , pixel-by-pixel methods using gradients 40,41 or directional variance 42,43 , Gabor filter methods 28 , and Fourier transform methods 32,[44][45][46] .In particular, Fourier transform-based methods, including wavelet transforms, have shown promise in terms of computational efficiency irrespective of image complexity, computation speed, and reduced sensitivity to noise and fiber crimping 32 .
While numerous methods exist for two-dimensional (2D) measurements of fiber distribution, the transition to 3D analyses raises open technical questions.For example, fibers with a large elevation with respect to the image plane are not captured because their cross-section is quasi-circular and the direction cannot be recognized by the 2D Fourier transform 47 .Also, existing 2D algorithms either focus exclusively on in-plane or out-of-plane measurements, which limits the comprehensive understanding of fiber dispersion.Even if both the in-plane and out-of-plane fiber distributions are measured, the architecture cannot be combined to produce a 3D distribution since it is not possible to reconstruct the 3D distribution from the distribution about two perpendicular directions without information about their covariance.
Few algorithms are explicitly designed for full 3D measurements, leaving a significant gap in the field 48 .Liu et al. 13 used a pixel-by-pixel method to extract a spatial description of collagen fibers, but then adopted 3D directional variance, an averaged metric to assess fiber concentration levels, and they provided separate in-plane and out-of-plane fiber distributions.Lau et al. 47 used a 3D Fourier transform method, but the fiber distribution is computed by dividing the 3D image into regions of interests (ROIs) and finding the overall fiber orientation of each region.The fiber distribution within each region is not considered, so information is lost when a region contains two or more fiber families.In addition, the number of measurements is limited to the number of regions of interest.
In this work, we present a novel 3D Fourier transform-based method to characterize the FOD.Our approach provides a complete 3D description of fiber dispersion within an image, facilitates the identification of individual fiber families, and enables quantitative analysis using fitted bivariate von Mises probability density function (PDF) parameters.The algorithm has fast computation times even with large-sized images and is therefore extremely practical for research and clinical applications.To validate the algorithm, 3D images containing artificial fibers are generated to assess its robustness and precision in computing in-plane and out-of-plane distribution parameters.The performance of the algorithm is then verified for mouse and human skin collagen fibers obtained from SHG images.

Discrete fiber orientation distribution
Fibers in biphasic materials are usually dispersed with a non-uniform orientation distribution characterized by one or more fiber families, each having a preferential direction and generally a non-symmetrical dispersion with different in-plane and out-of-plane concentrations.For example, biological tissues such as skin or arterial walls exhibit two distinct fiber families, together with an out-of-plane concentration that is more pronounced compared to the in-plane concentration 49 .In contrast, in 3D printed fiber-reinforced hydrogels, the fibers are dispersed axially along a single preferential direction, while the concentration mainly depends on the diameter of the deposition nozzle and the printing speed 3,9 .
The three-dimensional fiber organization can be described by a FOD function that returns the normalized amount of fibers along all directions in the angular domain of a hemisphere, where θ and φ are the azimuthal and elevation angles with respect to the in-plane and out-of-plane directions of a single fiber.The FOD reveals the architecture of the fiber network: the location of a local maximum in the FOD identifies the mean fiber direction, while the standard deviations σ θθ and σ φφ provide a measure of the in-plane and out-of-plane concentrations.A lower standard deviation corresponds to a higher fiber alignment along the mean direction.
We propose an algorithm that takes advantage of the directional information carried by the 3D discrete Fourier transform to obtain a discrete fiber orientation distribution (dFOD).First, the spectrum of a threedimensional grayscale image of the fiber network is filtered with a set of filters to obtain a raw dFOD d(θ, φ) (Materials and Methods, Eq. ( 6)).Then, we remove the interference due to the reciprocal overlaps of the filters in the frequency domain with a deconvolution step to obtain the actual (deconvoluted) dFOD d ′ (θ, φ) .Finally, we estimate the parameters of the fiber orientation distribution by fitting a bivariate von Mises PDF to the deconvoluted dFOD (Materials and Methods, Eqs. ( 2) and ( 3)).For multiple fiber families, the parameters of each n-th family are obtained by fitting a combination of bivariate von Mises PDFs ρ = n ν f ,n ρ n , where ν f ,n is the volume fraction of the family.The parameters of the distribution include the azimuthal and elevation angles α n and β n of the mean fiber direction, the rolling angle γ n of the family about the mean direction, and the in-plane and out-of-plane concentration parameters a n , b n (see Fig. 9 in Materials and Methods).
Figure 1 illustrates an application of the algorithm to an artificially generated fiber network.A fiber dispersion image ( 256 × 256 × 256 voxels) with N = 1 000 randomly generated fibers of diameter t = 3 voxels was gener- ated according to a prescribed von Mises distribution ( n = 1 ) with mean orientation angles α = β = γ = 0 • and concentrations a = 0.5 , b = 5 .The raw and deconvoluted dFODs are shown in Fig. 1c,d.In particular, the latter show that the algorithm is able to precisely recover the prescribed von Mises distribution ρ ( R 2 = 0.981).

Algorithm calibration and validation
We employed artificial fiber stacks with prescribed distributions to calibrate the algorithm and assess its robustness and precision (see Materials and Methods).

Calibration
The shape of the deconvoluted dFOD d ′ (θ, φ) is controlled by a power parameter q applied to the spectrum of the fiber image before the filtering process (Eq.( 6)).Since the spectrum power parameter affects the standard deviations of d ′ (θ, φ) , we analyzed the relative errors �σ θθ and �σ φφ between the standard deviations of the prescribed and the measured distributions.For this purpose, we generated four different single-fiber family distributions representing four scenarios of high/low in-plane and out-of-plane concentrations: Case 1: low concentration in both planes ( a = 0.5 , b = 0.5 ); Case 2: high out-of-plane concentration ( a = 0.5 , b = 5 ); Case 3: high in-plane concentration ( a = 5 , b = 0.5 ); Case 4: high concentration in both planes ( a = 5 , b = 5 ).In total we assumed α = β = γ = 0 • and generated 10 three-dimensional images ( 256 × 256 × 256 voxels) each with N = 6 000 fibers of diameter t = 3 voxels.
Figure 2a shows contour plots of the prescribed von Mises PDFs.The relative error �σ θθ is represented in Fig. 2b.For both Case 1 and 2, the error seems small and increases slightly with q, always staying below 6% .In Case 3, on the other hand, the error drops quickly from 50% to less than 10% without significant variations between q = 2.4 and q = 2.8 .In Case 4, the error decreases significantly with increasing values of the spectrum power parameter q and reaches its minimum at q = 2.8 .Regarding the relative error �σ φφ in the elevation angle, the results are shown in Fig. 2c.In Cases 1 and 3, the error increases with q and always stays below 5% .In Case 2, the errors are quite high, except for q = 2.4 , as in Case 4, where we observed a decreasing trend for q.We also found that errors are higher when concentrations are high along each direction.This is because the standard deviation is smaller at high concentrations and the relative error is more sensitive to small differences in the measured standard deviation since the standard deviation of the prescribed distribution tends to zero.

Robustness
We tested the ability of the algorithm to accurately measure dFOD, independently of the fiber number N and diameter t in the 3D image.The results are also compared in terms of fiber density δ , defined as the average of the voxel intensity over the entire 3D image volume, where 0% corresponds to an empty volume and 100% to completely bright voxels.We generated three-dimensional images ( 256 × 256 × 256 voxels) with increasing numbers of fibers and same diameter t = 3 voxels.Representative distributions with parameters a = 0.5 , b = 5 , α = β = γ = 0 • ( n = 1 ) were generated containing N = 1 000 , 2 000 , 5 000 , 10 000 , 20 000 , 50 000 , 75 000 , 100 000 fibers, corresponding to fiber densities between δ = 0.72% and δ = 28.23%(Fig. 3a).For each case 10 three-dimensional images are generated.
As above, we assessed the algorithm performance by computing the relative errors in the standard deviations along the azimuthal and elevation directions (Fig. 3b).It was found that the error �σ θθ depends slightly on the fiber number N and generally reduces at higher N.This trend might change at higher N because for N → ∞ ( δ → 100% ) fibers cannot be distinguished and no distribution can be measured.However, a volume with N = 100 000 fibers ( δ = 28.23% ) represents a density beyond that observed in real collagen fiber tomographies, as shown in the following applications to skin tissue ( δ = 16.97% for human skin, δ = 16.31% for mouse skin).With regard to the elevation angle, the error �σ φφ appears to be practically insensitive to N and shows small values, except for N = 100 000.
To assess the performance of the algorithm in terms of the fibers diameter, we analyzed artificial fiber stacks with the same number of fibers N = 2 000 and increasing fiber diameters of t = 3 , 5, 7, 9, 11 and 13 voxels, corresponding to densities from δ = 1.28% to 64.80% .For each diameter we generated 10 three-dimensional images ( 256 × 256 × 256 voxels) with same dispersion parameters, a = 0.5 , b = 5 , α = β = γ = 0 • ( n = 1 ).Representative images for each diameter are shown in Fig. 4a.As illustrated in Fig. 4b, the algorithm appears to be precise and insensitive to the fiber diameter along the azimuthal direction and predicts the true standard deviation with an error �σ θθ smaller than 3.5% , with no significant variations.On the other hand, the error along the elevation direction is small only for fiber diameters between 3 and 7 voxels ( �σ φφ < 8.7% ) and increases thereafter to about 50% for t = 13 voxels.
Interestingly, the behavior observed in these two analyses is similar and is directly related to fiber density.In fact, by comparing the errors in Figs.3b and 4b with respect to the density δ , it follows that �σ θθ and �σ φφ are correlated.The errors �σ θθ are small in in both cases and almost insensitive to the fiber density, while �σ φφ is small only for densities lower than δ ≈ 25% ( N = 75 000 , t = 3 voxels; N = 2000 , t = 7 voxels) and increases significantly for densities δ > 28% ( N = 100 000 , t = 3 voxels; N = 2000 , t = 9 voxels).This suggests that the algorithm only responds to fiber density and not individually to the number of fibers or their diameter.However, the difference between �σ θθ and �σ φφ is probably due to the small prescribed standard deviation σ φφ (higher concentration) along the elevation direction, increasing the relative error even though the absolute error is small.www.nature.com/scientificreports/

Precision
The precision of the algorithm is measured by its ability to estimate the parameters of the in-plane and out-ofplane distributions.For each of the combinations of parameters described below, a series of 10 three-dimensional images ( 256 × 256 × 256 voxels) was generated with N = 6 000 fibers of t = 3 voxel diameter.
For the in-plane parameters, we explored all combinations of mean angles α = 0 • , 30 • , 60 • , 90 • with con- centrations a = 0 , 0.5, 1, 2, 5, while the out-of-plane angle β = 0 • and concentration b = 1 remained fixed.For reasons of symmetry, negative angles are not considered.Errors between the estimated and true parameters are evaluated by the difference for the in-plane concentration and the absolute difference for the mean in-plane angles.The results are summarized in Fig. 5.
The estimate of the true concentration a t seems almost insensitive to the true mean fiber orientation α t , except for a t = 5 where the algorithm underestimates the true concentration with an average error of −0.745 (Fig. 5a).In absolute terms, the error increases with increasing a t .The opposite behavior can be observed for the estimation error of the mean angle (Fig. 5b).The isotropic case a t = 0 is not shown because the mean angle in not relevant in this case.Thereby, the precision improves with concentration a t since the uncertainty on the Letters must be compared among the same cases.www.nature.com/scientificreports/angular peak location reduces for more concentrated dFODs.Errors are independent of the true mean fiber orientation α t , show no significant variations (except for three means) and remain small for all values of a t and α t , with a maximum of 1.8 • .
For the out-of-plane precision assessment, we analyzed all combinations of mean angles β = 0 • , 30 • , 60 • , 90 • with concentrations b = 0 , 0.5, 1, 2, 5 while the in-plane angle α = 0 • and concentration a = 1 remained fixed.For reasons of symmetry, we have not considered negative angles.The results regarding the errors between the estimated and the true parameters are summarized in Fig. 6.
The estimated out-of-plane concentration parameter (Fig. 6a) shows a moderate dependence on the angle β t , with low values and no significant differences between 0 • and 60 • for concentrations b t ≤ 2 .However, for b t = 5 the concentration is underestimated at lower angles, but for β t = 60 • the error tends to decrease and for β t = 90 • it becomes positive.Regarding the mean out-of-plane angle (Fig. 6b), the error shows almost no dependence on β t , with significant variations between the different concentrations only at β t = 90 • .Similar to the in-plane case, the errors are generally low with a maximum of 1.65 • and tend to decrease with increasing concentration b t .

Application to skin tomography
We applied the algorithm to compute the dFOD and relative distribution parameters from SHG tomography of human and mouse skin tissue (see Materials and Methods).

Human skin
Figure 7a shows a tomography of human skin collagen fibers constructed from a sequence of SHG images acquired in a volume of 465 µ m × 465 µ m × 116.25 µ m (fiber density δ = 16.97% ).A detailed animated view of the tomography is shown in the Supplementary Video 1 available online.In contrast to the artificial straight fibers considered in the previous sections, collagen fibers in the skin are wavy and compacted in bundles of different sizes and diameters 50 .The bundles are in turn interwoven, making it difficult to distinguish them within the volume.However, this does not affect the accuracy of the measured dFOD as the algorithm can detect any fibrous unit within the specified diameters.In addition, wavy fibers can be viewed as a continuous sequence of smaller straight fibers.
The global dFOD, reported in Fig. 7b, is obtained by dividing the volume into 16 cubic ROIs with a size of M = 256 voxels ( 116.25 µm), as specified in the Materials and Methods section, and summing the respective 16 dFODs.The distribution has an angular resolution of 5.8 • ( D = 31 ) and is normalized such that the volume under the surface is one.As expected from experimental evidence 51 , the concentration along φ is higher than along θ , which means that almost all fibers lie in a co-planar plane x-y with the skin surface.Although the distribution shows only one main peak at around θ = −60 • and φ = 0 • , qualitatively indicating the presence of a single fiber family, a combination of two bivariate von Mises PDFs ( n = 2 ) (wireframe in Fig. 7b) is adopted to improve the quality of the fit.Note that the algorithm may underestimate the two out-of-plane concentrations b i , because both values are outside the tested range in the out-of-plane precision analyses.
In Fig. 7c, the planar distribution computed with the 3D discrete Fourier transform-based algorithm is compared with that obtained using a classic 2D discrete Fourier transform-based approach 32,44,45,52 .With our algorithm (gray histogram), the planar distribution is obtained by summing the dFOD along the direction φ shown in Fig. 7b, while the 2D case (red histogram) is realized by adding the individual distributions of all in-plane x-y slices of the tomography, each derived with the 2D algorithm (angular resolution of 1 • ).The two distributions are qualitatively similar with both peaks located around θ = −55 • and showing comparable con- centrations.This observation is also confirmed by the parameter set obtained from fitting the 2D von Mises distribution (red solid curve).Means not sharing uppercase letters differ significantly by the Tukey-test at the 5% significance level.Letters must be compared within the same concentration b t .

Mouse skin
Figure 8a shows a tomography of mouse skin collagen fibers constructed from a sequence of SHG images acquired in a volume of 465 µ m × 465 µ m × 116.25 µ m (fiber density δ = 16.31% ).Collagen fibers in mouse skin are finer compared to human fibers and warp around hair follicles, appearing as cavities in the tomography.A detailed animated view of the tomography is shown in the Supplementary Video 2 available online.
As for human skin, the global dFOD is calculated by dividing the volume into 16 cubic ROIs with a size of M = 256 voxels ( 116.25 µ m) and summarizing the respective 16 dFODs (schemes of the ROIs are omitted here).The distribution shown in Fig. 8b appears to be more clustered along the horizontal plane x-y with respect to the human skin, showing two distinct peaks at about θ = −65 • and θ = 60 • ( φ = 0 • ).A combination of two bivariate von Mises PDFs ( n = 2 ) is adopted for the fitting (wireframe in Fig. 8b).
The comparison between the planar distributions obtained with the 2D and 3D discrete Fourier transformbased algorithms is reported in Fig. 8c.Here the two peaks observed in the distribution constructed from the 3D dFOD in Fig. 8b are also observed in the 2D discrete Fourier transform-based distribution, although the trend seems smoother.In this case, the in-plane parameters present noticeable differences between the two algorithms.This is likely due to the of the 2D algorithm to detect fibers with high elevation, where quasi-circular fiber cross-section is detected as isotropicly distributed in the x-y plane.With the 2D approach, the concentration parameters are lower and the volume fractions appear more unbalanced.In addition, the mean in-plane fiber directions are closer to the x direction ( θ = 0 • ).However, despite the noticeable differences, the two collagen families are well captured by both algorithms.

Discussion
We proposed a discrete Fourier transform-based algorithm that extends the classical 2D discrete Fourier transform-based methods of measuring fiber orientation to the 3D case 45,49,53,54 .The filter used to identify the orientation of fibers in the 3D space is obtained by rotating a 2D band-pass filter (wedge filter) about its axis.However, in the 3D case, an extra step is required before finding the actual dFOD d ′ (θ, φ) .In particular, a deconvolution process is applied to clean each directional measurement of the interference from all other directions.Such interference is caused by the intersections of the 3D filters pivoting around the fixed frequency point (u, v, w) = (M/2, M/2, M/2) in the frequency domain, where M is the size of the 3D cubic image.Unlike the 2D case, where wedge filters can be shaped without overlap, filters in a 3D framework always intersect, even in case they degenerate to discs.Therefore, deconvolution is a necessary step to mathematically obtain a faithful representation of the actual fiber orientation distribution.
The algorithm is based on three parameters that need to be adjusted: the cut-off frequencies f min and f max and a spectrum power parameter q (see Materials and Methods).The two frequencies are derived directly from the examined maximum and minimum fiber diameters, respectively, while the parameter q is calibrated to match the output dFOD with the true distribution of artificial fiber stacks.This parameter behaves similarly to the one proposed by Polzer et al. 45 , with the difference that in their case the entire distribution d(θ) is raised instead of the spectrum |ĝ| .The three-parameter method gives satisfying results, but more sophisticated filtering techniques developed for 2D algorithms can also be integrated into our method.For example, Witte et al. 46 proposed an adaptive filtering technique based on the propagation of uncertainties that excludes spectrum magnitudes not carrying fiber-related information, thus, eliminating the need to define specific cut-off frequencies based on fiber diameters.
Given the fundamental importance of q in the proposed method, we performed a calibration analysis to assess its influence on the concentration parameters a and b of the von Mises distribution.Ideally, the errors �σ θθ and �σ φφ of the measured standard deviations should be independent on the specific dispersion analyzed and dependent only on q.In other words, there should be a unique value of q for which the errors �σ θθ and �σ φφ are minimized simultaneously and independently of the true fiber distribution.In practice, however, the errors also depend on the specific fiber dispersion considered.From Fig. 2b,c we can deduce that for increasing concentrations the algorithm tends to produce a smoother dFOD than the true one.This suggests that a larger spectrum power parameter q might be necessary to increase the sharpness of the measured distribution.
This behavior is also reflected in the Figs.5a and 6a, in which the true concentrations are underestimated for high concentrations ( a t = 5 , b t = 5 ).Their precision would improve if the value of q were increased, but on the other hand the algorithm would overestimate at lower concentrations.Since our goal was to get accurate estimates of the fiber dispersion parameters in most of the cases, we identified q = 2.4 as an acceptable compro- mise.In contrast to the concentration parameters, the true mean orientation angles α t and β t depend only on the distribution peaks and are not directly affected by q.As shown in Figs.5b and 6b, their precision generally increases with concentration as the uncertainty about the peak location is reduced for sharper distributions.Note that the precision of the parameters depends on the algorithm, while the nonlinear least squares fit ensures a reliable estimate of the dispersion parameters.The robustness analysis revealed that the reliability of our method is not affected by the number of fibers N and their diameter t, although they limiting factors may become larger in the out-of-plane measurements when fiber density δ is greater than about 25% .In order to understand the capabilities of the proposed method in detail, future analyses should deal with, e.g., the influence of noise 46 and the influence of the fiber waviness.
Despite the approximations and limitations discussed, the algorithm provides accurate quantitative descriptions of the spatial arrangement of complex fiber dispersions.In contrast to 2D approaches, where fibers that are strongly inclined from the x-y plane are difficult to detect due to their small cross-section, our method can consider all fibers with any orientation in 3D space.Differently from existing 3D methods 13,47 , which provide separate distributions along θ and φ , our algorithm provides a complete description of the spatial fiber distribu- tion over the entire 3D angular domain.This allows to identify the different fiber families in the volume and to compute for each family the specific set of parameters without ambiguities.In fact, when multiple fiber families are present, the different peaks cannot be combined together if only the separate in-plane and out-of-plane distributions are known, without information on how they co-vary in the two-dimensional angular domain.
In practice, as shown for human and mouse skin, this scenario rarely occurs in biological soft tissues where fibers are mainly aligned in the x-y plane (Figs.7b and 8b).In these cases, the fiber families and their in-plane parameters can be computed from the in-plane distribution, while the out-of-plane parameters can be extracted from the 2D distribution of a vertical tissue section.However, this is still an approximation since all the fiber families are assigned a unique concentration b 23 .Another benefit of full 3D distributions is the ability to measure the rolling angle γ , which allows for general rotations of the fiber family in space and compensates for slight planar misalignments of the samples during 3D imaging.
In addition, our method is fast and robust.To analyze stacks of 1024 × 1024 × 256 voxels ( 465 µm× 465 µ m × 116.25 µ m) the overall computation took 65 s for human skin (60.1 s for the computation of the raw dFOD and 4.9 s for the deconvolution), and 61 s for mouse skin (56.2 s for the computation of the raw dFOD and 5.1 s for the deconvolution, on a 2.20GHz CPU, 32GB RAM Desktop PC).Of course, the computation speed depends on the resolution of the dFOD, since the number of computations increases proportionally with the angular interval D 2 .Also, memory can be a limiting factor, since the free memory needed to store the filters increases proportionally to D 2 × M 3 .However, with a value of D = 31 adopted in the provided examples, distributions with a good angular resolution can be obtained in a reasonable computation time.
The robustness analysis showed that the algorithm reliably predicts the same FOD for stacks with a fiber density of up to 25% .This represents a major advantage over pixel-by-pixel methods, which are affected by ori- entation uncertainties in the pixels when two or more fibers overlap in very dense fiber dispersions 43 .In general, approaches based on discrete Fourier transforms are less affected by image complexity since fiber orientations are computed from the frequency domain, which is insensitive to fiber overlaps in the spatial domain 32 .It is important to note that the proposed method requires high-quality stacks with a z resolution comparable to that in the x-y plane, which are typically expensive and time-consuming to obtain.The SHG stack of human and mouse skin collagen in Figs.7a and 8a (shown after the vertical resampling) is acquired with our equipment at the highest possible z resolution ( �z = 0.57 µm), each requiring up to 2h to complete the acquisition.In addition, to reduce the blurring effect along the z direction due to the diffraction of the confocal microscope, we pre-processed the tomographies by applying three-dimensional deconvolution.This process is slow and computationally expensive, but recent deep learning techniques have shown promising results in deconvolution of microscope images 55 as well as in reducing noise 56 .Furthermore, since the algorithm is based on grayscale 3D images, any technique can be used to capture the stacks, making our method applicable to various imaging technologies, including x-ray computational tomography (xCT) 57 , optical coherence tomography (OCT), and ultrasound imaging 58,59 .
In summary, the presented novel algorithm is able to reliably quantify the fiber network obtained from 3D images by determining the dispersion parameters of each individual fiber family.These parameters can be used to inform mechanical models of soft fiber-reinforced materials and biological tissues that account for nonsymmetrical fiber dispersion.

Orientation and distribution of fibers
A general fiber orientation distribution can be assumed as a combination of a finite number of fiber families, each of which is described by a PDF ρ(N) providing the normalized angular density of the fibers in the direction of the unit vector N .We introduce two different bases: a general basis {E a } a=1,2,3 , and a principal basis of the fibers {M , M ip , M op } , where M represents the mean fiber direction, (M, M ip ) is the plane containing the mean direction of the fiber family and M op is the out-of-plane normal (Fig. 9).With respect to the principal frame, the fiber unit vector N is written as where and are the azimuthal and elevation angles of the fiber unit vector N with respect to the principal frame (Fig. 9a).
The PDF can be decomposed as a bivariate distribution of the form ρ(N) = ρ(�, �) = ρ ip (�)ρ op (�) , where ρ ip and ρ op represent the in-plane and out-of-plane distribution functions, respectively.In terms of the global frame, the general expression of the PDF is as follows where θ and φ represent the azimuthal and elevation angles of the unit vector N(θ, φ) in the global frame (Fig. 9b).
We consider a rigid rotation of the fiber families in Euclidean space by introducing the Tait-Bryan angles α , β , γ .The first two angles denote the azimuthal and elevation angles of the mean direction M and γ the rotation of the principal fiber frame by M (Fig. 9c).In general, the functions �(θ, φ) and �(θ, φ) in Eq. ( 2) are implicit and depend on the three angles α , β , γ .Explicit expressions are only available for some simple cases, such as β = γ = 0 for which � = θ − α and � = φ.
The functions ρ ip (�) and ρ op (�) are provided by two π-periodic von Mises distributions of the form 23 where a and b are the in-plane and out-of-plane concentration parameters, respectively, I 0 (x) = 1 π π 0 exp(x cos t)dt is the modified Bessel function of the first kind of order zero, and erf is the error function given by erf Large values of a and b correspond to a fiber distribution with a high degree of alignment along the mean fiber direction M , while for a, b → 0 an isotropic distribution is obtained.This choice for the in-plane and out-of-plane functions satisfies the symmetry requirement ρ(N) = ρ(−N) , or equivalently ρ(θ, φ) = ρ(θ + 180 • , −φ) , and the normalization condition over the unit sphere S 2 , 1 4π S 2 ρ(N)d� = 1 , required in mechanical models 23,24 .

3D discrete Fourier transform-based algorithm for measuring fiber orientation distributions
We have developed a novel method based on the 3D discrete Fourier transform in combination with a so-called funnel filtering to extract directional data from three-dimensional images.The algorithm is implemented in a custom Matlab code.

Raw discrete fiber orientation distribution
Let g(x, y, z), g : N 3 → R , be the signal in the discrete spatial domain representing a M × M × M image (voxel volume).We assume that the x, y, z axes are aligned with the global basis vectors E 1 , E 2 , E 3 , so that the definition of the azimuthal and elevation angles θ and φ is the same in both frames.If g(x, y, z) depicts a bundle of straight fibers all oriented along a certain direction (θ, φ) (Fig. 10a), then its spectrum ĝ(u, v, w) , ĝ : where the dominant magnitudes |ĝ(u, v, w)| are scattered in a plane that pass through (u, v, w) = (M/2, M/2, M/2) and perpendicular to the frequency direction (θ, φ) (yellow voxels in Fig. 10b).Note that the values of the spec- trum are shifted such that the lower frequencies occupy the central position, with frequency f = 0 located at (u, v, w) = (M/2, M/2, M/2) .Gibbs artifacts characterized by high magnitudes distributed along the u, v, w axes are removed by tapering g(x, y, z) with a 3D Tukey window with 40% cosine before the transform.Note that the orientation of the plane in the frequency domain is invariant with changes in fiber position in the spatial domain.Therefore, the sum of the moduli |ĝ(u, v, w)| filtered from this plane is a measure of the amount of straight fibers oriented towards (θ, φ) .Due to the aliasing filter plane, however, this measure reacts sensitively to small variations in the sampled direction (θ, φ) .To reduce this sensitivity, we perform the com- putation over a finite solid angle rather than the exact direction, so the filter corresponds to a volume obtained by enveloping all planes of all directions within the solid angle.The volume is still aliased along its boundary (Fig. 10c), but it represents only a lower fraction of the total number of frequencies filtered.Therefore, it is less sensitive to small changes in the sampled direction.To exclude the part of the spectrum that does not contain relevant fiber information filter is reduced to include magnitudes only between minimum and maximum cut-off frequencies f min and f max .Assuming a conical solid angle with the amplitude ε and the mean direction (θ, φ) the funnel-shaped filter is defined by where ũ, ṽ, w R w and R v are the rotation matrices around the w and v axes, respectively (Fig. 10b).
In Eq. ( 6) we introduced an exponent q to correct the distribution.Usually, in 2D discrete Fourier transformbased algorithms 53,54,60 , the exponent is set to q = 2 , but other values can be used to achieve better results 45 .The calibration analysis is used to identify the optimal value of the parameter to match the dFOD to the prescribed distribution of artificial fiber images.To ensure a thorough measurement of all fibers in the discrete angle set, we assume that the cone opening corresponds to the angular resolution ε = �θ = �φ .Assuming an image size of M = 256 voxels, which is suitable for most biological applications, and without assumptions about the investigated fiber diameter, we set the cut-off frequencies to f min = 4 and f max = 43 to exclude high frequencies associated with noise and low frequencies associated with exposition variations within the image 54 .Given the relationship f = M/(2t) , this frequency range corresponds to fiber diameters between t = 32 and t = 3 voxels.However, this does not necessarily mean that fibers larger than 32 voxels and smaller than 3 voxels will not be detected, since the information conveyed by the discrete Fourier transform about a specific fiber diameter t ′ is spread in a band of about f ′ ± 0.1f ′53 .
An essential aspect of the proposed method is that the signal g and accordingly the spectrum ĝ correspond to cubes with the same size M × M × M .Different sizes along the three dimensions of g and ĝ would result in an inaccurate dFOD since the higher number of frequencies along the dominant direction would make a larger contribution with respect to other directions in the Eq.(6).
Another aspect to consider is computational efficiency.For computational reasons, all filters are computed before the analysis, which requires free memory proportional to the product between the number of directions D 2 and the number of voxels in the filter M 3 .This might become too demanding when analyzing high-resolution 3D images ( M ≈ 1024 ).To limit the amount of memory required for the computations while maintaining a cubic shape for ĝ , the 3D image is divided into ROIs of 256 × 256 × 256 voxels.Then the raw dFOD d(θ, φ) is computed as the sum of the raw dFODs of each individual ROI.The sum does not require any weighting, since ĝ already takes into account the luminance, i.e. the overall amount of fibers within the ROI.If the dimensions are not a multiple of 256, the 3D image can be divided into cubic ROIs of different size (possibly close to 256), which can be scaled up or down to M = 256 by volumetric interpolation before the 3D discrete Fourier transform.

Deconvolution of the dFOD
The distribution d(θ, φ) needs to be deconvoluted to obtain the actual dFOD d ′ (θ, φ) .Because the filters overlap in the frequency domain, the measured amount of fibers along a generic direction I = (θ i , φ j ) contains a per- centage of the actual amount of fibers from all other directions J = (θ k , φ l ) that is proportional to the number of shared magnitudes between filters along I and J directions.
Mathematically, we can write the raw dFOD in vector form by the linear relationship d = Kd ′ , where d ′ denotes the unknown deconvoluted dFOD in vector form.The symmetric convolution kernel matrix K IJ rep- resents the normalized number of common magnitudes between filters I and J. Since the distributions can be arbitrarily rearranged from matrix to vector, we use the index rule I = (j − 1)D + i .We obtain the deconvoluted vector d ′ using the iteratively constrained Tikhonov-Miller algorithm, which gives the optimal non-negative solution to Here the second term depends on the discrete Jacobian operator matrix L and the regularization parameter ∈ [0, 1] which was introduced to avoid noise amplification in the solution 61 .This method is required because direct matrix inversion would provide noisy and possibly non-physical (negative values) results.A step-size of 0.1 and = 10 −3 is used for the iterations.Once the solution is found, the deconvoluted vector d ′ is rearranged again as a matrix using the aforementioned reordering rule, providing the sought dFOD d ′ (θ, φ).

Fiber orientation distribution parameters
To estimate the parameters of the fiber orientation distribution, the function ρ(θ, φ) = ρ(θ, φ) cos(φ) is fitted to the deconvoluted dFOD d ′ (θ, φ) , where ρ(θ, φ) is defined in Eqs. ( 2) and (3).The cosine of the elevation angle φ accounts for the surface area d on the unit sphere ( d� = cos φdθdφ ).Note that d ′ (θ, φ) represents the normalized quantity of fibers along each discrete direction, while the PDF ρ(θ, φ) n needs to be multiplied by the surface element area to obtain a normalized fiber quantity.For the multi-fiber distribution, a linear combination ρ = n ν f ,n ρ n is used, where ρ n is the function with respect to the n-th fiber family and ν f ,n its volume fraction, with n ν f ,n = 1.
The operation is performed in Matlab using the integrated nonlinear least squares function lsqnonlin.To assess the quality of the fit, we use the coefficient of determination R Vol:.( 1234567890 www.nature.com/scientificreports/

Statistical analysis
A two-way ANOVA analysis is used to calibrate the algorithm and to analyze the precision of the estimated parameters of the in-plane and out-of-plane fiber orientation distribution.Instead, a one-way ANOVA is used to evaluate the robustness of the algorithm.For all analyses, a post-hoc Tukey HSD (honestly significant difference) test is used to assess pairwise differences between groups.The results are considered significant at the 5% level.Due to the high number of pairwise comparisons, we use the compact letter display with uppercase letters to report the results.The data analysis is carried out employing the Real Statistics Resource Pack software 62 .

Artificial fiber dispersion volume
Using a custom Matlab code, we generated artificial grayscale 3D images with a number N of straight fibers with known fiber orientation distribution within a volume of M × M × M voxels.The midpoint (x, y, z) k of each fiber, k = 1, . . ., N , is chosen randomly according to a three-dimensional uniform distribution, while the orientation angles (θ, φ) k are sampled from the PDF given in Eq. ( 2) using a rejection method.The sam- pling is performed by generating uniformly distributed random points within the three-dimensional domain D = {(θ, φ, r)| − 90 • ≤ θ ≤ 90 • , −90 • ≤ φ ≤ 90 • , 0 ≤ r ≤ max (ρ)} and selecting only those points that fall under the function ρ(θ, φ).
Once the fiber positions (x, y, z) k and orientations (θ, φ) k are computed, a 3D binary image of M × M × M voxels is generated with an intensity Y = 1 for voxels belonging to a fiber and Y = 0 elsewhere.All fibers are of the same aspect ratio of more than 10, in order to enable a meaningful evaluation of the FOD using the discrete Fourier transform methods 53,54 .The image is then smoothed using a Gaussian kernel with a deviation of 0.7 voxel to reduce the sharpness of fibers that might introduce artifacts in the spectrum.This operation converts the image from binary to the grayscale image.Since smoothing reduces the maximum intensity Y max , all voxels are scaled by a factor Y −1 max to restore the maximum intensity to 1.

Collagen fiber tomography
Images of collagen fibers were acquired from human skin samples harvested from the abdominal region during a routine surgical procedure.The study was approved by the Regional Committee for Medical and Health Research Ethics (Project ID: 474249).All examinations were performed according to the rules for the investigation of human subjects set out in the Declaration of Helsinki.All study participants provided written informed consent.
In addition, we analyzed ex vivo skin samples from the dorsal region of mice.The skin tissue was obtained from mice used for in vivo studies in another project.After sacrifice by cervical dislocation under isoflurane anesthesia performed manually by trained personnel, the tissue was donated to the present project for ex vivo studies.The animal laboratory (Department of Comparative Medicine at Norwegian University of Science and Technology) is approved for animal tests by the Norwegian Food and Animal Safety Authority in document VSID 3506.In this document the laboratory's systems for animal care and ethics are approved according to the relevant national and EU regulations.
Samples were cleaned from adipose tissue and stored at −28 • within 2h of harvest.Prior to SHG imaging, tissue was thawed at room temperature ( 22 • ) and prepared according to the SeeDBp protocol 63 , consisting of a fixation step in 4% paraformaldehyde for 12h and 6 optical clearing steps in fructose solution in 0.1×PBS of increasing concentration.Specifically, the samples were incubated for 4h each in a 20% , 40% , 60% w/v solution, then for 12h each in a 80% , 100% w/v solution and finally for 24h in a 80.2% w/w solution (SeeDB solution).All steps took place at 25 • C.This technique allows biological tissues to be cleared without significant morphologi- cal changes, even in the case of fibrous tissue 64 .After completing the clearing process, the samples were placed in a press-to-seal silicone isolator (CoverWell™ Imaging Chambers, Grace BIO-Labs, Oregon, USA) filled with SeeDBp solution and sealed with two rectangular glass coverslips at the top and bottom.
For SHG imaging of collagen fibers, we used the confocal multiphoton microscope Leica TCS SP8 (Leica Microsystem, Germany) with a Leica HCX IRAPO 25× , NA 0.95 water objective with a working distance of 2.4 mm.The second harmonic of collagen is induced using a multiphoton laser source tuned at 890 nm (Chameleon Ultra I; Coherent Corp., Saxonburg, PA, United States) and the signal emitted at 445 nm is detected in the forward and backward directions.The images were acquired on a square target of 465 µ m × 465 µ m with a x-y resolution of 0.454 µm/px every 0.57 µ m (minimum vertical step of the microscope) in the z direction, for a total of ∼ 116.25 µ m scanned thickness.To compensate for blurring effects due to light diffraction from the confocal microscope, especially along the z direction, the 3D image was pre-processed in Fiji 65 using the open source DeconvolutionLab2 plugin for deconvolution microscopy 61 .In particular, we used the Richardson-Lucy algorithm with total-variation regularization (regularization parameter = 10 −3 , 30 iterations), adopting the Gibson and Lanny 3D optical model for the point spread function 66,67 .Then, to match the z spacing with the x-y resolution, the 3D image was resampled vertically using the imresize3 Matlab function with a stretch factor of �z/�x = 0.57/0.454≈ 1.26.

Data availibility
The data that supports the results within this paper are available from the corresponding author upon reasonable request.

5 Figure 5 .
Figure 5. Algorithm precision for in-plane parameters.Errors between the estimated (subscript e ) and true (subscript t ) in-plane parameters (mean and standard deviation of 10 images): (a) error in the in-plane concentration a = a e − a t ; (b) error in the in-plane angle �α = |α e − α t | (isotropic case a t = 0 omitted).Means not sharing uppercase letters differ significantly by the Tukey-test at the 5% significance level.Letters must be compared within the same concentration a t .

5 Figure 6 .
Figure 6.Algorithm precision for out-of-plane parameters.Errors between the estimated (subscript e ) and true (subscript t ) out-of-plane parameters (mean and standard deviation of 10 images): (a) error in the out-of-plane concentration b = b e − b t ; (b) error in the out-of-plane angle �β = |β e − β t | (isotropic case b t = 0 omitted).Means not sharing uppercase letters differ significantly by the Tukey-test at the 5% significance level.Letters must be compared within the same concentration b t .

Figure 9 .
Figure 9. Schematic representation of the orientation of the unit fiber N (in black), the global basis {E a } a=1,2,3 (in blue) and the principal basis {M , M ip , M op } (in red): (a) unit vector N in the principal frame; (b) unit vector N in the global basis; (c) rotation of the principal basis with respect to the global basis using the Tait-Bryan angles.Dashed vectors represent M ip and M op before the rotation about M.

Figure 10 .
Figure 10.Illustration of the 3D discrete Fourier transform algorithm: (a) signal g(x, y, z) of a 256 × 256 × 256 representative image in the discrete spatial domain with N = 14 fibers aligned to (θ, φ) = (60 • , 60 • ) ; (b) center- shifted spectrum |ĝ(u, v, w)| in the frequency domain, with a sketch of the overlapping funnel filter sampling the frequencies associated with the fibers with orientations within the cone, shown in (a), with mean direction (θ, φ) = (60 • , 60 • ) and opening ε ; (c) cross section of the filter in the local v-w plane at u = 128 .The filter is axisymmetric about the direction (θ, φ) (note that the symmetry axis shown is slightly skewed with respect to the v-w plane); (d) raw dFOD (discrete fiber orientation distribution) of the representative fiber dispersion shown in (a) using q = 2.4 , and D = 31 , corresponding to an angular resolution of �θ = �φ = 5.8 • ( ε is taken to be equal to �θ = �φ).