Stochastic modelling of 3D fiber structures imaged with X-ray microtomography

Many products incorporate into their design fibrous material with particular levels of permeability as a way to control the retention and flow of liquid. The production and experimental testing of these materials can be expensive and time consuming, particularly if it needs to be optimised to a desired level of absorbency. We consider a parametric virtual fiber model as a replacement for the real material to facilitate studying the relationship between structure and properties in a cheaper and more convenient manner. 3D image data sets of a sample fibrous material are obtained using X-ray microtomography and the individual fibers isolated. The segmented fibers are used to estimate the parameters of a 3D stochastic model for generating softcore virtual fiber structures. We use several spatial measures to show the consistency between the real and virtual structures, and demonstrate with lattice Boltzmann simulations that our virtual structure has good agreement with respect to the permeability of the physical material.


Introduction
Many products incorporate fibrous materials with carefully chosen permeability and absorption levels as a way to control liquid transport and retention in the product. The key to understanding the flow of liquid through the material is knowledge of the internal microstructure and the geometric properties of the fibers. Analysis of the performance of such materials and the transport of fluid through the structure has been considered by a number of authors [1][2][3].
Performing such analyses of the fiber structure requires the material to be well imaged and visualised. X-ray microtomography is a common method for visualizing such fiber materials, providing high resolution 3D image data sets that capture the geometry of the fibers within the structure [4]. This image data can be segmented to obtain a binary structure, which can be used to simulate fluid flow or other similar properties. A drawback of this approach, however, is that running simulations on different materials requires each to be imaged separately, which can be expensive and time consuming. This becomes more so if we wish to design a new material, requiring repeated production and imaging of different samples.
One possible solution to this challenge is to create a stochastic model that generates 3D virtual fiber structures intended to mimic the properties of the material in question. The advantage of a parametric model is that we can control the characteristics of the fibers, such as fiber orientation and curvature, by changing the values of the model parameters. Such an approach affords significant flexibility, enabling us to describe a variety of different fiber structures. Stochastic modelling of material structures in general is an extensively researched field of study [5,6], and modelling of fiber structures includes methods for non-woven fiber materials [7,8], composite fiber materials [9,10], and hardcore non-overlapping fiber structures [11].
To illustrate the methodology, a sample of fibrous material was imaged using X-ray microtomography and analysed to visualise its internal structure. Within the material, fibers are bound together to a certain extent to stabilise the structure. We therefore consider a softcore fiber model where the interpenetration of fibers is intended to mimic the fiber-to-fiber binding points. The use of stochastically generated random walks to create the fiber skeleton has been considered by a number of authors [11,12], and we use this approach in our method to control the bending of the fiber.
To fit our stochastic model to the data, we develop a method to estimate the parameters of the model from a sample of isolated fibers, and test the accuracy of the method against virtual fibers with known parameters. Identification of fibers in 2D images is well developed, but analysis of 3D structures requires greater care given the more complex geometry of the structure [13][14][15][16].
As our virtual structure is intended to mimic the physical material, we use three particular spatial measures to compare the two structures: fiber-to-fiber contact distance, the fiber K function, and intrusion porosimetry. We also run mass transport simulations on each structure using a lattice Boltzmann method [17], to demonstrate that our stochastically generated structure can be used to estimate the permeability of the real fiber material.

Image acquisition -X-ray microtomography
A sample of commercial fiber material provided by Essity (Gothenburg, Sweden) was scanned using X-ray microtomography, a nondestructive 3D imaging technique capable of revealing the internal structure of materials. The sample in question was a highly porous nonwoven structure used for absorbent products. A cut 10 mm diameter disk of fiber material was scanned in an up-right position and in local geometry (i.e., only the central portion of the sample was imaged and reconstructed) using a ZEISS Xradia Versa XRM520 system (Carl Zeiss, Germany) at the 4D Imaging Lab at Lund University. Due to the sample being scanned in local geometry, some fibers that are primarily outside the FoV may still appear in the images as small, divided fibers. As the length of these divided fibers is small, they do not contribute to the estimation method we will describe below.
The following scanning parameters were used; source-voltage: 40 kV, source power: 3 W, exposure time per projection: 1.25 s, total number of projections: 1001, total sample scanning time, including reference and alignment images: 1 h, 4x optic with a source-to-sample distance: 18.48 mm and sample-to-detector distance: 31.41 mm, Fieldof-View (FoV): 2.525 × 2.525 mm 2 (camera binning 4), and effective isotropic pixel size: 5 μm. The linear X-ray attenuation coefficient (LAC) [cm − 1 ] was reconstructed on a 32-bit grey level scale, with cubic voxels of dimension 5 μm, using the Zeiss reconstruction software.

