Tessellation-based stochastic modelling of 3D coating structures imaged with FIB-SEM tomography

To facilitate printing, coatings are typically applied to paperboard used for packaging to provide a good surface for application. To optimise the performance of the coating, it is important to understand the relationship be-tween the microstructure of the material and its mass transport properties. In this work, three samples of paperboard coating are imaged using combined focused ion beam and scanning electron microscope (FIB-SEM) tomography data appropriately segmented to characterise the internal microstructure. These images are used to inform a parametric, tessellation-based stochastic three-dimensional model intended to mimic the irregular geometry of the particles that can be seen in the coating. Parameters for the model are estimated from the FIB-SEM image data, and we demonstrate good agreement between the real and virtual structures both in terms of geometrical measures and mass transport properties. The development of this model facilitates exploration of the relationship between the structure and its properties.


Introduction
Coatings consisting of clay, calcium carbonate, and a latex binding material are typically applied to paperboard used for packaging to provide a good surface for printing.The aim is to keep the ink particles inside the coating and to transport the solvent through the paperboard.
To optimise the performance of the coating, it is therefore important to understand and characterise the relationship between the microstructure of the material and its mass transport properties [1,2].
Analysis of the properties of the coating structure requires the material to be well imaged and characterised.However, given the relatively small size of the pores, standard approaches such as X-ray computed tomography are not viable as they will not give us the resolution we require.To address this, we instead used combined focused ion beam and scanning electron microscope (FIB-SEM) tomography [3][4][5], whereby a three-dimensional volume is imaged as a series of thin twodimensional cross section images acquired sequentially.By gradually milling away thin slices of the sample using the FIB, a series of planar cross sections can be imaged using the SEM, resulting in high resolution two-dimensional slice images that enable us to visualise the individual particles that comprise the coating and together form a threedimensional representation of the material.
After processing the three-dimensional image data to segment the structure into a binary array and identify the pore space, simulations of mass transport can be performed directly on the imaged structure.However, preparing and imaging samples with different properties can be very time consuming and expensive.A solution, therefore, is to develop a realistic parametric model and generate a three-dimensional virtual coating structure that is intended to mimic the properties of the material.Exploring virtual structures by characterising their geometry and simulating their mass transport properties is much cheaper and faster than performing real experiments, and a good complement to the real image data.Stochastic modelling of material structures is an extensively researched field of study [6], and has been used in, among others, the generation of fiber structures [7,8], foam materials [9,10], and the design of batteries [11].
In this work, three samples of paperboard coating were imaged using FIB-SEM tomography to visualise the internal microstructure, with the data then segmented to classify each pixel as solid or pore.These threedimensional structures were analysed to give estimates of the parameters of the model [12].The virtual structure is generated in three steps.First, a sphere packing is generated informed by the analysis of the real material.We then generate a Laguerre tessellation around these spheres to capture the fractured geometry of the particles that comprise the coating.Finally, we build our pore structure around this skeleton using an appropriately smooth Gaussian random field.The real and virtual structures are compared using several statistical measures [13], together with mass transport simulations through each structure [14], and we demonstrate good agreement with all three imaged samples.The model can also be used to explore the parameter space and simulate different material structures to understand the impact these parameters have on the structure and its properties.
The paper is organised as follows: first, we present details of the imaging and segmentation and show the segmented images for the three coating samples we will use for our analysis.Then, we give details of each step of the stochastic model and present our characterisation analysis of the structures from which we estimate the parameters of the model.Finally, we detail the statistical measures we use to compare the real and virtual structures.

Paperboard sample imaging
We consider three samples of paperboard coating with different ratios of pigment and binder, which are specified at the creation of the material.The paperboard in question is a liquid packaging board used to make milk packages, and has only a minor impact on the coating itself.The pigment used is Hydrocarb 90 (HC90), a commercial calcium carbonate powder with relatively small particles, and the binder a latexbased binding material.The samples were imaged using combined focused ion beam and scanning electron microscope (FIB-SEM) tomography.In this technique, thin slices of material are sequentially removed using a focused ion beam and each exposed surface is then imaged using the electron beam of the scanning electron microscope.The high resolution images of the slices are combined to form a three-dimensional structure of the specimen material.Prior to imaging, small pieces of paperboard (around 1 cm 2 ) were cut out of a larger sheet using a scalpel and mounted onto aluminium stubs using adhesive carbon tape.The edges of the pieces were additionally fixed using conductive silver paint.The surface of each sample was sputter-coated with a thin layer of palladium for increased electrical conductivity using a Emitech K550X sputter coater.
The FIB-SEM imaging was performed using a Tescan GAIA3 FIB-SEM with a gallium ion column.An ion beam energy of 30 keV and electron beam energy of 1.55 keV were used.The pixel size in each image was 10 nm, and the thickness of each removed slice was 20 nm.A backscattered electron detector was used for imaging, and prior to reconstruction of the three-dimensional structure the images were aligned using Fiji and tomviz to compensate for sample movement [15,16].The material details for each sample are shown in Table 1, and the output from one slice for each coating sample is shown in Fig. 1(a).From observing these slices, we can see the fractured and irregular geometry of the particles that comprise the coating.1, which also includes the data for each sample.

