Dorian Depriester

Sabine Rolland Du Roscoat
email: sabine.rollandduroscoat@3sr-grenoble.fr

Laurent Orgéas

Christian Geindreau

Benjamin Levrard

Florian Brémond

Modelling the physical behaviour of fibrous materials still remains a great challenge because it requires to evaluate the inner structure of the different phases at the phase scale (fibre or matrix) and the at constituent scale (fibre). X-ray Computed Tomography (CT) imaging can help to characterize and to model these structures, since it allows separating the phases, based on the grey level of CT scans. However, once the fibrous phase has been isolated, automatically

Individual fibre separation in 3D fibrous materials imaged by X-ray tomography separating the fibres from each other is still very challenging. This work aims at proposing a method which allows separating the fibres and localizing the fibre-fibre contacts for various fibres geometries, that is: straight or woven fibres, with circular or non circular cross sections, in a way that is independent of the fibres orientations. This method uses the local orientation of the structure formed by the fibrous phase and then introduces the misorientation angle. The threshold of this angle is the only parameter required to separate the fibres. This paper investigates the efficiency of the proposed algorithm in various conditions, for instance by changing the image resolution or the fibre tortuosity on synthetic images. Finally, the proposed algorithm is applied to real images or samples made up of synthetic solid fibres. Ω orph Subset of voxels in Ω f which are not allocated to any fibre.

P Considered voxel in Ω f . ψ Azimuth angle.

r Mean fibre radius.

θ Elevation angle (latitude).

t i i-th computation direction for dc P (t i ).

ū Mean orientation of a single fibre.

u I P Minor axis of the inertia matrix I P .

subscripts: Voxel index (P ) or fibre index (k)

superscripts: Computational direction index (i) or eigenvalue indices (I, II and III)

Introduction

It is well known that the macroscopic behaviour of fibrous materials is inherited from their microscopic composition and microstructure [START_REF] Catalin | Mechanics of Random Fiber Networks: Structure-Properties Relation[END_REF] . Their physical properties depend on the volume fraction of fibres but also on their geometries and 3D spatial repartitions (orientation and fibre-fibre contacts). Their microstructures have been studied for the last 20 years using X-ray micro-tomography 2,3,4 . This technique allows imaging the inner structure of materials that are opaque to visible light and gives valuable information about the 3D spatial distribution of each phase. Nevertheless, the analysis of fibre-fibre contacts and the geometrical characterization of the fibres is still challenging because it requires to numerically separate the fibres from Computed Tomography (CT) images.

For solid circular fibres, recursive erosions can be applied for separation; then, the connected regions can be labelled prior to dilation [START_REF] Masse | Influence of relative density on the architecture and mechanical behaviour of a steel metallic wool[END_REF][START_REF] Rosini | In situ statistical measurement of local mor-phology in carbon-epoxy composites using synchrotron x-ray computed tomography[END_REF][START_REF] Arthur Hessman | Microstructural analysis of short glass fiber reinforced thermoplastics based on x-ray micro-computed tomography[END_REF] . Another method consists in first extracting the centerlines of the fibres, also called skeletonization algorithm [START_REF] John | The image processing handbook[END_REF] . Usually, the skeleton is defined by a series of control points and a piecewise-defined and smooth function (e.g. splines). Once the skeleton is computed, the main branches can be used to define each fibre. The skeletonization can be made manually [START_REF] Elg Wernersson | Characterisations of fibre networks in paper using micro computed tomography images[END_REF] or automatically [START_REF] Requena | 3d-quantification of the distribution of continuous fibres in unidirectionally reinforced composites[END_REF][START_REF] Latil | Towards the 3d in situ characterisation of deformation micro-mechanisms within a compressed bundle of fibres[END_REF] . In the general case, one would use the 3D euclidean distance map [START_REF] Bonnassie | A new method for analyzing local shape in three-dimensional images based on medial axis transformation[END_REF][START_REF] Orgéas | In-plane conduction of polymer composite plates reinforced with architectured networks of copper fibres[END_REF] or homotopic thinning [START_REF] Lux | 3d characterization of wood based fibrous materials: An application[END_REF][START_REF] Tan | Analysis of tomography images of bonded fibre networks to measure distributions of fibre segment length and fibre orientation[END_REF][START_REF] Yang | Three-dimensional image analysis of fibrous materials[END_REF] to skeletonize the fibrous phase. However, skeleton-based methods are restricted to low volume fraction and fibres exhibiting nearly isotropic cross sections.

For hollowed fibres, such as wood pulp fibres, specific algorithm may be used, usually based on the lumen tracking [START_REF] Svensson | Using distance transform based algorithms for extracting measures of the fiber network in volume images of paper[END_REF][START_REF] Malmberg | Measurement of fibre-fibre contact in three-dimensional images of fibrous materials obtained from x-ray synchrotron microtomography[END_REF] .

In more complex fibrous media, other methods are based on the local orientation of the structure. Some authors [START_REF] Krause | Determination of the fibre orientation in composites using the structure tensor and local x-ray transform[END_REF][START_REF] Naouar | Meso-scale {FE} analyses of textile composite reinforcement deformation based on x-ray computed tomography[END_REF] use the grey-levelled images for evaluating a structure tensor, computed from the convolution of gaussian distributions in the vicinity of each considered voxel. Eberhardt and Clarke 21 used the chord lengths whereas Naouar et al. [START_REF] Naouar | 3d composite reinforcement meso f.e. analyses based on x-ray computed tomography[END_REF] computed the Grey Level Co-occurrence matrix (GLC Matrix) [START_REF] Robert | Image segmentation techniques[END_REF] . Miletić et al. [START_REF] Miletić | Gradient-based fibre detection method on 3d micro-ct tomographic image for defining fibre orientation bias in ultra-high-performance concrete[END_REF] used the greyscale gradient of the images to track the surface of each fibre.