Image processing and visualisation
The reconstructed data set was processed using standard image processing methods available in ImageJ/Fiji [18], including the application of a 3D median filter of pixel width 9 × 9 × 9 to remove noise and to facilitate the consecutive image segmentation using grey-level thresholding. The segmented image mask was processed using the Ori-entationJ function [19,20], which is based on evaluation of the gradient structure tensor and creates a colour-coded image based on the orientation of the fibers in 3D space. 3D renderings of the images are shown in Fig. 1.
Observing the sample at different perspectives, we see that the fibers travel from edge to edge; we can therefore assume the fibers are infinitely long, and thus do not need to consider a length distribution for the fibers. Individual fibers tend to remain oriented in one particular planar direction, and the fibers also demonstrate significant degrees of curvature and bending, in some cases curving back on themselves or becoming intertwined with other fibers.

Stochastic model
We construct our 3D fiber structure by first generating a set of nodes using a random walk, through which we plot a smooth curve that comprises our fiber skeleton. This skeleton is then dilated at a specified radius to give cylindrically shaped fibers. This process of generating fibers is repeated until the solid volume fraction matches a target value.

Random walk
We generate our virtual structure on a three-dimensional domain of size L X × L Y × L Z . As our virtual fibers will initially be generated in continuous space, the specific choice of values is primarily only important for the scale of the dimension dependent parameters we will introduce below.
Each random walk is characterised by an initial startpoint s 0 which is uniformly distributed in the domain. To build the fiber skeleton, we evolve two random walks in different directions from this startpoint until they hit the domain boundary. We denote the m-th point of the first generated random walk as s A m and the n-th point of the second generated random walk as s B n ; thus In defining the next (xy)-coordinate in our random walk, we use a von Mises distribution to sample the angle for each step [21,22], with probability density function where I 0 (κ) is the modified Bessel function of order 0. The parameters of the von Mises distribution are the mean direction μ and the concentration κ, with the latter controlling the global bending tendency of the fiber; taking a larger value for κ will result in straighter fibers.
To generate the first point s A 1 , since we do not have an existing direction to use, we assume the angle θ is uniformly distributed on the circle and thus take μ = 0 and κ = 0. For the remaining points in the first random walk, we use the specified concentration value κ and define the mean direction using the angle between the previous node points: .
We can then sample our bending angle from the density function (2) and calculate the step distances for the x-and y-coordinates of the next node point as with given step distance r. A diagram of these steps is shown in Fig. 2.
For the z-coordinate, we recall that we saw in our visualisation of the fiber sample in Fig. 1 that the fibers tend to be oriented parallel to the (xy)-plane, with 77% of the fibers having a total elevation change of less than 20 • , and thus the change in the z-coordinate of the node points along the curve will be minimal. We therefore sample the change in zcoordinate as z step ∈ N (0, σ 2 ), assuming a mean of zero and given variance σ 2 .
The generation of new node points for the first random walk is continued until a node is generated outside of the specified domain limits. For the second random walk, we calculate the mean direction for the first node point s B 1 using the distance between the points s A 1 and s 0 from the first random walk: The fraction is inverted in this calculation compared to (3) since we are travelling from the point s A 1 to s 0 . We can then once again calculate the (xy)-step distances for the new node point. These steps are shown in Fig. 3. We continue to generate subsequent points as in the steps detailed above, thereby giving us our two branching random walks, which we stitch together to create a single random walk (see Fig. 4).