Segmentation
Segmentation, i.e. identification of the porous network, is a challenge for FIB-SEM data and requires more advanced methods.For this reason, segmentation of the image data was performed in the same manner as in [4], and which we describe briefly here.Manual segmentation was performed for each of the data sets in 50 randomly selected two-dimensional regions (of size 128 × 128 pixels) by an expert.In order to extract intensity information at different spatial scales, a set of Gaussian smoothing filters were applied to the data (commonly referred to as linear scale-space features).Gaussian filters with standard deviations σ = 1, 2, 4, 8, 16, 32, 64, and 128 pixels were applied, and the original data (corresponding to σ = 0) was also used.These filters were applied not only to the slice to be segmented, but also to 5 adjacent slices in each direction.For each data set, 30 regions were used for training, 10 for validation, and 10 for testing.
A random forest classifier was trained to classify individual pixels as either solid or pore, and its parameters were optimised with respect to the classification performance on the validation data.After developing the algorithm based on the manually segmented part of the data, it was applied to the full data sets, and the final segmentation was smoothed to reduce classification noise.The classification accuracies on the test data for the data sets A,B, and C were 91.9%,91.9%,and 92.7%.The code was developed in-house using Matlab [17].One slice from each of the segmented structures is shown in Fig. 1(b), with the full threedimensional reconstruction of the coating structure shown in Fig. 1(c).

Tessellation model
The construction of the virtual three-dimensional coating structure is performed in three steps; first, the cell structure representing the solid material is generated using a random sphere packing method.The radii of the spheres are sampled from an estimate of the cell-size distribution of the physical material.A Laguerre tessellation is then generated around the cells as the pore skeleton.To control the pore geometry, a Gaussian random field is generated and then scaled to match the estimated pore-size distribution.Finally, the rescaled field is intersected with the pore skeleton to determine the pore structure.This intersection is thresholded in such a way as to give the desired porosity.

Cell structure
The cell structure is built using a sphere packing method.Each sphere is generated on a domain of size (nx, ny, nz) with a uniformly distributed midpoint, while the diameter of the given sphere is determined by sampling from the cell-size distribution estimated from the physical sample.After the new sphere is generated, it is checked against all existing spheres to ensure they do not overlap; where an overlap is detected, the midpoint is re-sampled and the overlap checked once again.This process continues until the sphere is placed in the domain without any overlap or a maximum number of locations tested is reached; for our model, we consider 1000 locations, after which the sphere is discarded.Each sphere is given a label i, which will be used in estimating the tessellation edges.
The maximum number of spheres to generate should be reflective of both the size of the domain and the cell-size distribution; if we increase the domain size the number of spheres should increase, and if the distribution is weighted more towards larger cells, the number should decrease.For our approach, we specify an average cell diameter d and radius r = d/2, and then calculate the maximum number of spheres as: rounded to the nearest integer.In generating the cell structures for our coating samples, we will consider four choices for the average cell diameter d: the mean, the median, the 75th quartile, and half the maximum diameter.An example sphere packing is shown in Fig. 2.