In binarized images, orientation characterization is usually based on the morphologies of the objects. Sandau and Ohser 25 proposed to evaluate the local orientation from the oriented chord-lengths.

Some authors used the last method to compute either the local inertia matrices [START_REF] Altendorf | 3d directional mathematical morphology for analysis of fiber orientations[END_REF] or the structure tensor [START_REF] Viguié | Finding fibres and their contacts within 3d images of disordered fibrous media[END_REF] . Miettinen et al. [START_REF] Miettinen | Non-destructive automatic determination of aspect ratio and cross-sectional properties of fibres[END_REF] also used the gaussian convolution to compute the structure tensor from binarized images.

In each method resulting in a structure tensor or an inertia matrix, the local orientation can be evaluated by solving the eigenvalue problem. Once the local orientations are characterized, the separation of the fibres is usually made based on the local changes of these orientations. Naouar et al. [START_REF] Naouar | 3d composite reinforcement meso f.e. analyses based on x-ray computed tomography[END_REF] used the homogeneity parameter [START_REF] Haralick | Textural features for image classification[END_REF] to separate warp and weft yarns in CT grey-level images of a 3D textile. Viguié et al. [START_REF] Viguié | Finding fibres and their contacts within 3d images of disordered fibrous media[END_REF] evaluated the local change (called "mean angular deviation") along the horizontal plane of this axis on binarized images. Thus, it was possible to separate the fibre in CT images by thresholding the mean angular deviation. After separation, the fibres were fully reconstructed using a recursive anisotropic dilation (dilation along the mean fibre direction). This method lead to valuable results on various media but it was still restricted to networks with planar fibre orientations.

It is worth mentioning that recent advances in Artificial Intelligence (AI) have brought novel techniques for automatic detection of shapes or patterns in 2D or 3D images; still, only few authors [START_REF] Hong | Effective segmentation of short fibers in glass fiber reinforced concrete's x-ray images using deep learning technology[END_REF][START_REF] Ivanov | Extraction and quantification of features in xct datasets of fibre reinforced polymers using machine learning techniques[END_REF] used convolution neural network to seperate fibres in CT images. This is because of the large amount of data required for training the neural network.

In this work, the separation criterion proposed by Viguié et al. [START_REF] Viguié | Finding fibres and their contacts within 3d images of disordered fibrous media[END_REF] was revisited for fully 3D oriented medium. Then, the dilation procedure used by the last authors was also improved to make it valid for either straight or woven fibres with 3D local orientation. We first introduce the separation criterion, the reconstruction algorithm and some metrics applied to the fibres. In a second task, we investigate the influence of the parameters used for the separation and the reconstruction of the fibres on the resulting geometry. Then, the proposed algorithm is applied to different media in order to evaluate its ability to separate the fibres in different conditions.

Method 2.1 Example of fibrous medium

In this section, the fibrous medium illustrated in Figure 1 and a mean length of about 10 mm, as detailed elsewhere [START_REF] Orgéas | In-plane conduction of polymer composite plates reinforced with architectured networks of copper fibres[END_REF] . Consider that CT image as a 3D regular grid of points (voxels), which belong either to the matrix phase or the fibrous phase. Let Ω f be the subset of voxels corresponding to the fibrous phase.

The 3D images are stacks of 2D images (slices) where e 1 and e 2 denote the directions associated to each slice, whereas e 3 denotes the stacking direction.Since the proposed approach is independent of the voxel size, all dimensions shown in figures are given in voxels throughout this paper.

Description of the overall algorithm

Figure 2 represents the method to identify the fibres and the fibre-fibre contacts from a binarized CT image, that is a 3D image where the fibrous phase is already identified. This algorithm is composed of three steps: step 1 refers to the separation of fibres, step 2 refers to the reconstruction of the fibres and step 3 deals with the analyses of contacts and fibres geometries. where n t denotes the number of investigated directions t i .

Local orientation

In P ∈ Ω f , Altendorf and Jeulin 26 define the inertia matrix I P such that:

I P = 1 n t nt i=1 dc P t i 2 M t i (2)
with: where t i denotes the euclidean norm of t i .

M t i = 1 t i 2       t i 2 2 + t i 3 2 -t i 1 t i 2 -t i 1 t i 3 -t i 1 t i 2 t i 1 2 + t i 3 2 -t i 2 t i 3 -t i 1 t i 3 -t i 2 t i 3 t i 1 2 + t i 2 2       (3)
Let I I P ≤I II P ≤I III P be the eigenvalues of I P and u I P , u II P , u II P the corresponding normalized eigenvectors.

For a given fibre, its orientation is characterized by its lowest inertia momentum. In other words, the minor axis of I P (denoted u I P above) gives the local orientation of the fibrous phase in P .

Figure 2b illustrates the local orientations u I P (white arrows) computed in the medium introduced in Figure 1. It is clear that, away from the contact zones, local orientations appear to be fairly aligned with the macroscopic orientation of each fibre. On the contrary, contact zones (e.g., with coordinates around (35, 95) in Figure 2b) are characterized by large changes in the local orientations.

Misorientation angle

For a given voxel P ∈ Ω f , let N P be the subset of voxels in Ω f and in contact with P , considering its 26-connected neighbourhood. Let α be the misorientation angle, expressed in radian-per-voxel, so that: αP = 1 δ

Q∈N P u I P , u I Q PQ with: δ = P ∈N P 1 PQ (4)
where u I P , u I Q denotes the bearing angle between u I P and u I Q . This definition of the misorientation angle is a 3D generalization of the mean angular deviation introduced by Viguié et al. [START_REF] Viguié | Finding fibres and their contacts within 3d images of disordered fibrous media[END_REF] . Here, the angular values are weighted by the inverse of the distance from the central point.