Bézier curves
We generate our fiber skeleton using a Bézier curve, a method for generating smooth parametric curves through a set of control points P 0 , …, P n , where n is the order of the curve [23,24]. The resulting curve is given by the equation   In our case, we will use quadratic (n = 2) Bézier curves for the first and last node point pairs and cubic (n = 3) Bézier curves for the inner pairs. Here, B(t) will be a vector of (xyz)-coordinates. The points P 0 and P n will be each pair of neighbouring node points from our random walk, while the remaining control points will be defined between these points and provide directional information for the curve; in the following, we will refer to these as support points.
We specify the support points at each node point using the distance between the previous and subsequent node points; support points will therefore only be generated for the inner points of the random walk. This distance is weighted by a curvature parameter c; the smaller the parameter is, the less curvature our Bézier curve will have. This parameter therefore enables us to capture and control the local bending tendency of the fibers, that is the bending of the fibers between the node points.
Considering the first three points (s 1 ,s 2 ,s 3 ), we can calculate the two support points for the inner node point s 2 as follows: The Bézier curve for the first pair of node points (s 1 , s 2 ) is thus given by Taking an example curvature value of c = 0.5, we can determine the support points for s 2 and the resultant Bézier curve, as shown in Fig. 5.
For the inner point s 3 , we can calculate the support points (s 3,1 ,s 3,2 ) in the same way as above with the inclusion of the point s 4 . We then have a set of four points, and hence we use a cubic Bézier curve: Applying this method to the remaining node points, we can build the full fiber skeleton, as shown in Fig. 6. We note that the Bézier curve between the final node points (s 8 , s 9 ) will be quadratic, as for the first segment.

Cylindrical fibers
In generating our cylindrical fibers, we first consider a discrete lattice grid with resolution nx × ny × nz and a corresponding binary matrix M that represents the generated structure. We discretise the points from our Bézier curves by first mapping the points from the original domain to the new ranges and then mapping these newly defined points to the closest grid point on the lattice. Elements of the matrix M with value 1 then indicate the presence of a fiber at the corresponding coordinate on the lattice grid.
The procedure thus far yields a fiber skeleton structure; the final step is to dilate the skeleton at a specified fiber radius, denoted by f r . Our method assumes that every element in the matrix is dilated at the same radius, and for simplicity we also assume that every fiber has the same radius. The analysis in the next section proves that this assumption is reasonable for the fibrous material sample we are considering.

Isolating fibers
To isolate the individual fibers in the material sample shown in Fig. 1, we first apply a skeletonisation to the segmented structure. Then, we define an indicator for each point on the skeleton corresponding to the number of fiber points in its 3 3 neighbourhood. We define a crossing or intersection point as any point on the skeleton for which the value of the indicator > 3, and all such identified points are then removed from the skeleton. From the resulting structure, the individual, nonintersecting fiber segments can then be identified using 26-connected pixel connectivity. Finally, to ensure we do not retain any artefacts of the skeletonisation, any segment shorter than a specified length is removed.
To determine the appropriate reconnection of the isolated segments and rebuild the fiber structure, we follow the approach detailed in [13,14], in which a weighting is calculated for each proposed reconnection based on the distance and angle between the two endpoints. We give a brief overview of this approach here, with a slight change made to the treatment of the relative angle between the endpoints.
For a given set of isolated segments {p 1 ,…,p n }, we want to determine all potential reconnections for the endpoint l i of a segment p i . Since the point cannot be connected to itself, we consider the set of possible endpoints {l 1 , …, l n }⧹{l i } and determine those within a sector of length l max and angle θ max . The proposed connection between the endpoints l i and l j is then characterised by the absolute length l ij of the connecting segment and the angle θ ij between the two endpoints. For the method detailed in [13], this angle is calculated as the mean of the absolute value of the angles θ(l i , l ij ) and θ(l ij , l j ) between the corresponding line segments; for our approach, we adjust this method and instead consider θ ij as the absolute difference between the angles θ i and θ j at each endpoint. Since the maximum difference between the two angles is π, our weight for the proposed connection between endpoints l i and l j therefore becomes Fig. 5. We specify the support points for s 2 with curvature c = 0.5 and plot the resulting Bézier curve through the first pair of node points. Fig. 6. We stitch together the individual Bézier curves through each of the node points of our random walk to build our full fiber skeleton. For our example fiber skeleton here, we have used a constant curvature value of c = 0.5.
The change to the second component of ω ij places greater weight on connections between closely parallel fiber segments, regardless of the angle of the connecting segment.

