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

Marouane Kabbej, Valérie Guillard, Hélène Angellier-Coussy, Valentin

Thoury-Monbrun, Nathalie Gontard, Laurent Orgéas, Sabine Rolland Du

Roscoat, Sébastien Gaucel

► To cite this version:

Marouane Kabbej, Valérie Guillard, Hélène Angellier-Coussy, Valentin Thoury-Monbrun, Nathalie Gontard, et al. Cascading (3D) reconstruction procedure of composite structures from microtomography data. MethodsX, 2023, 10, pp.102177. 10.1016/j.mex.2023.102177. hal-04094129

HAL Id: hal-04094129

https://hal.science/hal-04094129

Submitted on 10 May 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Cascading (3D) reconstruction procedure of composite structures from microtomography data

<figure><img src="image_2.png" /></figure>

Check for
updates

Marouane Kabbejª, Valérie Guillardª, Hélène Angellier-Coussyª,
Valentin Thoury-Monbrunª, Nathalie Gontardª, Laurent Orgéasᵇ, Sabine Rolland du
Roscoatᵇ, Sébastien Gaucelª,*

## ARTICLE INFO

**Method name:**

Cascading (3D) reconstruction procedure of
composite material

**Keywords:**

Particle morphology
(3D) Structure building
Heterogeneous size distribution
(3D) tri-phasic structure

## ABSTRACT

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.

<footnote>ªIATE, CIRAD, INRAE, Université Montpellier, Institut Agro, Montpellier 34060 France</footnote>

<footnote>ᵇCNRS, 3SR Lab, Université Grenoble Alpes, UMR 5521, Grenoble F-38000, France</footnote>

<footnote>* Corresponding author.
E-mail address: sebastien.gaucel@inrae.fr (S. Gaucel).</footnote>
Specifications table

<table><tr><td>Subject Area:</td><td>Materials Science</td></tr><tr><td>More specific subject area:</td><td>(3D) Numerical modelling</td></tr><tr><td>Method name:</td><td>Cascading (3D) reconstruction procedure of composite material</td></tr><tr><td>Name and reference of original method:</td><td>The procedure is an original one and combined different standard/published techniques and novel ones. The standard/published techniques are:<br/>- Tschopp MATLAB code [7]<br/>- (3D) Objects Counter" in ImageJ software [4]<br/>- MorphoLibJ plugin in ImageJ software [2]<br/>- Conversion expressions of Tait-Bryan angles into Euler angles [1]</td></tr><tr><td>Resource availability:</td><td>Link to dataverse:<br/>Size distribution of cellulose particles into PP/cellulose particles<br/>https://doi.org/10.57745/PDQNR6<br/>Descriptors of the cellulose particles detected in the (3D) microtomography images of composites PP/cellulose<br/>https://doi.org/10.57745/KXPWXK</td></tr></table>

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³/voxel. From one (3D) tomographic image of real dimension of 1844 × 1403 × 298 µm³ (length x depth x height), 12 cuboid samples of same height, 250 µm, were extracted from different zones (Fig. 1a) 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. 1b). This allows the detection and identification of particles. Afterwards, observed particles were assimilated to ellipsoids (Fig. 1c) 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. 1c-1) and orientation in the composite material, given by three orientation angles ($\psi, \theta, \phi$) (Fig. 1c-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 \in [0, L_x] \times [0, Ly] \times [0, L_z]$, 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 $\varphi_{p, target}$ 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

$$|\varphi_{p, target} - \varphi_p| \le 0.01\%v/v.$$

Then, these particles are sequentially positioned and orientated in the periodic RVE following the sequential approach described below:
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Images analysis steps of the (3D) microtomography: (a) Extraction of 12 cuboid samples from the global image of dimension of 1844 × 1403 × 298 µm³;(b) Binarization of each of the 12 images of dimension reported in 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.</figcaption></figure>

i. The particle's position (i.e. centre coordinates $x_c, y_c, z_c$) 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. 1c-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 ($z=0$ and $z=L_z$) 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. 3a).

### 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. 3b).

### Method validation

**Step 1.** To test the cascading (3D) reconstruction procedure proposed here, a composite of cellulose (= the particles or filler), dispersed into Polypropylene (= the matrix) with a filler volume fraction of $φ_p = 2.96\%v/v$ was considered ([5]a). The (3D)
<figure><img src="image_4.png" /><figcaption>Fig. 2. Numerical process to generate the (3D) microstructure of the RVE performed in MATLAB.</figcaption></figure>

tomographic image of dimension of 1844 × 1403 × 298 µm³ was transformed into 12 cores cuboid samples cumulating a total of 21,570 particles (Fig. 1a) (see S1 in Support Material).

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, $\varphi_p = 2.96\%v/v$.

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. 1c). The ellipsoids of inertia, determined by the MorphoLibJ plugin are characterised by the following nine coefficients,
* Coordinates of the centroid $(x_c, y_c, z_c)$