Figure 2b illustrates the misorientation map (colormap), evidencing the previous observations (e.g. coordinates around (35, 95)). In the investigated media, this misorientation angle is typically larger than 5 • near contacts. On the opposite, the misorientation angle is close to 0 • away from the contact zones.

Thresholding and labelling

Based on the previous considerations, one can separate fibres in contact using a threshold value ᾱth on α, as illustrated in Figure 2c.

Once the fibres are separated, each one can be considered as a set of voxels, defined as a 3D connected region. Therefore, the fibres can be labelled using a standard flood-fill algorithm [START_REF] Soille | Morphological image analysis: principles and applications[END_REF] , associating a unique integer to each connected region, as illustrated in Figure 2d. In the present paper,

(a) Ω orph (i) D (F 1 ) D (F 2 ) F (i) 1 F (i) 2 (b) F (i+1) 1 F (i+1) 2 Ω f mult Ω orph (i+1)
Figure 3: Schematic representation of the boolean operations used for dilation: two fibres are represented (F 1 and F 2 ); the dilation is restricted to Ω orph , as defined in ( 5). (a) at the i-th dilation step and (b) results at the next one. The set Ω f mult represents the points that belong to both the dilated objects (introduced in section 2.4.2).

F k and n f denote the subset of voxels which belong to the fibre labelled k and the number of identified fibres, respectively.

When thresholding α, some very small disconnected regions may arise. In order to remove these unrealistic outliers from the labelled regions, regions smaller than a given number of voxels were discarded. The number of discarded regions and their associated volumes mainly depend on the 125 image resolution and the number of contacts.

Step 2: Fibre reconstruction

As illustrated in Figure 2d, at this stage of the algorithm, about 30% of the voxels belonging to the fibrous phase are not affected to a labelled fibre. Let Ω orph be this set, consisting in n orph orphan voxels. To allocate these voxels to a fibre, recursive dilations 27 on each identified fibre k are applied:

F k (i+1) = F k (i) ∪ D F k (i) ∩ Ω orph (i) ( 5
)
where i is the iteration step and D is the dilation operator, described hereafter. Figure 3 illustrates the previous equation as a Venn diagram. This operation is iterated until Ω orph is stabilized. Algorithm 1

gives the pseudo-code associated to this procedure. A dilation operation is defined for a given kernel.

Algorithm 1 Iterative fibre growth: this operation is performed until the number of elements in Ω orph (n orph ) converges.

Ω orph ← orphan voxels n orph ← card (Ω orph ) while n orph = constant do for all k = 1 to n f do dilation of F k (eq. 5) end for Ω orph ← orphan voxels n orph ← card (Ω orph ) end while

Isotropic and longitudinal dilations

The first kernel is an isotropic kernel constituted of 6-connected neighbourhood. The effect of using the isotropic kernel for dilation is schematically illustrated in Figure 4a. the ones of longitudinal dilations, mainly because of the size of structuring elements used for the dilations. On the contrary, at the end of the process, the longitudinal dilation leaves some orphan voxels whereas n orph converges toward 0 when using the isotropic dilations. The differences between the two kinds of dilation are discussed in section 3.2.

Figure 5 also illustrates the results from longitudinal dilation after separation on the example medium (introduced in Figure 1) in terms of geometry. According to this figure, it is clear that the 150 procedure quickly converges toward the final geometry as n orph converges toward its final value.

Multiple points

When recursively applying the dilation operation, some points may be allocated to several fibres at once, as schematically represented in Figure 3 and depicted in Figure 6. Let Ω f mult be this set of multiple points, as defined below:

Ω f mult = n f k=1 n f k =k+1 F k ∩ F k . ( 6
)
Figure 6: Multiple points resulting from the recursive dilations on copper fibres (1471 points).

The proposed method to allocate these multiple-points in a deterministic manner is the "leastmisorientations method". To ascribe these voxels to a single fibre, the principle is first to temporary 155 allocate the multiple-points to one fibre, then to recompute the local orientations and the misorientations (as defined in section 2.3.2) in the considered fibre. The same is done with all candidate fibres, allowing to evaluate its resulting misorientations for each fibre. Finally, this voxel is allocated to the fibre leading to its least misorientation. In other words, multiple points are allocated to the fibres in a way that maximizes the straightness of each fibre. This procedure is summed up in Algorithm 2.

Algorithm 2 Least-misorientation method. A P ,k denotes the tabular data, containing the misorientation value in P assuming that P ∈ F k .

for all k = 1 to n f do F tmp ← F k ∪ Ω f mult // Consider that Ω f mult belongs to F k Compute α, considering only F tmp for all P ∈ Ω f mult do A P ,k ← αP //
Store the resulting misorientation at P end for end for for all P ∈ Ω f mult do K ← arg min k (A P ,k ) // P ∈ F K minimizes the misorientation F K ← F K ∪ {P } // Add P to Fibre K end for 2. 5 Step 3: Geometrical analyses

Mean orientation of a fibre

The method proposed by Woodcock and Naylor 34 was used to evaluate the mean orientation of each fibre. Let A k be the mean fabric tensor of fibre F k , so that:

A k = 1 n k P ∈F k u I P ⊗ u I P ( 7
)
where n k is the number of voxels in F k and ⊗ denotes the tensorial product. The largest eigenvalue of A k is an estimate of the distribution variance (comprised between 1 /3 and 1) and its associated eigenvector, denoted ūk below, is an estimate of the mean orientation of F k 19 .

165

Mean radius