Parameter estimation 2.5.1. Estimating random walk nodes
We identify the initial node point estimate of an isolated fiber segment using a measure of local curvature at each discretisation point. We fit an osculating circle to each point of the fiber following the approach detailed in [25] with order m, and define the local point curvature as the inverse of the radius. In general, m will depend on the step distance of the fibers (see below), but for our purposes we will consider a fixed order m = 20. We then take the element with the highest point curvature as the first node point estimate (see Fig. 7). This approach ensures we do not miss parts of the fiber with high bending tendency.
To determine the remaining nodes, we identify points backwards and forwards along the fiber that are a given step distance away from the previous node estimate, starting from our initial node. The step distance is chosen based on the scale of the domain and the average length of the fibers. In testing the accuracy of this approach in Section 3 against virtual fibers with the exact step distance used to generate the fibers, it was found that node points would often be missed due to the discretisation of the fiber; therefore, in implementing the method, we label an element as a node point if it is within a small margin of the given step distance.

Estimating fiber curvature
We estimate the curvature of the fiber segment by fitting a series of Bézier curves through the identified node points for a given range of curvature values. We calculate the overlap of this curve with the initial fiber by comparing the number of points intersected by both, and then pick the curvature value that gives the greatest overlap. To ensure a better fit of our curves, we include two additional node points at the start and end of the fiber. Through testing against example virtual fibers, it was found that placing these additional nodes at the first and last fiber elements was sufficient.

Estimation method accuracy
The accuracy of each component of our parameter estimation method was tested by generating virtual fiber samples using our sto-chastic model. For each test, the value of one individual parameter was varied at a time within a given range while the remaining two parameters were fixed at a default value. We also used a fixed step distance of r = 100. Since we know the exact parameter values used to generate the virtual fibers, we can directly compare the accuracy of our estimation, and changing just one parameter at a time ensures any effects on the accuracy are solely the result of the parameter we want to measure. The parameter ranges and default values are shown in Table 1.
For estimating the mean and variance of the z-step, we fit a normal distribution to the values using the inbuilt Matlab function fitdist [26]. For estimating the von Mises concentration parameter κ, we follow the maximum likelihood approach detailed in [27]. Finally, for estimating the curvature c, we calculate the mean, median, and mode for each test value, and compare how each statistical measure performs against the exact value. For testing each parameter combination, we generate 1000 virtual fibers and exclude those that are not of a suitable length to estimate the parameters due to an insufficient number of node points; in our testing, around 20-25% were found to be unsuitable. The comparisons to the exact values are shown in Figs. 8-10.
The results for the estimation of the mean, variance, and concentration all compare well with the exact values, though in each case the accuracy of the estimation method decreases as the parameter values increase. This drop off is, however, within reasonable limits given the ranges we are considering for our test parameters. For the estimation of the curvature value, each of the three statistical measures we use accurately perform well up to around c = 0.5. However, above this point Fig. 7. We plot the local point curvature for an example virtual fiber in the region of the peak value, highlighted in red. Here, vector element refers to a unique point in the ordered list of grid points after the fiber has been discretized.

Table 1
Ranges for the parameters used in testing the method. the measures become significantly less accurate and also inconsistent with one another; the two outliers in the curvature estimation using the modal value are a result of the individual estimates having peaks at c = 0, which in some cases are the largest peak and thus incorrectly chosen as the estimate. Therefore, in running our estimation method against real structures, we will consider all three measures together with the full histogram of estimates to determine which value to use.

