Cascading (3D) reconstruction procedure of composite structures from microtomography data

Reconstruction of three-dimensional (3D) structure from experimental image acquisition (e.g., from micro computed tomography data) is very useful in composite material science. Composite considered are characterized by a dispersion of particles in a continuous phase. Many properties of the composite (e.g., mass transfer properties) depend on its structural assembly. A reliable prediction of these properties requires to well represent this structure and especially, the region at the vicinity of the dispersed phase. (3D) structure generation must thus permit to (1) simplify the real composite structure observed to make it compatible with further modelling tasks (e.g., meshing constraints in finite elements methods, computation time) and (2) keep enough representativeness of the structure of the specimen to produce reliable numerical predictions. This article describes an innovative, cascading (3D) reconstruction procedure of composite material from microtomography data.• First step of this pipeline is the extraction of relevant structural markers from microtomography images using image analysis.• Second step is the modelling of the distribution of the structural markers selected (statistical laws).• Third and final step is the reconstruction of the (3D) structures based on the pre-determined distribution laws in a RVE (representative volume element) of the composite.

Reconstruction of three-dimensional (3D) structure from experimental image acquisition (e.g., from micro computed tomography data) is very useful in composite material science. Composite considered are characterized by a dispersion of particles in a continuous phase. Many properties of the composite (e.g., mass transfer properties) depend on its structural assembly. A reliable prediction of these properties requires to well represent this structure and especially, the region at the vicinity of the dispersed phase. (3D) structure generation must thus permit to (1) simplify the real composite structure observed to make it compatible with further modelling tasks (e.g., meshing constraints in finite elements methods, computation time) and (2) keep enough representativeness of the structure of the specimen to produce reliable numerical predictions. This article describes an innovative, cascading (3D) reconstruction procedure of composite material from microtomography data.
• First step of this pipeline is the extraction of relevant structural markers from microtomography images using image analysis. • Second step is the modelling of the distribution of the structural markers selected (statistical laws). • Third and final step is the reconstruction of the (3D) structures based on the pre-determined distribution laws in a RVE (representative volume element) of the composite.

Method details
Step 1. (3D) microtomography images analysis For this work, we used (3D) grayscale tomographic image of the specimen analysed. Tomographic images used were obtained using X-ray microtomography on the ID19 beamline at the ESRF in Grenoble with a resolution of 0.65 μm 3 /voxel. From one (3D) tomographic image of real dimension of 1844 × 1403 × 298 μm 3 (length x depth x height), 12 cuboid samples of same height, 250 μm, were extracted from different zones ( Fig. 1 a) to be further analysed (precise dimensions of the twelve samples are given in Table S1.1 of suppl. material).
The analysis time of the global image is considerable (several hours or even days), this is the reason we have chosen to analyse 12 cuboid samples of the global image that reduce considerably the computation time (about 10-15 min for each sample: about 3 h for all the 12 samples). Note that images analysis was performed in a DELL computer with Intel Xeon E-2176 M Processor (2.7 GHz) and 32 Gb of Ram.
Using the option "(3D) Objects Counter'' in ImageJ software [4] , segmentation (binarization) was first applied to these 12 cuboid samples by changing grayscale images to binary images: 0 for the matrix and 1 for the fillers ( Fig. 1 b). This allows the detection and identification of particles. Afterwards, observed particles were assimilated to ellipsoids ( Fig. 1 c) whose shape descriptors were measured using MorphoLibJ plugin in ImageJ software [2 , 3] .
The distribution of these descriptors is the basis of the generated (3D) structures in this work.
In the present work, the particles are considered as ellipsoids and are characterised by morphological parameters, the three axes a 1 , a 2 and a 3 (a 1 > a 2 > a 3 ) ( Fig. 1 c-1) and orientation in the composite material, given by three orientation angles ( , , ) ( Fig. 1 c-2) which are further detailed.
Step 2. Structural parameters distribution laws Distribution of descriptors identified in step 1 is then modelled by using a statistical distribution law. The best fit obtained by the statistical tools in Matlab library permits to select the appropriate distribution law.
Step 3. (3D) Structure generation The MATLAB code developed for generating our (3D) structures was based on Tschopp MATLAB code [7] . The code generates (3D) microstructures composed of a population of non-overlapping ellipsoidal particles heterogeneously distributed in size and orientation, within a periodic RVE.
Organisational chart summarising the structure of the generation algorithm could be found in Fig. 2 . Bi-phasic or tri-phasic systems (i.e., obtained by adding a third compartment surrounding the dispersed particles) can be generated with this generation procedure.