* Axes $a_1, a_2, a_3$, corresponding to the length of the largest, medium and smallest axis, respectively ($a'_1 > a'_2 > a'_3$)

* Tait-Bryan angles $(\psi, \theta, \phi)$ for orientation with x-y-z rotation convention (extrinsic rotations)

<figure><img src="image_5.png" /><figcaption>Fig. 3. Representation of the RVE showing structures with size-heterogeneous ellipsoidal particles: (a) bi-phasic structures and (b) tri-phasic structures of interphase thickness of 2 µm (purple colour), both generated using the algorithm proposed in that work. The structures are generated for particles volume fraction corresponding to true composite materials ($\phi_p = 2.96 - 6.06 - 12.67 - 19.91\%v/v$).</figcaption></figure>

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 ($a_1', a_2', a_3'$) were normalised and converted as follows:

$$a_i = a_i' \left( \frac{V_{\text{real}}}{V_{\text{ell}}} \right)^{\frac{1}{3}} \quad \text{for } i \in \{1, 2, 3\} \qquad (1)$$

where $V_{\text{real}}$ and $V_{\text{ell}}$ are the volume of the real particle and the ellipsoid of inertia ($V_{\text{ell}} = \pi/6 \times a_1' \times a_2' \times a_3'$), respectively.

MorphoLibJ provides Tait-Bryan angles $(\psi, \theta, \phi)$ 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 $\psi$ about x-axis, the second is by the elevation angle $\theta$ about y-axis and the third is by the azimuth angle $\phi$ about z-axis.

The variation intervals of the angles $\psi, \theta$ and $\phi$ 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 $(\psi, \theta, \phi)$. Then, we used some properties of symmetry to reduce the variation intervals of both the angles $\psi$ and $\phi$ to [0°,180°], so that each particles orientation corresponds to a single combination of the triplet $(\psi, \theta, \phi)$.

For the upcoming numerical simulations, it is necessary to convert the Tait-Bryan angles $(\psi, \theta, \phi)$ obtained from the MorphoLibJ plugin of ImageJ to the classical Euler angles $(\gamma, \beta, \alpha)$ 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 $\gamma$ about the z-axis, the second is by nutation angle $\beta$ about the x-axis and the third is by precession angle $\alpha$ 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]).

The experimental particles sizes ($a_1, a_2, a_3$) found from particles' analysis were $a_1=11.57 \pm 9.87$, $a_2=6.83 \pm 4.77$ and $a_3=4.35 \pm 2.43$ µm (mean values of all particles).

**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:

$$f(x) = \frac{1}{N \cdot h} \sum_{i=1}^{N} K\left(\frac{x - x_i}{h}\right) \qquad (2)$$
<figure><img src="image_6.png" /><figcaption>Fig. 4. Scatter plots of axes a₁>a₂>a₃ 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₁ and a₂ (column a), a₁ and a₃ (column b) and a₂ and a₃ (column c) are presented. Raw data available at https://doi.org/10.57745/PDQNR6.</figcaption></figure>

where x₁, x₂, ..., xₙ 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):

$$K(x) = \frac{1}{2\sqrt{\pi}} e^{-\frac{1}{2}x^2} \qquad (3)$$