Analysis of material sample
The method detailed in Section 2.4 for isolating individual fibers was applied to the segmented image data sets of the fibrous material sample. Initial analysis identified 1655 individual segments, which, after reconnection, resulted in 420 individual fibers with an average length of 190 elements and the largest of length 1558. A comparison of the original segmented image and the individual isolated fiber skeletons is shown in Fig. 11.
An initial analysis of the structure was run using the software MIST, a program for the visualisation and characterisation of 3D geometries [28]. A pore-size distribution analysis was run on the inverted structure to determine an appropriate fiber radius of 2 voxels, consistent with the average fiber diameter, and the solid volume fraction of the fiber structure (after removing areas of open space at the top and bottom of the structure) was found to be 2.5%. The dimensions of the reduced structure, which were used in generating our virtual fiber structure, are 505 × 497 × 331 voxels.
For the estimation of the remaining model parameters, we ran the isolated fibers against the method detailed in Section 2.5. Given that the average size of the isolated fibers was 190 elements, we used a specified step distance r = 50, as a reasonable compromise between maximising the number of usable fibers and matching their relative scale. The outputs from applying the parameter estimation method to the real structure are given in Table 2.
In total, of the 420 isolated fibers, 218 were found to be of sufficient length to be used for the estimation. This is less than in our simulation study, when around 75-80% were suitable. However, given that the isolated fibers are less well-defined, the number of fibers is large enough for our purposes. The mean value and variance for the z-step both being close to zero is consistent with our assumption that the fibers in the material stay relatively flat in the z-axis, and the value of the concentration κ is in the range expected for fibers with a high degree of bending tendency.
For the fiber curvature, while the mean and median values were close to one another, the modal value was significantly different. These values indicate a peak at c = 0, with the remaining estimates more spread out among our test values. This suggests that using a single curvature value for all virtual fibers is inconsistent with the actual material; a potential extension would therefore be to sample the curvature value from a given distribution. In generating our virtual 3D fiber structure, we will use the mean value, and thus take c = 0.3956. A visualisation of the structure is shown in Fig. 12; the structure was postprocessed after being generated to apply some smoothing to the fibers, though the solid volume fraction of the original structure was maintained. The total computation time, including both generation and postprocessing, was around 30 min locally on a personal computer.   10. We consider three statistical measures for estimating the curvature parameter c. We see that each does well for smaller values, particularly taking the modal value, but there is a large drop off in accuracy for larger values.

Comparison of structures
We compare our physical and virtual structures using three spatial measures: fiber-to-fiber contact distance, the fiber K function, and intrusion porosimetry. These particular measures were chosen to enable a comparison of a range of aspects of the structures. Since our main focus is on fluid flow through the structures, we also run mass transport simulations on each structure and compare the permeability.

Fiber-to-fiber contact distance
The fiber-to-fiber contact distance measures the distribution of distances between fiber intersection points. The contact distance is a good spatial measure for comparison to the fibrous material sample, which has a high degree of intersection and layering of fibers.
We first skeletonise the fiber structure and isolate the crossing points in the skeleton as we did in Section 2.4. Since there can be a cluster of points identified as crossings at each intersection, we calculate the centroid of each cluster. We then generate a network over the fiber skeleton, using the centroids as the nodes and where the edges of the network connect each node to its closest neighbour. After removing duplicate connections, we calculate the length of each edge and generate a distribution of contact distances.
Comparison of the distribution of contact distance for the real and virtual fiber structures is shown in Fig. 13. We see that the measure compares well between the two structures, though there is a slight shift in the virtual structure toward longer contact distances. Fig. 11. Comparison between the orientation colour-coded fibers of the physical material and the corresponding skeletons of the individual fibers after they have been isolated and reconnected using the method presented above. The dimension of the extracted volume-of-interest is 2.5 × 2.5 × 1.4 mm 3 .

Table 2
Output values from the estimation method.   13. Comparison of fiber-to-fiber contact distances for the real (blue) and virtual (red) structures shows good consistency between the two structures, with the virtual structure having a slight shift toward longer contact distances.