Consider a voxel P ∈ Ω f and its parent fibre F k with mean orientation ūk . For slightly woven fibres, Altendorf and Jeulin [START_REF] Altendorf | 3d directional mathematical morphology for analysis of fiber orientations[END_REF] give the following equation to evaluate the local radius of the fibre F k in P ,

given by the half-chord length measured along t i :

∀P ∈ F k , r P t i = dc P t i sin ϕ i k (8)
with ϕ i k = (ū k , t i ). The more parallel t i and ūk are, the less the evaluation of r P (t i ) is relevant.

Hence, the evaluation of the mean radius was made by weighting the local radii based on the sine of

ϕ i k : ∀P ∈ F k , rP = 1 i sin ϕ i k nt i=1 r P t i sin ϕ i k . (9)
Let rk be the mean radius of F k . It comes:

rk = 1 n k P ∈F k rP . (10)

Contact identification

Let F i and F j be the volumes associated to two distinct fibres i and j. The set of voxels which belong to the fibre i and are in contact with the fibre j can be defined as follows:

C i,j = F i ∩ D iso (F j ) . ( 11
)
3 Choice of the parameters of the algorithm

The objective of this section is to propose guidelines to use the algorithm summarised in section 2 and to gauge the errors made on the identification of the fibres in terms of geometry (length and diameter) and location (position and orientation). For that purpose, the following methodology is applied: (i) a synthetic medium, representative of one of the various fibre networks mentioned in the introduction, is generated, (ii) the influence of one parameter on the fibre identification is gauged and (iii) guidelines are provided.

Computational directions

Example medium

In order to investigate the effect of increasing the number of computational directions of the chord length (denoted n t in this paper) on the computation of the mean orientation, a "hedgehog-like" fibrous medium was numerically generated, as depicted in Figure 7a. In order to avoid the influence of contacts on the evaluation of local orientations, each fibre of this hedgehog was not in contact to each other. These fibres have a nominal length of 35 voxels and two nominal radii are investigated, respectively 3 and 5 voxels. Fibres were generated for several different orientations, resulting in 133 individual fibres, as illustrated in Figure 7b.

According to Altendorf and Jeulin 26 , the computed directions (u I P ) undergo an attraction toward each direction t i . As a result, the set of investigated directions must be balanced so that it leads to an isotropic bias. In other words, the sum of the local inertia matrices, as defined in (3), must be proportional to the identity matrix. Different sets of computed directions were used for evaluating the chord lengths in the hedgehog, resulting in n t ranging from 13 to 61. These sets are schematically presented in Appendix B (see Figure 23).

ψ=0 • ψ=90 • θ=0 • 30 • 60 • e 1 e 2 ψ=0 • ψ=90 • θ=0 • 30 • 60 •

Mean orientations

The mean orientation of each fibre of the hedgehog was evaluated according to the procedure detailed in section 2.5.1 with different values of n t . Then, the calculated mean orientations were compared to the theoretical ones, as summed up in Figure 8. In this figure, (ψ, θ) denotes the spherical coordinates (ψ stands for the azimuth angle from e 1 , whereas θ denotes the elevation angles toward e 3 ). In each case, the orientation of fibres lying along the frame directions appear to be accurately evaluated.

Since these directions are parallel to the reference frame (see Figure 23), this is in agreement with the theory introduced by Altendorf and Jeulin 26 . On the contrary, the other directions appear to lead to wrong estimates, specially with n t = 13. With n t = 25, all directions characterized by θ = 30 • or 60 • appear to be well predicted but that with θ = 15, 45 and 75 • lead to largest errors. n t = 37 gives better results for the later orientations without worsening the evaluation of other orientations. When n t = 49 and 61, it is clear that all orientations (except the ones lying along the frame directions) are not accurately evaluated. For the special case θ = 0 • , some orientations appear to be 180 • off the theoretical ones; this is due to round-off errors and the ambiguous azimuth angle for an horizontal fibre ((ψ, θ = 0) ≡ (ψ + π, θ = 0) for non-signed orientations).

Figure 9a represents the mean and the largest angular errors between the computed orientations (a)

ψ=0 • ψ=90 • e 1 e 2 ψ=0 • ψ=90 • (b) ψ=0 • ψ=90 • e 1 e 2 ψ=0 • ψ=90 • (c) ψ=0 • ψ=90 • e 1 e 2 ψ=0 • ψ=90 • (d) ψ=0 • ψ=90 • e 1 e 2 ψ=0 • ψ=90 • (e) ψ=0 • ψ=90 • e 1 e 2 ψ=0 • ψ=90 • Figure 8:
Pole figure of the mean orientations of each fibre of the hedgehog (Figure 7), illustrating the influence of the number of computational directions of the chord lengths: (a) n t = 13, (b) n t = 25, (c) n t = 37, (d) n t = 49 and (e) n t = 61. Blue crosses ( ) represent the computed orientations whereas the red circles ( ) give the theoretical ones, as a comparison. For the sake of readability, the θ axis is not shown (see Fig. 7b for details).

and the theoretical ones for each value of n t . This plot confirms that the evaluation of the orientations of the hedgehog's fibres is better with n t = 37 than with n t = 49 and 61. The same procedure was used on another "hedgehog-like" medium with a nominal radius of 5 vx, evidencing the same trend, as illustrated in Figure 9a (red bars). Since the computation of the inertia matrices runs in linear time with n t , using n t = 37 appears to be the best compromise between execution time and precision.

Preliminary tests, not described in this article, were carried out with hedgehog of various steps. As a result, n t = 37 appears to produce the lowest error regardless the spacing of the hedgehog; thus, all analyses detailed below were done using this value.

Mean radii