Bi-phasic system
The bi-phasic composite structure is generated in a cuboid shaped representative volume element (RVE) defined by x,y,z ∈[0, where z is the direction orthogonal to the permeation flux for a composite membrane. The RVE is supposed periodic along its vertical faces, to represent an infinite repetitive structure along x and y axis.
Structure generation requires the RVE size and the target particles volume fraction , in% volume/volume (%v/v) as inputs. The geometric parameters of a population of particles, a 1 , a 2 and a 3 , are randomly generated using distribution law determined in step 2 and satisfying the target volume fraction of particles, using the criterion Then, these particles are sequentially positioned and orientated in the periodic RVE following the sequential approach described below:  Table S1.1 of suppl. material. Here is an example of binarization of the sample N°5. In the binary image, the matrix and the particles are in white and black colour respectively; (c) Modelling a particle as an ellipsoid in (3D) space. (c-1) a 1 : largest axis, a 2 : medium axis, a 3 : smallest axis (c-2) : roll angle, : elevation angle, : azimuth angle. All descriptors obtained from the 12 images analysed are available at https://doi.org/10.57745/KXPWXK .
i. The particle's position (i.e. centre coordinates , , ) is randomly generated using uniform distribution. Note that the ellipsoidal particles are sorted by decreasing largest axis a 1 , large particles being positioned first. At the same time, the orientation angles, i.e. roll , elevation and azimuth angles ( Fig. 1 c-2), are randomly generated using the appropriate distribution identified in step 2. ii. The non-overlapping of the particle with the horizontal faces of the RVE ( = 0 and = ) and with the existing particles is tested. If the non-overlapping tests are successful, the particle is added to the structure and then the next particle is considered for steps i-ii. If at least one of the non-overlapping conditions is not satisfied, then a new position is drawn, orientation being unchanged, and step ii is performed again on the updated particle. This last sequence is repeated until the particle is added to the structure. It should be noted that if a particle intercepts one of the vertical faces of the RVE then, the particle section outside from the RVE is shifted to the opposite face, in order to ensure the periodicity of the RVE ( Fig. 3 a).

Tri-phasic system
The tri-phasic composite structure is built by considering an interphase of fixed thickness around each particle of the already generated bi-phasic structures. The interphase volume is thus different for each particle keeping constant the initial particle volume. For the interphase compartment, the periodicity at vertical faces of the RVE still applies. For particles close to the upper (resp lower) face of the RVE, the interphase may intercept the boundary. In that case, the intercept with the exterior of the RVE is removed ( Fig. 3 b).
To detect the particles, different thresholding values of the segmentation (binarization) were tested in all the 12 cores of cuboid samples. To validate it, volume fraction of the particles on the binarized cuboids were calculated and compared to the experimental one. Globally, the more the thresholding increases, the more the volume fraction of the particles decreases. In the present case, the thresholding of 105 was selected as it permits to achieve a calculated volume fraction very close to the experimental one, = 2 . 96% ∕ . Before analysing the 12 cores of cuboid samples to extract the geometrical descriptors, the particles intercepting the edges of the cores were removed. Also, the particles with less than 30 voxels were rejected, to avoid meshing errors in the upcoming numerical simulations, reducing the total number of particles from 21,570 to 12,812. These rejected particles represent 40.60% in number and 0.47% in volume of the 21,570 particles, respectively. Despite their large number, these tiny, rejected particles represent less than 1% of the volume of all the particles, which is low. Theoretically, these tiny particles may be fragments of cellulose which, due to their large number, generate a large volume of interphase. However, no such particles were visible in 2D images of cellulose particles as observed and described in [6] .
Finally, geometrical descriptors of the 12,812 remaining particles were measured using MorphoLibJ plugin in ImageJ software (see Section 6.2 of the ref. [2] ) assimilating particles to ellipsoids ( Fig. 1 c). The ellipsoids of inertia, determined by the MorphoLibJ plugin are characterised by the following nine coefficients, Ellipsoid of inertia has the same first-order and second-order moments than the studied real particle i.e. same centre of inertia and orientation but not the same volume. To preserve the real volume, the axes obtained with the MorphoLibJ plugin ( 1 ′ , 2 ′ , 3 ′ ) were normalised and converted as follows: where and are the volume of the real particle and the ellipsoid of inertia ( = ∕6 × 1 ′ × 2 ′ × 3 ′ ) , respectively. MorphoLibJ provides Tait-Bryan angles ( , , ) with x-y-z extrinsic rotation convention. It corresponds to a succession of three rotations about the space-fixed principal axes , namely the x, y and z axes. The first rotation is by the roll angle about x -axis, the second is by the elevation angle about y -axis and the third is by the azimuth angle about z -axis.
The variation intervals of the angles , and considered by the MorphoLibJ plugin are [ − 180°,180°], [ − 90°,90°] and [ − 180°,180°] respectively. With these intervals of angles, we encounter cases where particles are orientated in a similar way but with different combinations of the triplet ( , , ) . Then, we used some properties of symmetry to reduce the variation intervals of both the angles and to [0°,180°], so that each particles orientation corresponds to a single combination of the triplet ( , , ) .
For the upcoming numerical simulations, it is necessary to convert the Tait-Bryan angles ( , , ) obtained from the MorphoLibJ plugin of ImageJ to the classical Euler angles ( , , ) with z-x-z extrinsic rotation convention that COMSOL software uses as orientation inputs of particles. The used convention of the Euler angles corresponds to a succession of three rotations about the space-fixed principal axes (extrinsic rotations), namely the x and z axes (z-x-z rotation convention). The first rotation is by intrinsic angle about the z-axis, the second is by nutation angle about the x-axis and the third is by precession angle about the z-axis. The conversion expressions are detailed in the section S2 of suppl. material (see also Table 1 in page 19 of the reference of [1] ).
Step 2 . To randomly generate particle descriptors selected in the previous step (i.e., sizes and orientation angles), their distribution must be first modelled. In the current case study, a non-parametric Kernel distribution model (proposed by MATLAB) was used to generate descriptors' distribution as it was identified as the best law to describe the experimental distribution for the selected descriptors.
Kernel distribution is a non-parametric representation of the probability density function (PDF) of a random variable x. A kernel distribution is defined by a smoothing function K(x) and a bandwidth value h, which control the smoothness of the resulting density curve: Fig. 4. Scatter plots of axes a 1 > a 2 > a 3 of ellipsoid-shaped particles from three different sources: (grey bullets) values from 12,812 particles analysed in a X-ray tomographic image, (red bullets) values from 15,000 generated particles considering a single Kernel distribution for the ratio a1/a2, (blue bullets) values from 15,000 generated particles considering four Kernel distributions for the ratio a1/a2, defined in four intervals for a1 described by dashed red lines in I-a. The correlations between a 1 and a 2 (column a), a 1 and a 3 (column b) and a 2 and a 3 (column c) are presented. Raw data available at https://doi.org/10.57745/PDQNR6 . where x 1 , x 2 , …, x n are random samples from an unknown distribution, N is the sample size, K(x) is the kernel smoothing function, and h is the bandwidth. K is often chosen as the density of a standard Gaussian function (mean μ= 0; Std. dev. = 1): The experimental distribution of the largest axis a 1 , the axes ratios a 1 /a 2 and a 2 /a 3 were fitted using the truncated Kernel law. The use of triplet (a 1 , a 1 /a 2 , a 2 /a 3 ), instead of direct fitting of (a 1 , a 2 , a 3 ), was motivated by the dependency observed between axis in data from image analysis. To preserve the correlation between a 1 , a 2 and a 3 on the whole a 1 value range it was necessary to use four Kernel distributions for a 1 /a 2 defined in four separate zones of a 1 (zones bounded by red dashed lines in Fig. 4 -Ia): a 1 ≤ 20 μm; 20 ⟨ a 1 ≤ 40 μm; 40 < a 1 ≤ 60 μm and a 1 ⟩ 60 μm. Fig. 4 shows the real correlations between the axes a 1 , a 2 and a3, for the 12,812 particles (grey bullets) analysed from the 12 cuboid samples extracted from the tomographic image. The red bullets present an example of generated axes values for 15,000 particles when a 1 /a 2 are generated according to a single Kernel distribution defined on the entire a 1 zone while the blue bullets present an example of generated axes values for 15,000 particles where a 1 /a 2 is generated according to four Kernel distributions defined in four separate zones of a 1 . The experimental correlations between a 1 , a 2 and a 3 (grey bullets) are better described in the case where a 1 /a 2 is generated in four separate zones of a 1 (blue bullets) contrary to the case where a 1 /a 2 is generated on the entire Fig. 5. Comparison between the experimental distribution and the fitted truncated Kernel distribution for (a) a 1 and (c) a 2 /a 3 obtained from 12,812 particles. Comparison between experimental distribution and the fitted truncated Kernel distribution for a1/a2 in the zone a 1 ≤ 20 μm (b-1), 20 < a 1 ≤ 40 μm (b-2), 40 < a1 ≤ 60 μm (b-3) and a 1 > 60 μm (b-4) obtained from 10,945, 1596, 219 and 52 particles, respectively. a 1 zone (red bullets). As shown in Fig. 4 , the preservation of the correlation between a 1 & a 2 (column a) and a 2 & a 3 (column c) tends to also preserve the correlation between a 1 & a 3 (column b). Fig. 5 presents the comparison between the experimental distribution and the fitted truncated Kernel distribution for the axis a 1 (a) and a 2 /a 3 (c). The same comparison is performed for a 1 /a 2 (b) in the zone a 1 ≤ 20 μm (b-1); 20 < a 1 ≤ 40 μm (b-2); 40 < a 1 ≤ 60 μm (b-3) and a 1 > 60 μm (b-4).
The Euler angles ( , , ) were randomly generated according to the truncated Kernel distribution fitted to their real distribution as shown in Fig. 6 .
For other simulation needs, the experimental distribution of the diameter d of the equivalent spheres, in volume, of the real particles could be also fitted using the truncated Kernel law (see Fig. S3.1 of suppl. material).
Step 3 . Based on the Kernel distribution previously identified, about 10 (3D) microstructures ( 50 × 50 × 100 3 ) of the composite Polypropylene/cellulose containing size-heterogeneous ellipsoidal particles with = 2 . 96% were generated using the Higher volume fraction (up to 19.91%) could be also generated without too much impact on computation time ( Fig. 4 ).
Structures with size-heterogeneous spherical particles could be also generated ( Fig. S3.2 of suppl. material).

Funding
This project has received funding from the European Union's Horizon 2020 research and innovation program under grant agreement No 773375. 3SR is part of the labex Tec 21 (ANR-11-LBX-0030) and of Institute Carnot Polynat (ANR-16-CARN-0025-01).

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.

Data availability
I have shared the link toward raw data via DOI (see in the manuscript)

Supplementary materials
Supplementary material associated with this article can be found, in the online version, at doi: 10.1016/j.mex.2023.102177 .