Pore skeleton
The skeleton for the pore structure is defined by generating a tessellation around the cell structure.The tessellation divides the threedimensional domain into a set of regions, each surrounding a particular sphere; the individual regions of the tessellation reflect the cells of the structure, and the edges of the tessellation between cells will constitute the pore skeleton.Tessellations are a common approach for the modelling of grain-based microstructures [18,19], and for our model we use a Laguerre tessellation, which has previously been used to generate virtual structures of other materials [20,10].Unlike a standard Voronoi tessellation, where the tessellation is taken just with respect to the distance to the centre of each sphere, the Laguerre tessellation uses a weighted distance that also takes into account the radius of each sphere.Using this approach for our model therefore allows us to reflect the cell distribution in our tessellation.
To generate the tessellation, we first consider a discrete lattice grid with resolution (nx, ny, nz) and a corresponding matrix M 1 that represents the tessellation regions.For every grid point x i , we calculate a weighted distance to each sphere; for a given sphere m with centre C m and radius R m , the weighted distance ω i,m is   The value of the corresponding element of M 1 is then the label of the sphere with the smallest weight.
To determine the edges of the tessellation, we consider a second binary matrix M 2 .For each element of M 2 , we compare the value of the corresponding element of M 1 to each of its neighbouring elements; if the element values are not all identical (i.e.there are at least two spheres to which the points are associated) then the element is taken as a tessellation edge and the corresponding element of M 2 given a value of 1.
Elements where these values are all identical are given a value of 0. The result of this process applied to the example cell structure is shown in Fig. 3.

Pore geometry
The geometry of the pore structure is determined using a combination of an appropriately rescaled Gaussian random field (GRF) and a distance transform to the pore skeleton.The size of the pore geometry is controlled in such a way as to match the porosity of the physical structure.We begin by generating an initial three-dimensional Gaussian random field on the domain following the approach detailed in [21], and which we give a brief summary of here.To generate a GRF G(x),x ∈ R 3 , with mean zero and covariance function Ψ(x,y), we use the fact that the covariance function can be written as where γ(p) is the spectral density of the GRF and 〈⋅, ⋅〉 the inner product.
Our field is generated with resolution N = (nx, ny, nz) and we specify a white noise parameter δ, which controls the dispersiveness of the field; increasing δ causes the points of the field to become more dispersed, whilst a smaller value causes the points to accumulate into just a few regions.To generate a structure with length scale parameter L = (L x ,L y , L z ) = δ/N, let FFT and FFT − 1 denote, respectively, the forward and inverse three dimensional Fast Fourier Transforms.We first generate an array W whose elements are independent and normally distributed with mean zero and standard deviation δ − 3 , and then compute FFT(W).We define the Fourier space grid as p = (p 1 , p 2 , p 3 ), where and similar for p 2 and p 3 .Following the approach in [21,22], we specify the spectral density on the grid as with smoothness parameters (k, n, l) = (1, 2, 2), and compute We then obtain the GRF as G = FFT − 1 (U).
To obtain a random field for which the values match the desired pore size distribution, we introduce a rescaling in the following manner; first, we fit a normal distribution to the values of the GRF using the inbuilt Matlab function fitdist, giving us the estimated parameters μ GRF and σ GRF .We then calculate the rescaled random field (RF) as where μ pore and σ pore are estimated from the physical pore structure.We generate the final pore structure by combining the rescaled random field with a distance transform D of the pore skeleton, which is computed using the Matlab function bwdist applied to the tessellation grid matrix M 2 .As the distance transform will give a value of 0 along the tessellation, we apply a small adjustment dx that will then control how strongly the pore structure adheres to the tessellation; the bigger the adjustment, the less it will adhere to the pore skeleton.We calculate the distance matrix as which is then thresholded at a specific percentile to give the desired porosity.

Structure characterisation
The stochastic model we use to generate our virtual coating structures requires estimation of the distribution of cell and pore radius of the physical material; for our model, we assume the material is statistically isotropic and hence the radii are direction-invariant.To estimate these parameters, the segmented coating structures are analysed using the software MIST, a program for the visualisation and characterisation of three-dimensional geometries [12].
The characterisation method fits to each point segmented as solid or pore the maximum diameter sphere that fits completely with the given structure; for our analysis, we assume an open boundary condition along the edges of the domain.To provide the appropriate details for the application of the characterisation method, a visualisation of the diameter values for sample A are shown in Fig. 4, and the values for the distributions can be seen below in Table 3.
We fit a lognormal distribution, with parameters μ and σ and probability density function to our size estimates, and there are two approaches we consider to determine the parameters of the distribution; for the first, we fit a lognormal distribution to the histogram output from our analysis in MIST using the maximum likelihood method.This approach ensures the best fit, but requires histogram data to be provided.For the second, the parameters are estimated using the mean (λ) and mode (ν) diameter values from the output data.Equating these values as they correspond to the distribution, we should therefore have that: which we can rearrange to give μ = ln(λ 2 ) + ln(ν) 3 The choice of mean and mode diameter value for performing the fitting is to ensure we capture the peak at the latter, but the values we get are not necessarily the best fit to the data.The calculated values of μ and σ for sample A using each approach are shown in Table 2.
We can see that for our considered sample, the two approaches give relatively similar values.Since in our case we have the histogram data for each sample, we will use the maximum likelihood method; the fitted lognormal curves for sample A using this approach are shown in red in Fig. 4.However, in instances where no such data is available, such as using the model to generate a theoretical virtual coating, using the desired mean and mode diameter is a reasonable alternative.