Since the orientation of each fibre is well computed, the evaluation of the mean radius (detailed in section 2.5.2) can be done. With n t = 37, the computation of the mean radii of the two aforementioned hedgehogs lead to the distributions given in Figure 9b. For the fibres with a nominal radius of 3 vx, the mean value of r was 3.41 vx (standard deviation: 0.19); for the fibres with a nominal radius of 5 vx, the mean value of r was 4.94 vx (standard deviation: 0.24). Therefore, it appears that the radii are well computed in both cases, leading to slight overestimation for smaller value.

Structuring element

Fibrous networks made up of straight fibres

In order to investigate the influence of the structuring element for the dilation algorithm (i.e. D = D iso vs. D = D in ( 5)) on the resulting geometries, we considered a medium, characterized by two straight fibres with circular cross-section, as illustrated in Figure 10a. These fibres were separated using the separation procedure proposed in this paper, as illustrated in Figure 10b.

After separation of the fibres (Figure 10b), they appear to be well defined, except near the contact zone. As a result, fibres obtained using isotropic dilation are really similar to the ones obtained using longitudinal dilation, as illustrated in Figure 11. Details of the 2D sections illustrate the differences in terms of cross section of each fibre. The isotropic dilation leads to the penetration of the green fibre into the blue one.

Fibrous networks made up of woven fibres

Now let us consider three wavy fibres, with circular cross-section, as illustrated in Figure 12a. They were separated according to the same procedure, leading to the results presented in Figure 12b. The results obtained using isotropic and longitudinal dilation are illustrated in Figure 12c and Figure 12d, respectively. Thus, the differences are evidenced by the remaining orphan voxels (see shaded regions), even after the recursive dilations, when using the longitudinal dilation. Figure 12d shows that the dilation of the green fibre results in a flattened fibre.

Synthesis

As a conclusion, the following statements can be made: longitudinal dilation is recommended for straight fibres because it retains the fibre morphology and leads to better evaluation of the contacts (avoids penetration); isotropic dilation is recommended for wavy fibres because it does not leave orphan voxels and does not flatten the fibres.

Influence of the geometry

In order to assess the robustness of the proposed algorithm against image quality and morphological parameters, two fibres with cylindrical geometry (with radius r = 10 vx), in contact to each other, were generated with varying nominal angle (denoted β below) between their axes. In addition, these geometries were generated for different values of the slenderness ratio λ (ratio of lenght to the diameter).

Then, all these geometries were down-sampled by a factor 2, 3 and 4. In the present section, s denotes the down-sampling factor (s = 1 corresponds to no down-sampling and s = n corresponds to retaining each n-th voxel on each direction). Among all these sets, three geometries are illustrated in Figure 13. Table 1 gives the resulting geometrical parameters (length and diameter) in each case. On each geometry, attempts were made to separate the two fibres using the method detailed in section 2.3. In each case, the largest thresholding value ᾱth leading to separation of the fibres was recorded, as plotted in Figure 14.

For λ ≥ 3, ᾱth appears to be an increasing function of the angle between the fibres. On the opposite, Table 1: Geometrical parameters (radius r and length , in voxels) of the generated fibres, depending on the down-sampling factor s and the slenderness ratio λ.

Down-sampling factor

0 • 10 • 20 • 30 • 40 • ᾱth s = 1 s = 2 0 • 10 • 20 • 30 • 45 • 60 • 90 • 0 • 10 • 20 • 30 • 40 • β ᾱth s = 3 0 • 10 • 20 • 30 • 45 • 60 • 90 • β s = 4 λ = 2 λ = 3 λ = 5 λ = 10
Figure 14: Largest values of ᾱth leading to the separation of the generated fibres, as a function of the angle between the two fibres β, the slenderness ratio λ and the down-sampling factor s. . ᾱth appears to randomly depend on β when λ = 2. This result indicates that the separation process somehow behaves randomly when one attempts to separate short fibres. As a conclusion, the present algorithm seems to work on fibres characterized by λ ≥ 3.

On higher down-sampling factor (s = 4), the separation process can not be performed with β < 20 • , whatever the slenderness ratio is. According to the geometrical parameters given in Table 1, it appears that a radius larger than 2.5 vx is required to ensure the ability of the present algorithm to separate the fibres.

As the proposed algorithm is based on the local change in fibre orientation, the separation of fibres nearly parallel to each other (β ∼ 0 • ) is not possible with the proposed algorithm. In this case, one should use other methods based on skeletonization 8,9,10,11 .

Wavy fibres: influence of the wrap angle

In order to investigate the ability of the algorithm to separate wavy fibres wrapped around each other, theoretical geometries were numerically generated, as presented in Figure 15. In each case, the geometry consisted in two fibres wrapped around each other, with a given wrap angle γ ∈ [0

• , 180 • ],
as detailed in Figure 15a. These geometries were generated for different values of the fibre radius and the wrap angle. In each case, the maximum value of ᾱth resulting in the separation of the two fibres are reported in Figure 15c.

When γ = 0 • , the values of ᾱth appear to be almost constant for a given radius, with mean values equal to 21.7 • , 17.8 • and 15.2 • if r = 2.5 vx, 5 vx and 10 vx, respectively. As a conclusion, this separation method is very efficient on various morphologies of wavy touching fibres regardless the wrap angle.

General guidelines to use the algorithm

Based on the previous results, the following guidelines can be provided to ensure the ability of the proposed approach to accurately separate and characterize the fibres:

n t = 37 investigated directions for the computations of chord lengths appears to be the best compromise between execution time and precision;

each fibre should have a diameter larger than 5 vx and a slenderness ratio larger than 3; longitudinal dilation should be used for straight fibres, since it keeps their straightness and avoids unrealistic inter-penetration; isotropic dilation should be used for wavy fibres because it leads to almost zero orphan voxels at the end of recursive dilatations;