The experimental distribution of the largest axis a₁, the axes ratios a₁/a₂ and a₂/a₃ were fitted using the truncated Kernel law. The use of triplet (a₁, a₁/a₂, a₂/a₃), instead of direct fitting of (a₁, a₂, a₃), was motivated by the dependency observed between axis in data from image analysis. To preserve the correlation between a₁, a₂ and a₃ on the whole a₁ value range it was necessary to use four Kernel distributions for a₁/a₂ defined in four separate zones of a₁ (zones bounded by red dashed lines in Fig. 4-Ia): a₁ ≤ 20 µm; 20 < a₁ ≤ 40 µm; 40 < a₁ ≤ 60 µm and a₁ > 60 µm.

Fig. 4 shows the real correlations between the axes a₁, a₂ and a₃, 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₁/a₂ are generated according to a single Kernel distribution defined on the entire a₁ zone while the blue bullets present an example of generated axes values for 15,000 particles where a₁/a₂ is generated according to four Kernel distributions defined in four separate zones of a₁. The experimental correlations between a₁, a₂ and a₃ (grey bullets) are better described in the case where a₁/a₂ is generated in four separate zones of a₁ (blue bullets) contrary to the case where a₁/a₂ is generated on the entire
<figure><img src="image_7.png" /><figcaption>**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 $a_1/a_2$ in the zone $a_1 \le 20 \ \mu\text{m}$ (b-1), $20 < a_1 \le 40 \ \mu\text{m}$ (b-2), $40 < a_1 \le 60 \ \mu\text{m}$ (b-3) and $a_1 > 60 \ \mu\text{m}$ (b-4) obtained from 10,945, 1596, 219 and 52 particles, respectively.</figcaption></figure>

$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 \le 20 \ \mu\text{m}$ (b-1); $20 < a_1 \le 40 \ \mu\text{m}$ (b-2); $40 < a_1 \le 60 \ \mu\text{m}$ (b-3) and $a_1 > 60 \ \mu\text{m}$ (b-4).

The Euler angles $(\psi, \theta, \phi)$ 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 \times 50 \times 100 \ \mu\text{m}^3$) of the composite Polypropylene/cellulose containing size-heterogeneous ellipsoidal particles with $\phi_p = 2.96\%$ were generated using the
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Comparison between experimental distribution and the fitted truncated Kernel distribution for (a) roll angle ψ (b) elevation angle θ and (c) azimuth angle φ obtained from 12,812 particles.</figcaption></figure>

procedure described in Figs. 1 and 2. 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.

**References**

[1] P. Berner, Technical concepts: Orientation, rotation, velocity and acceleration, and the SRM, SEDRIS 64 (2008) https://docslib.org/doc/12571458/orientation-rotation-velocity-and-acceleration-and-the.

[2] D. Legland, I. Arganda-Carreras, MorphoLibj User Manual. https://imagej.net/imagej-wiki-static/MorphoLibJ, 2020.

[3] D. Legland, I. Arganda-Carreras, P. Andrey, MorphoLibJ: integrated library and plugins for mathematical morphology with ImageJ, Bioinformatics 32 (July) (2016) 3532–3534, doi:10.1093/bioinformatics/btw413.

[4] J. Schindelin, I. Arganda-Carreras, E. Frise, V. Kaynig, M. Longair, T. Pietzsch, S. Preibisch, C. Rueden, S. Saalfeld, B. Schmid, J.Y. Tinevez, D.J. White, V. Hartenstein, K. Eliceiri, P. Tomancak, A. Cardona, Fiji: an open-source platform for biological-image analysis, Nat. Methods 9 (7) (2012) 676–682, doi:10.1038/nmeth.2019.

[5] V. Thoury-Monbrun, Thèse, Université de Montpellier, Thèse, université de Montpellier, 2018.

[6] V. Thoury-Monbrun, H. Angellier-Coussy, V. Guillard, D. Legland, S. Gaucel, Impact of two-dimensional particle size distribution on estimation of water vapor diffusivity in micrometric size cellulose particles, Materials (Basel) 11 (9) (2018), doi:10.3390/ma11091712.

[7] M.A. Tschopp (2016). 3-D synthetic microstructure generation with ellipsoid particles. 3-D synthetic microstructure generation with ellipsoid particles. Report number: ARL-TN-0789 of Army Research Laboratory, September. doi:10.13140/RG.2.2.19120.66560.