Analysis of material samples
The characterisation method detailed in Section 2.4 was applied to the segmented image data sets for each of the three coating samples.Each three-dimensional structure was inspected visually before the analysis to identify any potential regions on the boundaries of the structures that might adversely affect the results, due to either issues with the imaging or the structure itself; for each sample, no such areas were found, and thus the analysis was run on the full structure in each instance.The results of the analysis are shown in Table 3.
As noted, we use maximum likelihood estimation to determine the parameters for the fitted lognormal distribution for each structure of each sample.The calculated values for μ and σ are shown in Table 4.
In specifying the number of spheres to generate for the cell structure, we noted in Section 2.3 that we consider four choices for the average cell diameter: the mean, the median, the 75th quartile value, and half the maximum diameter.Each choice was tested by generating a sphere packing with the given number of spheres and comparing the results to the desired cell-size distribution.From our testing, the 75th quartile and half maximum were found to give the best comparison; since the 75th quartile value is less sensitive to outlier maximum values in the distribution, we use this as our choice for specifying the maximum number of spheres.The number of spheres to generate for each sample using this choice is 492, 2557, and 1459, respectively; the comparison for this  choice for sample A to the desired distribution is shown in Fig. 5.
Finally, for the parameters δ and dx used in specifying the pore geometry, values of δ = 0.015, 0.02, 0.025, 0.03 and dx = 5, 10, 15, 20 were tested for each sample and compared visually to the physical structures.Values of δ = 0.025 and dx = 10 were found to be the best choices for all three samples; the biggest differences were seen when choosing the value at either extremes, and thus this combination is considered a good compromise.The generated structures using these parameter inputs are shown in Fig. 6, and all three compare well visually to the real structures shown in Fig. 1.

Comparison of structures
To demonstrate the accuracy of our stochastic model, we compare the physical and virtual structures using the following geometric measures: geodesic tortuosity, cell-size distribution, and pore-size distribution.Since our focus is on fluid flow through the material, we also run mass transport simulations on each structure and compare the permeability.These particular measures were chosen to give a good representation of the pore network of each structure and understand how it compares to the physical material.To give a more robust analysis and evaluate whether any statistical effects are included, we generate five virtual structures for each sample and then take the average.

Geodesic tortuosity
The geodesic tortuosity measures the connectivity of the coating material in a given direction.For a point p in the pore space, we define the geodesic tortuosity as where the geodesic path G(p) is the shortest path that lies completely in the pore space connecting the inlet and outlet of the given direction to p, and H the height of the structure in that direction.The measure has been shown to be a good predictor of diffusive transport rates through a structure [13].
The geodesic tortuosity through the real and virtual structures for each sample were calculated using MIST, and the results of the analysis are shown in Table 5.The third row, labelled effective porosity, refers to the largest volume of connected pore structure, and provides a good indication of how the dispersion of the pores matches between the real and virtual structures.
For all three samples, the connected porosity measure compared very well between the real and virtual structures, and this was expected since the pore structure is thresholded at the desired porosity.For the mean geodesic tortuosity, the comparison for sample A was very good, but less so for the other two samples; the relative error to the real structures was 4.66%, 17.5%, and 18.8%, respectively.These results are reasonable given that the connectivity of the pore structure was not specifically controlled in the model; including estimation of the pore network and its   connectivity within the model itself is considered as a topic for future research.