wavy fibres can also be separated, as long as parts of them are straight segments.

No guideline has been raised about the choice for ᾱth . Since it is used as a global parameter (depending on the type of fibrous network), one can simply find the adequate value by trial-and-error.

4 Application: fibre identification in various fibre networks

Considered media

In order to investigate the ability of the proposed algorithm to separate and reconstruct the fibres on complex media, CT images of various materials were used, as summed up in Figure 16. These media were chosen to investigate whether the proposed algorithm can be applied for fibre identification in The fibrous medium (a) was constituted of short pieces of nylon fibres of identical length [START_REF] Laurencin | 3d real-time and in situ characterisation of fibre kinematics in dilute non-newtonian fibre suspensions during confined and lubricated compression flow[END_REF] .

Medium (b) was constituted of short carbon fibres of different lengths embedded in a polymer matrix. Medium (c) was a Sheet Moulding Compound (SMC) made up of planar random glass fibre bundles [START_REF] Sentis | 3d in situ observations of the compressibility and pore transport in sheet moulding compounds during the early stages of compression moulding[END_REF] . Medium (d) was imaged during an in situ bias test of a woven textile (plain weave of metallic fibres with circular cross sections). Medium (e) was an artificial 3D image of a woven textile, where each yarn consisted in a series of slender filaments with circular cross section. This image was generated using the WeaveGeo module as part of the Geodict software [START_REF] Wiegmann | Geodict: A softwarecentered approach to the simulation of thin porous media and their properties[END_REF] .

Medium (a) was imaged on TOMCAT beamline at the Swiss Light Source (SLS). Medium (b) was imaged on ID19 beamline at the European Synchrotron Radiation Facility (ESRF). Media (c) and

(d) were imaged on a laboratory tomograph (3SR Lab, Grenoble, France, RX Solutions apparatus).

The structural properties of each medium are detailed in Table 2. In each case, the threshold value of the misorientation angle (ᾱ th ) used for separation and the kernel used for dilations (D ) are given in this table.

Straight short fibres

Figure 17a represents the fibres identified in the medium composed of nylon fibres (see Figure 16a).

1318 distinct fibres were identified.

The mean orientation of each fibre is presented in the Pole Figure 17b. In this figure, the diameter of each dot is related to the corresponding fibre length. One can see a near transverse isotropic orientation distribution, with almost all fibres nearly lying along the horizontal plane (x, y), with θ < 30 • .

In Figure 17a, one can see that if two fibres are in contact and nearly parallel to each other, they can not be separated. This is in accordance with the cases corresponding to the lower values of β in section 3.3. Figure 18b represents the mean orientation of each fibre (the dot size is related to the corresponding fibre length), emphasizing that the orientation distribution is fully 3D.

ψ=0 • ψ=90 • θ=0 • 30 • 60 • e 1 e 2 ψ=0 • ψ=90 • θ=0 • 30 • 60 • =200 vx

Polydisperse and randomly oriented short fibres

ψ=0 • ψ=90 • θ=0 • 30 • 60 • e 1 e 2 ψ=0 • ψ=90 • θ=0 • 30 • 60 • =200 vx

Sheet moulding compound

Figure 19 shows the identified fibres in the SMC. It appears that most of the fibres were identified despite the random shapes of their cross sections. From these data, the contacts were identified, leading to a coordination number (defined as the mean number of fibre-fibre contact per fibre) equal to 6.57. This value is similar to the ones reported in other SMCs [START_REF] Meyer | Direct bundle simulation approach for the compression molding process of sheet molding compound[END_REF] .

In Figure 19, it appears that the fibres are well-separated despite their random cross sections and their highly anisotropic cross sections. Nevertheless, some fibres have been broken up by the separation, since they appear discontinuous. These splits may be due to the low number of voxels in

Woven plain weave of dense yarns

Attempts were made to separate the tows from the binary image of a textile constituted of filament tows (introduced in Figure 16d), as depicted in Figure 21. After separation, it appears that each Nevertheless, some regions appear to be allocated to wrong fibres near contacts (see circles in Fig-

ure 21). Indeed, because of the filaments, each tow can be divided into subregions. The difficulty of separating such a medium is because it can be considered at different scales: the filament or the tow scale. In addition, because the filaments within each tow are nearly parallel, the proposed algorithm is unable to properly separate them. Thus, other techniques based on morphological analysis [START_REF] Masse | Influence of relative density on the architecture and mechanical behaviour of a steel metallic wool[END_REF][START_REF] Rosini | In situ statistical measurement of local mor-phology in carbon-epoxy composites using synchrotron x-ray computed tomography[END_REF][START_REF] Arthur Hessman | Microstructural analysis of short glass fiber reinforced thermoplastics based on x-ray micro-computed tomography[END_REF] or skeletonization [START_REF] John | The image processing handbook[END_REF][START_REF] Latil | Towards the 3d in situ characterisation of deformation micro-mechanisms within a compressed bundle of fibres[END_REF] may be used in such a context.

Performances

The computational time of the methods proposed in this paper are strongly dependent on the medium. For instance, the computation time of the chord lengths is linear with the number of voxel in Ω f , i.e. linear time with the volume fraction and the volume of the CT image. Since the Least-Misorientation method (detailed in section 2.4.2) requires to consider each fibre at a time, it runs in linear time with n f . This algorithm was developed in Matlab ® , using the Parallel Computing and Image Processing toolboxes. Hence, some routines were parallelized (e.g. local orientations and recursive dilations) in order to speed up the process. In addition, the dilation operations (used in sections 2.4 and 2.5.3) take advantage of bit packing 39 , making the corresponding algorithms very efficient.