Fiber K function
Ripley's K function [29] is a method from spatial statistics that can be used to determine the level of dispersion or clustering of a given point process. An analogous K function for fiber processes was presented in [30] and which we use for comparing our real and virtual fiber structures. The fiber K function we use is given by where Φ f is our fiber structure, A the domain of interest with volume |A|, Φ f (⋅) the number of fiber voxels in the given domain, and S(x, h) a sphere centred at x with radius h. We remove the need for edge-correction in the function by considering A to be a subdomain of our full structure, and we take a maximum distance of h max = 80 based on the largest fiber-to-fiber contact distance found in our previous analysis. Comparison of the fiber K function for each structure is shown in Fig. 14. The two structures show almost identical profiles, indicating that our virtual structure is a good representation of the real material. The profile itself is also consistent with our expectation, showing greater clustering of fiber points at greater distances.

Intrusion porosimetry
The intrusion porosimetry measures at every point p in the pore space the diameter of the largest sphere that can travel to p starting from a given plane. The measure provides an indication of the porous nature of the material, which is particularly important for the structures we are looking at given our interest in understanding fluid flow through the material. The porosimetry through both structures was calculated using MIST, and we consider a subdomain of the real structure to ensure all areas of open space that can affect the analysis are removed. The domain for the virtual structure was similarly reduced to match the size of the real structure, and doing so enables the intrusion porosimetry to be calculated through two non-intersecting areas to see whether there is any variation across the structure.
The results of the analysis are shown in Table 3 and visualised in Fig. 15, with points highlighted in darker blue indicating that a larger sphere was able to access the corresponding point p in the pore space; fibers are indicated in bright red. The measure shows good consistency between the real structure and the two samples from the virtual structure, with the average and standard deviation of the pore-size both being close in value in each instance. Visual inspection of the results also shows good comparison, with areas further away from the starting plane at the top of the structure having a lower value (towards the redder end of the scale) as expected. Consistency of the analysis through the two samples from the virtual structure also demonstrates the homogeneous   Fig. 15. Comparison of the intrusion porosimetry (in voxels) for the real structure (a) and the virtual structure samples (b) and (c); areas of darker blue indicate a larger sphere is able to access the corresponding point p in the pore space, with redder areas indicating a smaller accessible sphere and the fibers indicated in bright red. The line plot in (d) shows a comparison of the accessible volume fraction for a given probe diameter for each structure. The average pore-size, standard deviation, and maximum value for each structure is included in Table 3.
nature of the structure.

Mass transport simulations
Our final step is to run simulations of fluid flow through each structure to compute and compare the fluid permeability. These simulations were run using Gesualdo, an in-house software using the lattice Boltzmann method to compute mass transport properties of porous materials [1,17]. We consider a subdomain of the full structure as we did in calculating the intrusion porosimetry, and in each case the simulations were run until convergence was reached.
Comparison of the simulation results for each structure are shown in Table 4, and we can see that the estimates for the permeability match very well across all three simulations. The porosities for the virtual samples are consistent with expectations given our target solid volume fraction of 2.5%, with the slight differences being a result of considering a subdomain of the full structure.

Conclusion
A stochastic model for generating 3D fiber structures was constructed and fitted to an X-ray microtomography image data set of a sample of fibrous material used for absorbent products. A method for isolating fibers was applied to the physical structure, and the parameters for the model were estimated from these individual fiber segments. The estimation method was tested on virtual fibers with known parameters and shown to approximate the exact values to a reasonable degree of accuracy. The real and virtual fiber structures were compared using several measures, and the virtual structure was found to be a good approximation of the physical structure across every chosen measure. Fluid flow simulations were also run on each structure, and the permeability of the physical material was found to be approximated well by the virtual structure.
The benefit of being able to generate 3D fiber structures with the stochastic model is not just to be able to approximate a real material, but also to generate new and interesting structures. An interesting topic for future research would be to generate a range of fiber structures with different parameter combinations and analyse the effect of changing these parameters on the properties of the structure. This analysis can then be used to inform the future design of new materials. Additionally, while the method as detailed was developed for the specific structure considered, it could also be easily extended to consider a wider range of similar fiber materials. Some simple adjustments to the random walk can be used to analyse periodic structures or those with fibers that do not span the boundary of the structure as with the sample materials, and a more localised modelling of the fiber radius could be used to analyse materials with defects in the fibers.

Data availability
The raw/processed data required to reproduce these findings cannot be shared at this time due to legal or ethical reasons.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.