Cell and pore-size distribution
In generating the virtual structures, we estimated the cell and poresize distributions of the real materials to use as input parameters for the model.To ensure that the resulting generated structures do indeed match the properties of the real structures, we re-run our characterisation analysis on the virtual structures.The averaged distribution values of the virtual structures are given in Table 6, together with the relative errors to the real structures in Table 3.
The distribution values for the cell structures are reasonably close to those of the real material, with the biggest differences being in the 75th quartile values.The results for sample C had the biggest errors, underestimating each measure except for the maximum value, which also had a large error compared to the other two samples.The differences suggest improvements can be made in ensuring that the intended cell-size distribution is better maintained as subsequent steps in the model are implemented.We note, however, that this difference in maximum value for sample C was due to a single large cell in one of the virtual structures; using instead the 97.5th quartile value, the average falls to 169.39, in line with that of the real structure.
For the pore structure, the distribution values for the virtual structure show bigger differences compared to the real materials.The median values compare reasonably well, but the mode is notably larger for the virtual structures, and the maximum value is also significantly lower, particularly for sample C. Since the porosity for the real and virtual structures are the same, these results suggest the volume of the pore structure was more spread out for the virtual structure compared to the real material, and the model failed to capture larger areas of pore volume.As we used the pore-size distribution of the real material as input for the model, these results suggest that improvements could be made in ensuring the distribution is better maintained when we intersect the Gaussian random field with the tessellation.The distance transform could be estimated more directly from the sample itself, such as by generating a simplified network through the pore structure and measuring the variance of the pore geometry to this.These improvements are considered for future research.

Mass transport simulations
Our final step is to run fluid flow simulations through both the real and virtual structures to compare the fluid permeability.These simulations were run using Gesualdo, an in-house software that computes mass transport properties of porous materials using the lattice Boltzmann method [23,14].A pressure difference is applied to the structure in the specified direction, and the average steady-state velocity is then computed.From this, the permeability can be calculated using Darcy's law.In each case, the simulations were run until convergence was reached.
Comparison of the simulation results for each structure is shown in Table 7 and we can see that the estimates for the permeability match well for samples A and C, but less so for sample B, where the virtual structure overestimates the permeability of the real structure.These comparisons do not appear to correlate with the relative errors we saw in Table 6, and suggest that even if we over or underestimate the cell and pore-size distributions, our model still captures the fluid permeability of the samples to a reasonable accuracy.

Conclusion
A tessellation-based stochastic model for generating 3D coating structures was presented and fitted to FIB-SEM tomography data of three samples of coating material.Parameters for the model were estimated from the physical structures using existing characterisation methods.The geodesic tortuosity for the real and virtual coating structures for each sample was calculated, and comparison showed that the virtual structures were a good approximation.The characterisation analysis was also applied to the virtual structures, and though the distribution values for the cell structure were comparable, the values for the pore structure were less so and suggest improvements can be made to the model to estimate this better.Fluid flow simulations were also run on each structure, and the permeability of the virtual structures was found to be a good approximation for two of the physical samples, with an overestimation on the remaining sample.Further research will consider how the model, and in particular estimation of certain parameters, can be improved.
The presented stochastic model allows us to not only generate 3D coating structures that approximate real materials, but also to generate new structures and understand how changing certain parameters affects the properties of the structure.The construction of the model was achieved with this application in mind and in the absence of real materials to estimate the input parameters from.A topic for future research would therefore be to explore this application and test different ranges for the input parameters and analyse the generated structures.This analysis can then be used to inform the design of new materials in the future.

Fig. 1 .
Fig. 1.(a) Backscatter electron images of the first slice of each coating sample; (b) classification of each pixel in the segmentation as solid (white) or pore (black); (c) three-dimensional reconstruction of each sample, with cell particles in red and the pore network in blue.Each row in the figure corresponds to sample A, B, and C, respectively, in Table1, which also includes the data for each sample.

Fig. 2 .
Fig. 2. Example of a sphere packing used to generate the cell structure.

Fig. 3 .
Fig. 3. Discrete tessellation grid for the cell structure generated by the example sphere packing.
P.Townsend et al.

Fig. 4 .
Fig. 4. Visualisation of the maximum sphere diameter at each point for the (a) cell; and (b) pore structure for sample A, together with the resulting histogram of values.Darker red areas in the visualisation correspond to larger diameter values.The fitted lognormal distributions are shown in red, appropriately rescaled to match the value at the mode.

Fig. 5 .
Fig. 5. Comparison between the desired distribution and the sphere packing generated using the 75th quartile.

Fig. 6 .
Fig. 6.Virtual coating structures generated using parameters estimated from samples A, B, and C, respectively.The cell particles are shown in red and the pore network in blue.

Table 1
Structure properties for each sample.

Table 2
Estimated parameter values for each structure for sample A.

Table 3
Characterisation results for each sample; unless otherwise stated, each value is given in voxels.

Table 4
Estimated parameter values for each sample.

Table 5
Comparison of the geodesic tortuosity for each sample.