Among all media presented in section 4, the image of nylon fibres (Figure 17a) lead to the longest computational time (about 36 h on a 6-core CPU). Lastly, it is worth mentioning that this time could by highly reduced by using Graphics Processing Unit (GPU) computing.

General conclusion

We proposed a method to identify the solid fibres in various fibrous media with volume fraction up to 45%. Starting from binary 3D images, chord lenghts and inertia matrices, as introduced by Altendorf and Jeulin 26 , are computed first; then the misorientation angle is used as a criterion for separation. In order to fully reconstruct the fibres geometries, two dilation algorithms were studied (namely isotropic and longitudinal). As a result, this approach relies on three main parameters: the threshold value on misorientation, the number of considered neighbours for the computation of the chord lenghts, and the dilation operator to be used. All these parameters were analysed, leading to optimum values, depending on the properties of the fibrous network. This approach works well, even on complex media such as:

high volume fraction of fibres, random cross sections, wavy fibres.

The efficiency of the proposed method for separating the fibres was evaluated on elementary cases.

It was shown that a minimum radius of 5 vx and a slenderness ratio larger than 3 is required to ensure the fibre separation between straight fibres in contact to each other. The error made on the radii and the fibres orientations was estimated. For woven fibres, no restriction on the fibre radius was observed.

Since the fibres can be separated, the contact surface between them can be evaluated. In order to get more detailed informations about the contacts, one may mesh the corresponding surface (instead of using the regular grid, inherited from the voxels), for instance using the marching cube method [START_REF] Lorensen | Marching cubes: A high resolution 3d surface construction algorithm[END_REF] , based on the grey level of CT images.

In this paper, the fibres are either considered as wavy or perfectly straight. In a further work, the fibres may be defined piecewise, allowing to investigate the local orientations (such as bending) and local cross sections. This improvement may be in great interest for the dilation procedure. [2] John Summerscales. Non-destructive testing of fibre-reinforced plastics composites, volume 2.

Springer Science & Business Media, 1990.

[3] B London, Robert N Yancey, and JA Smith. High-resolution x-ray computed tomography of composite materials. Materials Evaluation, 48:604-608, 1990.

  Keywords. X-ray tomography, Fibre separation, Contact identification , 3D fibrous materials, Microstructural descriptors Notations α Misorientation angle. ᾱth Threshold value on α. dc P (t i ) Half-chord length in P along t i . D iso Isotropic dilatation operator. D Dilatation operator parallel to ū. e n n-th direction associated to the reference frame. F k Set of voxels belonging to the fibre k. I P Inertia matrix in P ∈ Ω f . λ Slenderness ratio. n t Number of investigated values of t i . n f Number of identified fibres. n orph Number of voxels in Ω orph . Ω f Fibrous phase. Ω f mult Multiple points, resulting from the dilatation operation.

Figure 1 :

 1 Figure 1: Example of a two-phased fibrous material imaged by X-ray tomography: 3-dimensional surface rendering of the binarized image. The fibrous phase is represented as a red solid volume whereas the matrix phase is set to transparent. Dimensions are expressed in voxels.

  90

Figure 2 :

 2 Figure 2: Scheme of the overall proposed algorithm to separate the fibres, reconstruct them, track the fibre-fibre contacts and characterize the fibre geometries. The first block (which aims to separate the fibrous phase from the matrix phase) is not detailed in this work. Despite the fact that all the steps are computed in 3D, the intermediate steps are illustrated in this figure as 2D slices (at z = 80) taken from the example medium (see Fig. 1): (a) threshold on the grey level image, (b) misorientation angle, (c) threshold on the misorientation and results from (d) labelling then (e) reconstruction.

Figure 4 :

 4 Figure 4: Schematic representation of applying either (a) isotropic or (b) longitudinal dilations on aset of voxels (represented in white). The result from dilation is represented in green. In each case, the structuring element[START_REF] Robert M Haralick | Image analysis using mathematical morphology[END_REF] is represented (left-hand side).

Figure 5 :

 5 Figure 5: Evolution of the fraction of orphan voxels among all the fibrous phase, as a function of the iteration step, when applying longitudinal or isotropic dilations.

Figure 7 :

 7 Figure 7: Numerically generated medium: (a) 3D view of the hedgehog-like fibrous phase and (b) pole figure of the mean orientation of each fibre (stereographic projection).

Figure 9 :

 9 Figure 9: Results from the geometrical analyses on the hedgehog: (a) evolution of the mean angular error (darker bars) as a function of n t , depending on the nominal fibre radius (r), lighter bars give the largest angular error; (b) distribution of the fibres' radii evaluated on the hedgehog (with n t = 37), depending on the nominal radius.

Figure 10 :

 10 Figure 10: Two straight cylindrical fibres (copper). (a) Fibrous phase and (b) separated fibres (ᾱ th =15 • ) before dilation.

Figure 11 :

 11 Figure 11: Results after (a) isotropic and (b) longitudinal dilations of the fibres presented in Figure 10. Details of the resulting cross sections are given by the 2D sections (at x = 63).

Figure 12 :

 12 Figure 12: Woven cylindrical fibres: (a) binarized fibrous phase, (b) result from separation with ᾱth =4 • and reconstruction using (c) isotropic or (d) longitudinal dilations. Shaded regions correspond to the remaining orphan voxels at the end of the recursive dilation procedure.

Figure 13 :

 13 Figure 13: Examples of the generated geometries for the investigation of the influence of the image parameters on the separation: (a) β = 3 • , λ = 2, not down-sampled (b) β = 20 • , λ = 5, s = 2 and (c) β = 90 • , λ = 10, s = 4.

Figure 15 :

 15 Figure 15: Separation of wavy touching fibres: (a) schematic view of the synthetic geometry, introducing the wrap angle and (b) example of geometry with γ = 150 • and r = 10 vx; highest values of ᾱth leading to the separation of the fibres wrap around each other, as functions of wrap angle mean fibre radius.

Figure 16 :

 16 Figure 16: 3D images of the investigated media: (a) nylon fibres in candle gel (b) carbon fibres in polymer matrix (c) Sheet Moulding Compound (d) woven textile with monofilament with circular cross-section (e) woven fabric (plain weave) with filament-based yarns (synthetic).

Figure 17 :

 17 Figure 17: Result from separation of nylon fibres: (a) 3D representation and (b) pole figure illustrating the mean orientation of each fibre; a unique colour is attributed to each fibre; dot size is related to the corresponding fibre length.

Figure 18a represents carbon

  Figure 18a represents carbon fibres (introduced in Figure16b) after separation, where 446 fibres were identified. From this figure, it appears that the present algorithm is able to separate randomly oriented fibres (3D distribution), while that proposed by Viguié et al.[START_REF] Viguié | Finding fibres and their contacts within 3d images of disordered fibrous media[END_REF] is only efficient if fibres show a near-planar distribution.

Figure 18 :

 18 Figure 18: Results from separation of carbon fibres: (a) 3D representation and (b) pole figure illustrating the mean orientation of each fibre; a unique colour is attributed to each fibre; the dot size is related to the corresponding fibre length.

Figure 19 :Figure 20

 1920 Figure 19: Results of separation of the fibres in a SMC. A unique colour is randomly attributed to each fibre.

Figure 20 :

 20 Figure 20: Illustration of the contact zones found in the woven fibres, shown in red (fibres are represented in transparency).

Figure 21 :

 21 Figure 21: Identified fibres after separation of the textile composed of filament yarns. The circles show locations where the contacts are not accurately defined.

Figure 22 :

 22 Figure22: Schematic representation of the method used for computing of the chord lengths d j along a given direction t i at a given location (red square): the white squares illustrate the fibrous phase nested in the matrix (black squares).

Figure 23 :

 23 Figure 23: Illustration of the investigated computational directions of the chord lengths from a given voxel: the reference voxel (where the chord lengths are computed) is represented in red and its considered neighbours are represented in blue. (a) n t = 13, (b) n t = 25, (c) n t = 37, (d) n t = 49 and (e) n t = 61.

  ). Consider a voxel P ∈ Ω f ; the distance from P to the matrix along t i , denoted d P (t i ) below, is computed using the method proposed by Altendorf and Jeulin[START_REF] Altendorf | 3d directional mathematical morphology for analysis of fiber orientations[END_REF] . Appendix A provides some details

	2.3 Step 1: Fibres Separation
	2.3.1 Computation of the chord lengths
	Let t i be a given direction and t i 1 , t i 2 and t i 3 its coordinates in the reference frame (t i = t i 1 e 1 + t i 2 e 2 +
	t i 3 e 3 95
	about this method.
	These distances can be used to compute the half-chord length in P , denoted dc P (t i ), as defined
	by Sandau and Ohser 25 :
	∀i = 1, ..., n t dc P t i =	d P (t i ) + d P (-t i ) 2	(1)

Table 2 :

 2 Properties of the media used as applications of the proposed approach, presented in Figure16, and corresponding parameters used for separation (ᾱ th ) and reconstruction (D ).

		Material	Volume fraction	Nominal diameter	ᾱth	D
	(a)	Nylon fibres	33.4 %	16 vx	4 • Longitudinal
	(b)	Carbon fibres	6.12 %	11 vx	4 • Longitudinal
	(c)	SMC	25.9 %	-	4 • Longitudinal
	(d)	Woven fibres	5.76 %	10 vx	10 •	Isotropic
	(e)	Woven tows of fibres	43.1 %	300 fib./tow, fib.: 4 vx 8 •	Isotropic
	(a)			(b)

The two types of kernels used in this paper, illustrated in Figure4, are described bellow.

Acknowledgments

This work was supported by the Michelin Corporation, 23 place Carmes Déchaux, 63000 Clermont Ferrand, FRANCE through a research contract, as part of the PHARE project. The 3SR Lab is part of LabEx Tec 21 -ANR-11-LABX-0030 and of Institut Carnot PolyNat (ANR16-CARN-0025).

Author's contributions

The contributions of the authors are the following: L.O, B.L. and S.R. conceived the study; D.D., L.O. and S.R. discussed the algorithms to de developed and the validation tests to be carried out to identify fibres and fibre-fibre contacts; D.D. developed the code; D.D. and S.R. wrote the manuscript; all the authors analysed the results and revised the manuscript before submission.

Appendix A: Computation of the chord lengths

In this paper, the method proposed by Altendorf and Jeulin 26 have been used to compute the chord lengths at voxel P , denoted d P here. Figure 22 schematically illustrates how d j (t i ) and d j (-t i ) are defined.

Appendix B: Investigated directions for computation of the chord lenghts

According to the definition of the half-chord lenght dc P (t i ), given in (1), it is clear that dc P (-t i ) = dc P (t i ). This means that the computation of I P based on the N-connected neighbourhood can be reduced to n t = N/2 directions. In addition, because of the way the chord lenghts are computed dc P (t i ) (see for instance Figure 22), it appears that: ∀k, dc P kt i = dc P t i .

In other words, if dc P (t i ) has already been computed, the computation of the half-chord length along any direction parallel to t i would be redundant for the evaluation of I P . Thus, all sets of 410 computational directions t i used for computing the chord lengths, as investigated in section 3.1, are presented in Figure 23.