A three-dimensional statistical model for imaged microstructures of porous polymer films

TL;DR

This paper introduces a 3D Gaussian random field model for porous polymer microstructures, enabling efficient simulation and fitting to microscopy images, and investigates how microstructure affects drug release.

Contribution

It develops a novel stochastic PDE-based model for 3D porous microstructures and provides an efficient MCMC fitting method to analyze microscopy data.

Findings

Model fits stationary microstructure regions well

Goodness-of-fit measures based on diffusion and geometry are effective

Model helps understand microstructure's impact on drug release

Abstract

A thresholded Gaussian random field model is developed for the microstructure of porous materials. Defining the random field as a solution to stochastic partial differential equation allows for flexible modelling of non-stationarities in the material and facilitates computationally efficient methods for simulation and model fitting. A Markov Chain Monte Carlo algorithm is developed and used to fit the model to three-dimensional confocal laser scanning microscopy images. The methods are applied to study a porous ethylcellulose/hydroxypropylcellulose polymer blend that is used as a coating to control drug release from pharmaceutical tablets. The aim is to investigate how mass transport through the material depends on the microstructure. We derive a number of goodness-of-fit measures based on numerically calculated diffusion through the material. These are used in combination with measures that characterize the geometry of the pore structure to assess model fit. The model is found to fit stationary parts of the material well.

Figures

Tables

Table 1. Table 1 : Differences in performance for simulating from 𝐰 ∣ 𝐬 , 𝜸 , u , 𝐲 conditional 𝐰 𝐬 𝜸 𝑢 𝐲 \mathbf{w}\mid\mathbf{s},\boldsymbol{\gamma},u,\mathbf{y} by using a PCG method with different preconditioners.Here d ​ i ​ a ​ g ​ ( 𝐐 ^ ) 𝑑 𝑖 𝑎 𝑔 ^ 𝐐 diag(\hat{{\boldsymbol{\mathrm{Q}}}}) is a diagonal matrix with the diagonal elements of 𝐐 ^ ^ 𝐐 \hat{{\boldsymbol{\mathrm{Q}}}} . The results shown were computed from 100 100 100 iterations of the MCMC algorithm for dataset H ​ P ​ C ​ 30 1 𝐻 𝑃 𝐶 subscript 30 1 HPC30_{1} , for which 𝐐 z ∈ ℝ 20 × 20 subscript 𝐐 𝑧 superscript ℝ 20 20 {\boldsymbol{\mathrm{Q}}}_{z}\in\mathbb{R}^{20\times 20} and 𝐐 𝐬 ∈ ℝ 6305 × 6305 subscript 𝐐 𝐬 superscript ℝ 6305 6305 {\boldsymbol{\mathrm{Q}}}_{{\boldsymbol{\mathrm{s}}}}\in\mathbb{R}^{6305\times 6305} , where one sample was generated in each iteration. When possible, the performance when using (5) to compute 𝐏 − 1 ​ 𝐯 = 𝐑 𝐏 − 1 ​ ( 𝐑 𝐏 − 1 ) ⊤ ​ 𝐯 superscript 𝐏 1 𝐯 superscript subscript 𝐑 𝐏 1 superscript superscript subscript 𝐑 𝐏 1 top 𝐯 {\boldsymbol{\mathrm{P}}}^{-1}{\boldsymbol{\mathrm{v}}}={\boldsymbol{\mathrm{R}}}_{{\boldsymbol{\mathrm{P}}}}^{-1}({\boldsymbol{\mathrm{R}}}_{{\boldsymbol{\mathrm{P}}}}^{-1})^{\top}{\boldsymbol{\mathrm{v}}} was compared with the performance when the full Cholesky factor 𝐑 𝐏 subscript 𝐑 𝐏 {\boldsymbol{\mathrm{R}}}_{{\boldsymbol{\mathrm{P}}}} was evaluated. The PCG CPU-time is the mean time that was spent in the PCG-routine simulating one vector from 𝐰 ∣ 𝐬 , 𝜸 , u , 𝐲 conditional 𝐰 𝐬 𝜸 𝑢 𝐲 \mathbf{w}\mid\mathbf{s},\boldsymbol{\gamma},u,\mathbf{y} ; the total CPU-time is the mean total time it took to simulate one vector, which includes the time spent in the PCG-routine, as well as the time it took to calculate the preconditioner and to generate 𝝃 𝝃 \boldsymbol{\xi} . The PCG-iterations are the mean number of iterations of the PCG-routine. The proportion of non-zero elements of the preconditioner is also indicated. A approximate minumum degree reordering (Amestoy, Davis, and Duff, 1996 ) was performed when it led to a reduction in computation time.

		CPU-time
Preconditioner $𝐏={𝐑}_{𝐏}^{\top }{𝐑}_{𝐏}$	Full Cholesky factor evaluated	PCG	Total	PCG iterations	Proportion non-zero elements
${𝐑}_{𝐏}=ichol\left({𝐐}_z\right)\otimes ichol\left({𝐐}_{𝐬}\right)$	(yes) no	(7.3) 0.9	(8.4) 1.0	(42) 42	0.003
${𝐑}_{𝐏}={𝐑}_z\otimes {𝐑}_{𝐬}$	(yes) no	(32.8) 2.9	(35.4) 2.9	(70) 70	0.004
$𝐏=diag\left(\widehat{𝐐}\right)$	-	9.3	9.3	2555	$8\cdot {10}^{-6}$
no preconditioner	-	19.4	19.4	5967	-
${𝐑}_{𝐏}=ichol\left({𝐐}_z\otimes {𝐐}_{𝐬}\right)$	yes	27.3	30.2	135	0.007
${𝐑}_{𝐏}=ichol\left(\widehat{𝐐}\right)$	yes	27.8	30.7	138	0.007

Table 2. Table 2 : Voxel sizes and number of voxels for all the samples.

Sample	Voxel size ($\mu$m3)	Number of voxels
$HPC{30}_1$	$0.23\times 0.23\times 0.20$	$74\times 74\times 20$
$HPC{30}_2$	$0.23\times 0.23\times 0.20$	$74\times 74\times 20$
$HPC{30}_3$	$0.15\times 0.15\times 0.20$	$101\times 76\times 25$
$HPC{30}_4$	$0.15\times 0.15\times 0.20$	$101\times 76\times 26$
$HPC{40}_1$	$0.15\times 0.15\times 0.10$	$81\times 81\times 36$
$HPC{40}_2$	$0.15\times 0.15\times 0.10$	$81\times 81\times 31$
$HPC{40}_3$	$0.15\times 0.15\times 0.10$	$91\times 91\times 26$
$HPC{40}_4$	$0.15\times 0.15\times 0.10$	$91\times 91\times 21$

Table 3. Table 3 : The oscillation parameters and range parameters for the model fitted to each sample. The volume fractions and surface areas for the pore structures obtained from the CLSM samples, as well as the mean values of 500 500 500 stochastic simulations from each fitted pore model with standard deviations in parentheses.

	Parameter estimates		Volume fraction		Surface area
Sample	(${\theta }_{𝐬}$, ${\theta }_z$)	(${\kappa }_{𝐬}^2$, ${\kappa }_z^2$)	Sample	Fitted model	Sample	Fitted model
$HPC{30}_1$	$(0.86,0.56)$	$(254,61)$	0.697	0.698 (0.004)	0.161	0.173 (0.004)
$HPC{30}_2$	$(0.84,0.51)$	$(219,78)$	0.694	0.694 (0.004)	0.150	0.170 (0.004)
$HPC{30}_3$	$(0.86,0.58)$	$(369,87)$	0.691	0.692 (0.004)	0.148	0.159 (0.003)
$HPC{30}_4$	$(0.88,0.56)$	$(801,112)$	0.684	0.685 (0.006)	0.210	0.217 (0.004)
$HPC{40}_1$	$(0.89,0.56)$	$(1088,271)$	0.610	0.609 (0.007)	0.254	0.271 (0.003)
$HPC{40}_2$	$(0.90,0.57)$	$(1077,272)$	0.615	0.621 (0.007)	0.248	0.269 (0.003)
$HPC{40}_3$	$(0.90,0.54)$	$(1532,319)$	0.618	0.616 (0.005)	0.291	0.306 (0.003)
$HPC{40}_4$	$(0.87,0.50)$	$(1825,384)$	0.619	0.606 (0.006)	0.345	0.358 (0.003)

Equations

$$y_i=\{\begin{array}{ll}1,& \text{ if X(si,zi)+εi≥u,}\\ 0,& \text{ otherwise,}\end{array}i=1,\dots ,m.$$

$$\tilde{\mathbf{y}}=F(\mathbf{y}),$$

$${\text{Cor}}_z(z_1,z_2)$$

$${\text{Cor}}_{\mathrm{s}}({\mathrm{s}}_1,{\mathrm{s}}_2)$$

$$\left({\kappa }_z^2\exp (\text{i}\pi {\theta }_z)-\frac{{\partial }^2}{\partial z^2}\right)\left({\kappa }_s^2\exp (\text{i}\pi {\theta }_s)-\Delta \right)\{\tau \left[X(\mathrm{s},z)+\text{i}{X}_2(\mathrm{s},z)\right]\}=W_1(\mathrm{s},z)+\text{i}{W}_2(\mathrm{s},z),$$

$${\mathrm{Q}}_{\star }({\theta }_{\star },{\kappa }_{\star })={\kappa }_{\star }^4{\mathrm{C}}_{\star }+2{\kappa }_{\star }^2\cos (\pi {\theta }_{\star }){\mathrm{G}}_{\star }+{\mathrm{G}}_{\star }{\mathrm{C}}_{\star }^{-1}{\mathrm{G}}_{\star },$$

$$\begin{array}{rl}{y}_i|\mathrm{w}& \sim Be(\Phi ({\mathrm{A}}_{i,\bullet }\mathrm{w}-u)),i=1,\dots ,m,\\ \mathrm{w}& \sim \mathsf{N}(0,\mathrm{Q}(\mathbf{\gamma }{)}^{-1}),\\ \mathbf{\gamma }& \sim \pi (\mathbf{\gamma }).\end{array}$$

$$\mathrm{C}\mathrm{v}=vec(({\mathrm{C}}_1({\mathrm{C}}_2\mathrm{V}{)}^T{)}^T),$$

$${\mathrm{C}}^{-1}\mathrm{v}=vec(({\mathrm{C}}_1^{-1}({\mathrm{C}}_2^{-1}\mathrm{V}{)}^T{)}^T),$$

$$E_{u_J,\alpha }^{+}=\underset{D\subset {\mathbb{R}}^3}{\mathrm{argmax}}\{|D|:P(D\subseteq D_{u_J}^{+})\ge 1-\alpha \}.$$

$$E_{u_J,\alpha }^{-}=\underset{D\subset {\mathbb{R}}^3}{\mathrm{argmax}}\{|D|:P(D\subseteq D_{u_J}^{-})\ge 1-\alpha \},$$

$${\mathrm{Q}}_z=\left[\begin{array}{cc}{\mathrm{Q}}_{z,EI}& {\mathrm{Q}}_{z,EE}\\ {\mathrm{Q}}_{z,II}& {\mathrm{Q}}_{z,IE}\end{array}\right],$$

$$\pi (\mathbf{s}\mid \mathbf{w},\mathbf{\gamma },u^{\prime },\mathbf{y})\propto {\phi }_{(\mathrm{A}\mathbf{w},{\sigma }^2\mathrm{I})}(\mathbf{s})\prod _{i:y_i=0}{1}_{(-\infty ,u^{\prime })}(s_i)\prod _{i:y_i=1}{1}_{[u^{\prime },\infty )}(s_i),$$

$$\pi (\tau \mid \mathbf{w},\mathbf{s},{\theta }_{\mathrm{s}}^{\prime },{\theta }_z^{\prime },{\kappa }_{\mathrm{s}}^{\prime },{\kappa }_z^{\prime },u,\mathbf{y})\propto \pi (\tau ){\phi }_{(0,{\mathrm{Q}}^{-1}(\mathbf{\gamma }))}(\mathbf{w})$$

$$\Xi \circ B=\{x\in \Xi :\exists z\text{ such that }x\in z+B\text{ and }z+B\subseteq \Xi \}.$$

$$S(\lambda ;x,\Xi ,B)=P(h(x,\Xi ,B)\ge \lambda \mid x\in \Xi ),$$

$$\hat{S}(\lambda )=\frac{\hat{P}(x\in \Xi \circ \lambda B)}{\hat{P}(x\in \Xi )}=\frac{v_d((\Xi \cap W)\circ \lambda B)/v_d(W)}{v_d(\Xi \cap W)/v_d(W)},$$

Full text

A three-dimensional statistical model for imaged microstructures of porous polymer films

Sandra Barman and David Bolin111Correspondence to: David Bolin, Department of Mathematical Sciences, University of Gothenburg, 41296 Gothenburg, Sweden. Tel: +46 31 772 53 75; email: david.bolin@chalmers.se

Department of Mathematical Sciences

Chalmers University of Technology and the University of Gothenburg

Abstract: A thresholded Gaussian random field model is developed for the microstructure of porous materials. Defining the random field as a solution to stochastic partial differential equation allows for flexible modelling of non-stationarities in the material and facilitates computationally efficient methods for simulation and model fitting. A Markov Chain Monte Carlo algorithm is developed and used to fit the model to three-dimensional confocal laser scanning microscopy images. The methods are applied to study a porous ethylcellulose/hydroxypropylcellulose polymer blend that is used as a coating to control drug release from pharmaceutical tablets. The aim is to investigate how mass transport through the material depends on the microstructure. We derive a number of goodness-of-fit measures based on numerically calculated diffusion through the material. These are used in combination with measures that characterize the geometry of the pore structure to assess model fit. The model is found to fit stationary parts of the material well.

Keywords: Porous media, Gaussian field, Gaussian Markov random field, Markov Chain Monte Carlo, model validation.

Introduction

Studies of porous media have applications in many areas, ranging from geophysics, energy, electrical and chemical engineering, composite material design, and biomedical and pharmaceutical science (Vafai, 2015; Torquato, 2002), to the oil industry (Blunt, Bijeljic, Dong, Gharbi, Iglauer, Mostaghimi, and Paluszny, 2013), and modelling of transport in polymer-electrolyte fuel cells (Weber et al., 2014).

The macroscopic properties—e.g. electrical conductivity, heat or mass transfer—depend not only on the properties of the material, but also on the microscopic geometry of the pore structure. It is therefore important to characterize the microstructure, and determine how it influences the macroscopic properties. We study the microstructure of a porous ethylcellulose/hydroxypropylcellulose (EC/HPC) polymer blend, which is used as a coating to control the drug release from pharmaceutical tablets. To create coatings with desired transport properties it is important to understand how the properties of the microstructure depend on the different manufacturing parameters, and how those properties influence mass transport—in this case diffusion—of the drug through the coating.

Statistical properties of the microstructure of porous media can be obtained from microscopy data, which then can be related to the macroscopic properties. However, producing the porous material and doing the microscopy imaging can be time-consuming and costly. Formulating a stochastic model for the pore structure, which can be fitted to microscopy images of the material, is therefore useful. It is easier and cheaper to control and change different microscopic properties of the stochastic simulations from the model than it would be to produce samples of the material with these properties. By determining the macroscopic properties of each stochastic simulation, the simulations can be used to understand how the microscopic properties influence the macroscopic properties.

Thresholded Gaussian random fields were among the first stochastic models used to reconstruct heterogeneous materials (Torquato, 2002, p. 295), and such models have been used to characterize the microstructure and its relation to macroscopic properties of porous media in e.g. Adler, Jacquin, and Quiblier (1990); Roberts and Teubner (1995); Mukherjee and Wang (2007). Other examples of stochastic models and methods used to make inferences about macroscopic properties of porous media are random set models (Hermann and Elsner, 2014), network models (Blunt, Jackson, Piri, and Valvatne, 2002; Gaiselmann, Neumann, Schmidt, Pecho, Hocker, and Holzer, 2014), process-based models (Malek, Eikerling, Wang, Navessin, and Liu, 2014), and simulated annealing-based reconstructions (Yeong and Torquato, 1998; Kim and Pitsch, 2009). These models are usually validated by comparing the statistical properties of the model with those of the imaged pore structures. Also macroscopic properties computed for the stochastic simulations from the model are often compared with those computed for the pore structure obtained from microscopy images, as well as with experimental results. These macroscopic properties can be estimated in different ways: by computing analytical bounds or by solving the governing equations approximately on simplified network models or in the pore structure (see e.g. Torquato (2002, Chapter 21); Mukherjee, Kang, and Wang (2011); Weber et al. (2014)).

In this work, we develop a new class of thresholded Gaussian random field models for the microstructure of porous materials. Instead of specifying the Gaussian random field through its mean and covariance function, we specify the field as a solution to a specific stochastic partial differential equation (SPDE). The reason for this is twofold: Firstly, the model can easily be extended to handle non-stationarities in the material by allowing for spatially varying parameters in the SPDE. Secondly, a computationally efficient representation of the model can be constructed by discretising the SPDE using a finite element method. Using this representation, we develop a Markov Chain Monte Carlo (MCMC) algorithm that can be used to fit the model to large three-dimensional confocal laser scanning microscopy images.

We also construct a number of goodness-of-fit measures based on numerically calculated diffusion through the material. These can be combined with measures that characterize the geometry of the pore structure to test the goodness-of-fit when the model is used in practice. We use the methods to study the microstructure of EC/HPC polymer blends with two different weight ratios. The polymer blend is used as coating which controls drug release from pharmaceutical tablets and the aim is to investigate how mass transport properties of the material depend on the microstructure. The model is fitted to stationary but anisotropic parts of confocal laser scanning microscopy (CLSM) images of the EC/HPC films and is found to fit the data well.

The article is structured as follows. In the next section, the data, model, estimation procedure, and validation methods are presented. After this, the results for the EC/HPC data are presented, and the article concludes with a discussion of the results and future work. The supplementary material for the article contains three appendices that gives additional details and results.

Materials and methods

CLSM images of porous EC/HPC polymer films

EC/HPC polymer blends are used as coatings to control the release of drug from tablets. The HPC is water soluble while the EC is not, and so the HPC will be leached out when the tablet is immersed in water, creating a porous structure through which drug can be transported. The microstructure of the EC/HPC polymer blend is determined by the polymers’ weight ratios and molecular weights, and by processing parameters such as the temperature and spraying rate used in the manufacturing of the coating (Marucci, Hjärtstam, Ragnarsson, Iselau, and Axelsson, 2009, 2013; Andersson, Hjärtsam, Stadig, von Corswant, and Larsson, 2013). We analyse CLSM images of two EC/HPC free films—i.e. films that have been sprayed onto a rotating drum—one with $30\%$ and one with $40\%$ HPC weight ratio. More details about the preparation and imaging of these films can be found in Häbel, Rajala, Boissier, Marucci, Schladitz, Redenbach, and Särkkä (2017). The EC/HPC material has previously been studied experimentally in Marucci, Hjärtstam, Ragnarsson, Iselau, and Axelsson (2009, 2013); Andersson, Hjärtsam, Stadig, von Corswant, and Larsson (2013) and characterized statistically in Häbel, Andersson, Olsson, Olsson, Larsson, and Särkkä (2016, 2017). The mass transport properties of the material have been investigated using numerically calculated diffusion in Gebäck, Marucci, Boissier, Arnehed, and Heintz (2015).

In the manufacturing of the films, the EC and HPC phase separate into regions rich in EC and regions rich in HPC. The phase separation continues until the solution has dried. Layers are sprayed onto the film one after the other. The drying of the film is slightly different at the top of the film compared to at the bottom, causing the average size of the pores to be smaller at the top than at the bottom. This is expected to have an effect on the pore shapes, and we observed that the pore shapes are slightly elongated in the direction perpendicular to the layers (depth), i.e. the pore shapes are anisotropic.

To fit our model, we selected parts of the CLSM image stacks which have approximately homogeneous pore sizes: four image samples $HPC30_{1},\dots,HPC30_{4}$ from the HPC $30\%$ film, and four image samples $HPC40_{1},\dots,HPC40_{4}$ from the HPC $40\%$ film. The eight sets of CLSM images were chosen to be representative of the variation in pore size in the middle part of the two films.

The intensity of the signal decreases with depth. To account for this, we first thresholded (binarized) the CLSM images with a threshold that was a function of depth. See Figure 1 for a comparison of binarized and original CLSM images. The threshold was determined by fitting a quadratic function to the $70\%$ and $60\%$ intensity value quantiles for the $30\%$ and $40\%$ HPC samples respectively, calculated for each layer individually. The resulting threshold decreases slowly with depth.

The pore structure model

In this section, we define the model for the pore structure. To simplify indexing, we define ${\boldsymbol{\mathrm{s}}}=(x,y)$ so that a location $(x,y,z)$ in the CLSM image is written as $({\boldsymbol{\mathrm{s}}},z)$, where $z$ corresponds to the depth of the film. The model for the noisy binarized CLSM image $\mathbf{y}$ is obtained by taking a Gaussian field $X({\boldsymbol{\mathrm{s}}},z)$ with additive noise and thresholding it, as

[TABLE]

Here $m$ is the number of voxels in the sample, $u$ is the threshold, and $\operatorname{\varepsilon}_{i}$ are independent standard Gaussian variables that capture measurement noise. The model for the pore structure $\tilde{\mathbf{y}}$ is obtained by smoothing $\mathbf{y}$ with a filter $F$,

[TABLE]

where $\tilde{y}_{i}=1$ means that there is a pore at voxel $i$, and $\tilde{y}_{i}=0$ that there is no pore. The filter that is used is a simple mean value filter, which takes the unweighted average in a 3-by-3-by-3 or a 5-by-5-by-5 neighbourhood depending on the pore sizes of the sample, combined with a global re-thresholding to obtain the correct volume fraction. See Figure 2 for an example of a sample from $X({\boldsymbol{\mathrm{s}}},z)$ and the corresponding model structures, obtained using the model for $X({\boldsymbol{\mathrm{s}}},z)$ that is defined in the next section.

Separable oscillating Gaussian Matérn fields

To capture the different shapes of pores in the depth of the film (corresponding to the z-direction) compared to the within the layers (the ${\boldsymbol{\mathrm{s}}}$-planes), the stochastic model should be anisotropic and allow for a different behaviour in the $z$-direction. We do this by using a Gaussian random field $X({\boldsymbol{\mathrm{s}}},z)$ with a separable covariance function given as $\text{Cov}(X({\boldsymbol{\mathrm{s}}}_{1},z_{1}),X({\boldsymbol{\mathrm{s}}}_{2},z_{2}))=\sigma^{2}\text{Cor}_{{\boldsymbol{\mathrm{s}}}}({\boldsymbol{\mathrm{s}}}_{1},{\boldsymbol{\mathrm{s}}}_{2})\text{Cor}_{z}(z_{1},z_{2})$, where $\sigma^{2}$ is the variance of the process. The anisotropy in the $z$-direction can then be controlled by letting the correlation function $\text{Cor}_{z}$ be different from the correlation function $\text{Cor}_{{\boldsymbol{\mathrm{s}}}}$.

The imaged pore structures have regularly spaced pores, as can be seen from Figure 1. To obtain a Gaussian field that is regular we need to use an oscillating covariance function. Covariance functions with differing oscillating strengths should be allowed for, since covariance functions estimated from different binarized CLSM image samples are oscillating to various degrees. We use the family of oscillating Matérn correlation functions

[TABLE]

introduced by Lindgren et al. (2011). Here $K_{0}$ is the modified Bessel function of the second kind and order zero. The positive parameters $\kappa_{*}>0$ scales the correlation functions, and so control the correlation ranges and the spacing between the correlation functions’ peaks. The spacing determines the typical distance between neighbouring peaks in the Gaussian field. The parameters $\theta_{*}\in[0,1)$ control the amount of oscillation in the correlation functions. For $\theta_{*}=0$, the functions are equal to the well-known regular Matérn correlation function, with smoothness parameter $\nu=1$ for ${\boldsymbol{\mathrm{s}}}$ and $\nu=3/2$ for $z$. The oscillation of the correlation functions is more pronounced the closer $\theta_{*}$ is to one, and the corresponding Gaussian field will therefore be more regular for higher values of $\theta_{*}$. Examples of the correlation functions $\text{Cor}_{z}$ and $\text{Cor}_{\boldsymbol{\mathrm{s}}}$ are shown in Figure 3. Figure 4 shows correlations $\text{Cor}_{z}\text{Cor}_{\boldsymbol{\mathrm{s}}}$ for the Gaussian field $X({\boldsymbol{\mathrm{s}}},z)$, with parameters taken from the model fitted to the dataset $HPC30_{1}$, where the oscillation parameters were $\theta_{{\boldsymbol{\mathrm{s}}}}=0.86$ and $\theta_{z}=0.56$. A stochastic simulation of this Gaussian field is shown in Figure 2(a). The correlation is symmetric in the ${\boldsymbol{\mathrm{s}}}$-plane since $\text{Cov}_{{\boldsymbol{\mathrm{s}}}}$ is isotropic, and has a different behaviour along the $z$-axis since the joint covariance is anisotropic.

Lindgren, Rue, and Lindström (2011) derived the oscillating Matérn covariance function as the covariance function of the solution to a certain SPDE. By similar arguments, one can show that $X({\boldsymbol{\mathrm{s}}},z)$ in our case can be represented as the solution to the SPDE

[TABLE]

where $\Delta=\frac{\operatorname{\partial}^{2}}{\operatorname{\partial}x^{2}}+\frac{\operatorname{\partial}^{2}}{\operatorname{\partial}y^{2}}$ is the Laplace operator, $\tau$ is a parameter that controls the variance, and $W_{1}$ and $W_{2}$ are independent Gaussian white noise fields. A stationary weak solution $(X,X_{2})$ to this SPDE has the property that the component fields $X$ and $X_{2}$ are real, independent, and have the same distribution. These fields are stationary and isotropic and we call them separable oscillating Gaussian Matérn fields.

There are two important reasons for representing $X$ using the SPDE (3). First, it allows us to define non-stationary extensions of the separable oscillating Gaussian Matérn field by allowing the parameters to be spatially varying. This will be used in future work to analyse larger samples of EC/HPC films where the stationarity assumption is problematic. Secondly, the representation makes it possible to approximate the Gaussian field using a Gaussian Markov random field (GMRF). A GMRF is a multivariate Gaussian random vector with a sparse precision matrix (inverse covariance matrix). As we will see later, the sparsity greatly reduces the computational cost of fitting the model to data.

The Gaussian Markov random field approximation

A GMRF approximation of the separable oscillating Gaussian Matérn field $X({\boldsymbol{\mathrm{s}}},z)$ is obtained by solving the SPDE (3) approximately using the finite element method, on a bounded domain $\Omega=\Omega_{{\boldsymbol{\mathrm{s}}}}\times\Omega_{z}\subset\mathbb{R}^{3}$ with Neumann boundary conditions. The idea is to approximate $X({\boldsymbol{\mathrm{s}}},z)$ by a basis expansion $X_{\text{FEM}}({\boldsymbol{\mathrm{s}}},z)=\sum_{i=1}^{n}w_{i}\varphi_{i}({\boldsymbol{\mathrm{s}}},z)$, where $\{\varphi_{i}\}$ is a suitable basis and ${\boldsymbol{\mathrm{w}}}$ are stochastic weights. Here a good choice of basis is the Kronecker basis $\{\varphi_{ij}({\boldsymbol{\mathrm{s}}},z)=\varphi_{{\boldsymbol{\mathrm{s}}},i}({\boldsymbol{\mathrm{s}}})\varphi_{z,j}(z)\}$, where $\{\varphi_{{\boldsymbol{\mathrm{s}}},i}\}$ and $\{\varphi_{z,i}\}$ are piecewise linear and continuous basis functions obtained by a triangulation of $\Omega_{{\boldsymbol{\mathrm{s}}}}$ and $\Omega_{z}$ respectively. The basis functions $\varphi_{{\boldsymbol{\mathrm{s}}},i}$ and $\varphi_{z,i}$ take the value one in node $i$ and zero in all other nodes in their respective triangulations. The distribution of the stochastic weights ${\boldsymbol{\mathrm{w}}}$ is derived using the Galerkin method (Lindgren, Rue, and Lindström, 2011). Straightforward calculations give that ${\boldsymbol{\mathrm{w}}}\sim N(0,{\boldsymbol{\mathrm{Q}}}({\boldsymbol{\mathrm{\gamma}}})^{-1})$, where ${\boldsymbol{\mathrm{\gamma}}}=(\theta_{s},\kappa_{s},\theta_{z},\kappa_{z},\tau)$ and ${\boldsymbol{\mathrm{Q}}}({\boldsymbol{\mathrm{\gamma}}})=\tau^{2}{\boldsymbol{\mathrm{Q}}}_{{\boldsymbol{\mathrm{s}}}}(\theta_{s},\kappa_{s})\otimes{\boldsymbol{\mathrm{Q}}}_{z}(\theta_{z},\kappa_{z})$. The Kronecker structure of the precision matrix ${\boldsymbol{\mathrm{Q}}}({\boldsymbol{\mathrm{\gamma}}})$ is a result of the separability of the covariance of $X({\boldsymbol{\mathrm{s}}},z)$, and of the choice of basis functions. The two matrices ${\boldsymbol{\mathrm{Q}}}_{{\boldsymbol{\mathrm{s}}}}$ and ${\boldsymbol{\mathrm{Q}}}_{z}$ are given by

[TABLE]

for $\star={\boldsymbol{\mathrm{s}}},z$, where ${\boldsymbol{\mathrm{C}}}_{\star}$ is a diagonal matrix with diagonal elements $C_{\star,ii}=\int_{\Omega}\varphi_{\star,i}(t)dt$ and ${\boldsymbol{\mathrm{G}}}_{\star}$ is a sparse matrix with elements $G_{\star,ij}=\int_{\Omega}\nabla\varphi_{\star,i}(t)\nabla\varphi_{\star,j}(t)dt$. Since the matrices ${\boldsymbol{\mathrm{C}}}_{\star}$ and ${\boldsymbol{\mathrm{G}}}_{\star}$ are sparse, so is ${\boldsymbol{\mathrm{Q}}}_{\star}$, and the weights ${\boldsymbol{\mathrm{w}}}$ is therefore a GMRF with the sparse precision matrix ${\boldsymbol{\mathrm{Q}}}({\boldsymbol{\mathrm{\gamma}}})$.

With this approach $X_{\text{FEM}}({\boldsymbol{\mathrm{s}}},z)$ is a piecewise linear approximation of $X({\boldsymbol{\mathrm{s}}},z)$. We refer to the stochastic weights ${\boldsymbol{\mathrm{w}}}$ as the separable oscillating Matérn GMRF. The value of $X_{\text{FEM}}({\boldsymbol{\mathrm{s}}},z)$ at the node of a basis function $\varphi_{ij}({\boldsymbol{\mathrm{s}}},z)$ is the corresponding weight. When the triangulation used to construct the basis functions is refined, the approximating field $X_{\text{FEM}}$ converges weakly to $X$ (Lindgren, Rue, and Lindström, 2011) .

The domain $\Omega$ is chosen to be slightly larger than the domain of the CLSM image to reduce the effect of the boundary conditions imposed on the SPDE. We use one basis function for each voxel in the sample, with the node of the basis function at the center of the voxel, as well as basis functions outside the domain of the sample. For further details, see Appendix A.

Markov Chain Monte Carlo model fitting

To fit the model to the noisy binarized CLSM image sample $\mathbf{y}$, the parameters ${\boldsymbol{\mathrm{\gamma}}}$ and the threshold $u$ have to be estimated. Evaluating the likelihood of the model is too computationally demanding, so we estimate the model using Markov Chain Monte Carlo (MCMC) in a Bayesian context instead of doing maximum likelihood estimation.

Let ${\boldsymbol{\mathrm{A}}}\in\mathbb{R}^{m\times n}$ be a matrix that maps the values of the GMRF ${\boldsymbol{\mathrm{w}}}$ to the values of the piecewise linear approximation $X_{\text{FEM}}({\boldsymbol{\mathrm{s}}},z)$ at the voxels in the CLSM image, where as before $m$ is the number of voxels and $n$ is the number of stochastic weights. Specifically, element $ij$ in ${\boldsymbol{\mathrm{A}}}$ is obtained by evaluating the $j$th basis function in the location of the $i$th voxel, so that ${\boldsymbol{\mathrm{A}}}{\boldsymbol{\mathrm{w}}}$ is a vector that contains the values for $X_{\text{FEM}}({\boldsymbol{\mathrm{s}}},z)$ at the voxels of the CLSM image.

Approximating $X({\boldsymbol{\mathrm{s}}},z)$ with $X_{\text{FEM}}({\boldsymbol{\mathrm{s}}},z)$, we have from model (1) that $y_{i}=1$ if ${\boldsymbol{\mathrm{A}}}_{i,\bullet}{\boldsymbol{\mathrm{w}}}+\epsilon_{i}\geq u$, and $y_{i}=0$ otherwise. This means that we can write the full model as

[TABLE]

Here $Be$ denotes the Bernoulli distribution, $\Phi$ denotes the distribution function of a standard normal random variable, and $\pi({\boldsymbol{\mathrm{\gamma}}})$ is a prior distribution for the parameters.

The MCMC algorithm we use to estimate the parameters is a blocked Metropolis-within-Gibbs sampler, see Appendix A for details. The algorithm is used to estimate the joint posterior density $\pi(\mathbf{w},\boldsymbol{\gamma},u\mid\mathbf{y})$ and the posterior expectation $\hat{{\boldsymbol{\mathrm{\gamma}}}}=E({\boldsymbol{\mathrm{\gamma}}}|\mathbf{y})$ is taken as a point estimate of the parameters. When simulating from the fitted model, we draw parameters from the posterior distribution of the parameters.

Sampling from the posterior distribution of ${\boldsymbol{\mathrm{w}}}$ is the most computationally demanding part of the MCMC algorithm. A method for reducing the complexity of the sampling is presented in Appendix A. The method takes advantage of the Kronecker structure of the precision matrix ${\boldsymbol{\mathrm{Q}}}({\boldsymbol{\mathrm{\gamma}}})=\tau^{2}{\boldsymbol{\mathrm{Q}}}_{{\boldsymbol{\mathrm{s}}}}\otimes{\boldsymbol{\mathrm{Q}}}_{z}$ and its Cholesky factor, and the sparsity of these matrices. The most important implication of the Kronecker structure, which we use to improve the computational efficiency of the method, is that matrix-vector operations can be computed efficiently as follows.

For a matrix ${\boldsymbol{\mathrm{C}}}={\boldsymbol{\mathrm{C}}}_{1}\otimes{\boldsymbol{\mathrm{C}}}_{2}$, where ${\boldsymbol{\mathrm{C}}}_{1}\in\mathbb{R}^{n_{1}\times n_{1}}$ and ${\boldsymbol{\mathrm{C}}}_{2}\in\mathbb{R}^{n_{2}\times n_{2}}$, we have that

[TABLE]

where $vec$ denotes the vectorization operator which maps a matrix to a vector by stacking the columns of the matrix, and ${\boldsymbol{\mathrm{V}}}$ is an $n_{2}\times n_{1}$ matrix containing the elements of ${\boldsymbol{\mathrm{v}}}$ with ${\boldsymbol{\mathrm{V}}}_{ij}={\boldsymbol{\mathrm{v}}}_{(j-1)n_{2}+i}$ (Buis and Dyksen, 1996). Using the same trick for the inverse we get

[TABLE]

using that ${\boldsymbol{\mathrm{C}}}^{-1}={\boldsymbol{\mathrm{C}}}_{1}^{-1}\otimes{\boldsymbol{\mathrm{C}}}_{2}^{-1}$.

For full matrices ${\boldsymbol{\mathrm{C}}}_{1}$ and ${\boldsymbol{\mathrm{C}}}_{2}$, computing ${\boldsymbol{\mathrm{C}}}^{-1}{\boldsymbol{\mathrm{v}}}$ without taking advantage of the Kronecker structure requires $\mathcal{O}((n_{1}n_{2})^{3})$ operations, whereas the computation $vec(({\boldsymbol{\mathrm{C}}}_{1}^{-1}({\boldsymbol{\mathrm{C}}}_{2}^{-1}{\boldsymbol{\mathrm{V}}})^{T})^{T})$ only requires $\mathcal{O}(n_{1}^{3}+n_{2}^{3})$ operations, which is a considerable reduction. In the method for posterior sampling, the matrix-vector operations (5) and (6) can be performed efficiently due to the sparsity of the Kronecker product matrices, further reducing the computational cost. Moreover, using (5) and (6), the MCMC algorithm can be implemented so that we never need to store full Kronecker products. These reductions in algorithmic complexity are crucial for making it possible to fit the model to images of realistic sizes. Further details are given in Appendix A.

Model validation

To assess the model fit we need to measure the similarities between the binarized CLSM image and the pore structure obtained by filtering the binarized CLSM image, and the stochastic simulations from the corresponding fitted models (4) and (2). We have chosen the following measures: the empirical covariance function, size distributions with respect to linear and spherical shapes, and measures related to numerically calculated diffusion. Other popular summarizing functions that could be applied as measures of goodness-of-fit are, e.g., the contact distribution functions, empty space functions and chord length distribution functions (see Chiu, Stoyan, Kendall, and Mecke, 2013, for a review). These are all related to the size distributions. Below we present the goodness-of-fit measures in more detail.

Covariance functions

As a first assessment of the model fit, we look at marginal covariance functions for both models (4) and (2). To obtain the theoretical covariance of model (4), we use the following result for a thresholded Gaussian vector: If $(X,Y)$ is a zero-mean bivariate Gaussian vector with unit variances and correlation $\rho\in(-1,1)$, then the covariance between the elements of the vector thresholded at levels $u_{X}$ and $u_{Y}$ respectively is $\int_{0}^{\rho}\varphi(u_{X},u_{Y};z)dz$, where $\varphi(\cdot,\cdot;z)$ is the density of a bivariate Gaussian vector with unit variances and correlation $z$ (Cramér and Leadbetter, 1967, pp. 26–27). The covariances of the model are obtained using that $(({\boldsymbol{\mathrm{A}}}_{i\cdot}{\boldsymbol{\mathrm{w}}}+\epsilon_{i})/\alpha_{i})_{i=1}^{n}$ is a zero-mean Gaussian vector with unit variances and correlations $\rho_{i,j}=({\boldsymbol{\mathrm{Q}}}(\boldsymbol{\gamma})^{-1})_{i,j}/(\alpha_{i}\alpha_{j})$, $i,j\in\{1,\dots,m\}$, where $\alpha_{i}^{2}=\text{Var}({\boldsymbol{\mathrm{A}}}_{i\cdot}{\boldsymbol{\mathrm{w}}}+\epsilon_{i})=({\boldsymbol{\mathrm{Q}}}(\boldsymbol{\gamma})^{-1})_{i,i}+\sigma^{2}$. Since the marginal distribution of the field evaluated in an ${\boldsymbol{\mathrm{s}}}$-plane is stationary and isotropic except for boundary effects, the variances $\alpha_{i}^{2}$ will be approximately constant, and the correlations $\rho_{i,j}$ will depend only on the distance between the voxels $i$ and $j$ for voxels lying in the same ${\boldsymbol{\mathrm{s}}}$-plane (as long as the voxels are far enough away from the boundary). Therefore we consider a single marginal covariance function $C_{{\boldsymbol{\mathrm{s}}}}$ for voxels lying in the same ${\boldsymbol{\mathrm{s}}}$-plane. Similarly, we consider a single marginal covariance function $C_{z}$ for voxels lying on the same $z$-line.

We compare the empirical marginal covariance functions $\hat{C}_{{\boldsymbol{\mathrm{s}}}}$ and $\hat{C}_{z}$ estimated from the binarized CLSM image $\mathbf{y}$, with envelopes estimated from stochastic simulations $\mathbf{y}_{\text{model}}$ from the fitted noise model (4), as well as with the theoretical covariance functions $C_{{\boldsymbol{\mathrm{s}}}}$ and $C_{z}$ of this model. We also compare the empirical marginal covariance functions estimated from the filtered CLSM pore structure $\tilde{\mathbf{y}}$, with envelopes estimated from stochastic simulations $\tilde{\mathbf{y}}_{\text{model}}$ from the fitted pore model (2). We use simultaneous $1-\alpha$-envelopes, defined as in Bolin and Lindgren (2016) so that $1-\alpha$ percent of the covariance functions estimated from stochastic simulations lie completely within the envelopes.

Size distributions

The morphological size distributions, or granulometries, were introduced in Matheron (1967) as a way of characterizing porous media and random closed sets, and is a common tool in image analysis (Soille, 2004). Size distributions provide a measure of size for sets—such as porous structures—that do not have well-defined shapes. Local sizes are measured with respect to simpler sets, so called structuring elements, such as spheres and line segments, and the set of voxels in the three-dimensional pore space of a pore structure $\tilde{\mathbf{y}}$ is seen as a realization of a random set $\Xi_{\text{pores}}=\{i:\tilde{y}_{i}=1\}$, which consists of those voxels $i$ where there is a pore.

To define the size distribution of a random set $\Xi$ with respect to a structuring element $B$, first define the local size $h(x,\Xi,B)$ of a point $x\in\Xi$ with respect $B$ as the largest value for which the rescaled structuring element $B$ can be translated so that it fits within the set and covers the point $x$. Figure 5 illustrates this local size concept, in this case for two-dimensional sets, showing local sizes with respect to three structuring elements $B$. The local sizes with respect to the circle (Panel (a)) are smaller in the middle of the ellipse compared to the local sizes with respect to the lines (Panels (b) and (c)), since it is possible to fit longer lines within the ellipse covering these pixels than it is possible to fit circles with those lines as radii.

The size distribution $S(\lambda;x,\Xi,B)=P(h(x,\Xi,B)\geq\lambda\mid x\in\Xi)$ gives the probability that a point $x\in\Xi$ has a local size with respect to $B$ which is greater than or equal to $\lambda$. We let $s(\lambda;x,\Xi,B)$) denote the corresponding size density. The size distribution and density do not depend on the point $x$ if the random set is stationary.

For the discretely indexed random set $\Xi_{\text{pores}}=\{i:\tilde{y}_{i}=1\}$, we let the distance between two neighbouring voxels be one. We estimate stationary size distributions $S(\lambda;\Xi_{\text{pores}},B)$, with lines aligned with the coordinate axes and the sphere as structuring elements. In addition, we also estimate these four size distributions for the random set $\Xi_{\text{matrix}}=\{i:\tilde{y}_{i}=0\}$, which consists of the voxels in the matrix (the solid part of the pore structure). The size distributions estimated from stochastic simulations $\tilde{\mathbf{y}}_{\text{model}}$ from the pore model (2) are compared to the size distributions estimated from the filtered CLSM pore structure $\tilde{\mathbf{y}}$, using simultaneous envelopes. For details about the estimation, see Appendix B.

Numerically calculated diffusion

The final aim is to use the developed model to analyze how mass transport—mainly diffusion—of drug through the EC/HPC coating depends on the microstructure. Because of this, we in this work use numerical calculations of diffusion to test the model fit. We compare the diffusion calculated in the filtered CLSM pore structure $\tilde{\mathbf{y}}$ with the diffusion calculated in stochastic simulations $\tilde{\mathbf{y}}_{\text{model}}$ from the fitted pore model (2).

The diffusion calculations are done using the software Gesualdo (Gebäck and Heintz, 2014), which provides a diffusive flux vector $\mathbf{J}(i)=(J_{x},J_{y},J_{z})(i)$ for each voxel $i$ in the pore space (see Figure 6 for an illustration). The diffusive flux $J_{x}(i)$ gives the amount of particles transported in steady state in the $x$-direction, per unit square area and second, $J_{y}(i)$ the amount transported in the $y$-direction, and $J_{z}(i)$ the amount transported in the $z$-direction. Since the $z$-direction is the direction of transport, the effective diffusion coefficient for the pore structure is defined as the average of $J_{z}$, taken over the whole structure including the matrix where $J_{z}$ is set to zero.

We compare both the effective diffusion coefficient and the diffusive flux. The comparison of the diffusive flux is done using 3D histograms of the diffusive flux vectors. Since a stochastic simulation $\tilde{\mathbf{y}}_{\text{model}}$ from the model is random, the corresponding flux histogram of the diffusive flux vectors $\{\mathbf{J}(i)\}_{i:\tilde{y}_{\text{model},i}=1}$ is a random histogram. To test if the random histograms for the simulated pore structures are similar to that from the CLSM pore structure, we test if the set where the model histogram values are usually high is included in the set where the CLSM histogram values are high, and if there are no high values in the CLSM histogram where the model histogram values are usually low. To do this formally we use the excursion set methodology by Bolin and Lindgren (2015).

The set where the model histogram values are usually high is characterised by a positive excursion set. This is defined as a maximal set in which, with a given probability $1-\alpha$, all values of the model histogram are above a given value $u_{J}$:

[TABLE]

Here $D_{u_{J}}^{+}=\{b:H(b)>u_{J}\}$ is the random set of 3D bins where the model histogram $H$ exceeds the value $u_{J}$. Similarly, we define the negative excursion set as a maximal set where for most model histograms all values are below the given value $u_{J}$:

[TABLE]

where $D_{u_{J}}^{-}=\{b:H(b)<u_{J}\}$. Set $E_{\text{CLSM},u_{J}}=\{b:H_{\text{CLSM}}(b)>u_{J}\}$, where $H_{\text{CLSM}}$ is the diffusive flux histogram of the CLSM pore structure. Selecting a high probability $1-\alpha$, we can then test if $E_{u_{J},\alpha}^{+}\subseteq E_{\text{CLSM},u_{J}}\subseteq(E_{u_{J},\alpha}^{-})^{c}$, where $E_{u_{J},\alpha}^{+}$ and $E_{u_{J},\alpha}^{-}$ are estimated from model simulations. If this inequality is satisfied, it implies that flux field through the observed CLSM pore structure is similar to flux fields obtained from model simulated pore structures.

Computation

The MCMC algorithm was implemented in Matlab (MATLAB, 2015) and so was the estimation of the size distributions. The triangulation of the domain was done using the R package INLA (Rue, Martino, and Chopin, 2009) and the 2D-covariance estimation was done using the R package RandomFields (Schlater, Malinowski, Menck, Oesting, and Strokorb, 2015). The excursion sets and simultaneous envelopes were estimated using the R package excursions (Bolin and Lindgren, 2016). The diffusive flux was calculated using the software Gesualdo (Gebäck and Heintz, 2014).

Results

In this section, detailed results for the model fitted to the CLSM image $HPC30_{1}$ are presented. The results for the other CLSM image samples are similar, see Appendix C. The MCMC algorithm was run for $2\cdot 10^{5}$ iterations, making sure that the chains reach stationarity. Each iteration took around $1$–$2$ seconds. About $85$–$90\%$ of each iteration was spent sampling the posterior distribution of the GMRF.

Regarding the prior densities for the parameters, we used uniform priors for the oscillation parameters, $\theta_{{\boldsymbol{\mathrm{s}}}},\theta_{z}\sim U(0,1)$, an improper uniform prior on the entire $\mathbb{R}$ for the threshold $u$, and gamma priors for the range and variance parameters, $\kappa_{{\boldsymbol{\mathrm{s}}}}^{2},\kappa_{z}^{2}\sim\Gamma(1,6\cdot 10^{-5}),\tau^{2}~{}\sim~{}\Gamma(1,5\cdot 10^{-3})$. Figure 7 shows these densities as well as the estimated posterior densities. Because of the vast amount of data, the posterior distributions of the parameters are insensitive to the priors and have quite small variances. The posterior distribution for the oscillation parameter $\theta_{z}$ has a higher variance than the oscillation parameter $\theta_{{\boldsymbol{\mathrm{s}}}}$, which is likely caused by the fact that there are fewer voxels and hence less information in the $z$-direction.

Figure 8 shows stochastic simulations from the models (4) and (2) with parameters chosen as the posterior expectation $E({\boldsymbol{\mathrm{\gamma}}}|\mathbf{y})$. Comparing these to the binarized CLSM image and the filtered CLSM pore structure, which are also shown in the figure, it seems like the model simulations are similar to the data.

Figures 9(a) and 9(b) show that for the fitted noise model (4), the empirical covariances for the binarized CLSM image goes above the model envelopes at small distances. The model marginal covariance functions were computed using the posterior mean of the parameters. The variances in the model and in the empirical variance of the CLSM image are almost the same, and hence this is caused by the variance of the GMRF in the fitted model (4) being too small compared to the variance of the noise. In other words, the signal to noise-ratio is underestimated in the fitted model. However, Figures 9(c) and 9(d) show that the empirical covariances of the CLSM pore structure falls within the model envelopes (although with a small margin for the $z$-line) for the fitted pore model (2). Hence the underestimation of the signal to noise-ratio does not have a big effect on the pore model.

The empirical size densities in the pore space and in the matrix, shown in Figure 10, tend to lie within the model envelopes, showing that the geometry of the CLSM pore structure and the stochastic simulations from model (2) are similar. Comparing the size densities with respect to the different structuring elements, there is room for longer lines along the $z$-axis than along the other two axes, and the size densities with respect to the lines are naturally larger than those for the sphere. Comparing the size densities of the pore space and the matrix, it is clear that we can fit larger structuring elements in the matrix than in the pore space. This is not surprising since we have a smaller volume fraction of pores than of matrix. The high values of the size density for the largest voxel in Panels (d) and (h) are caused by boundary effects. Since the size of a sample (the CLSM pore structure or stochastic simulations) is $74\times 74\times 20$, local sizes with respect to the line aligned with the $z$-axis cannot be larger than $10$. But if we think of the sample as being part of a larger pore structure, we cannot tell if the local sizes that are found to be $10$ in fact are larger. For further discussion of this, see Appendix B.

Figure 11 shows that the effective diffusion coefficient of the diffusion computed using Gesualdo in the CLSM pore structure lies well within the range of the effective diffusion coefficients computed in $1400$ stochastic simulations from model (2). The figure also shows the estimated negative and positive excursion sets $E_{0.01,0.01}^{-}$ and $E_{0.01,0.01}^{+}$, which were estimated from the model diffusive flux histograms corresponding to the same $1400$ stochastic pore structures. In the figure one can see that the inclusions $E_{0.01,0.01}^{+}\subset E_{\text{CLSM},0.01}\subset(E_{0.01,0.01}^{-})^{c}$ hold, which means that the set where the CLSM histogram takes values above $1\%$ includes the set where for most model histograms all values are above $1\%$. Also, no bins with values above $1\%$ in the CLSM histogram are included in the set where for most model histograms all values are below $1\%$. This indicates that the diffusive flux computed in the CLSM pore structure is similar to the diffusive flux computed in the stochastic pore structures. One can also note that the estimated excursion sets for the model are symmetric in $x$ and $y$, which is explained by the fact that the model is isotropic in the ${\boldsymbol{\mathrm{s}}}$-plane.

Discussion

We have formulated a parametric stochastic model for the pore structure of films of EC/HPC polymer blends, with the aim to analyze the mass transport through the films. The model is based on a separable oscillating Gaussian Matérn field, which is stationary and anisotropic. The field is obtained as the solution to an SPDE, which is solved using the finite element method giving an approximation of the Gaussian field by a GMRF. The model was fitted to CLSM images using an MCMC algorithm which takes advantage of the sparsity of the GMRF approximation and uses the Kronecker structure of its precision matrix in a, for this type of problem, new way to reduce the complexity of the algorithm. The computational efficiency of this model fitting algorithm allowed fitting the model to CLSM images of realistic sizes. We characterized the stochastic simulations from the model using the covariance and pore size distributions, and developed a measure using excursion sets to characterize the diffusive flux computed in the stochastic simulations. Using these measures we concluded that the model fits the CLSM images well.

A strength of the GMRF approximation is that it is not constrained to the triangulation that was used to estimate the parameters in the GMRF. Instead a finer triangulation can be used for the stochastic simulation from the model. As long as the same parameters are used, the model will be a linear approximation of the same continuously indexed Gaussian field. This is useful, since it is more computationally demanding to fit the model than to simulate from the model.

We use a simple filter to obtain the pore structure of model (2). A better alternative could be to use the posterior mean of the pore structure given the data, obtained from the MCMC algorithm, or to use a filter that e.g. takes the point spread function into account. However, the filter has no effect on the model fit since the model is fitted to model (4). Also, since the binarized CLSM data and stochastic simulations from model (4) are filtered in the same way, the filter should not have a large effect on the model validation either.

An important aspect to consider is what size of the observation window is sufficient to get good parameter estimates, i.e. how big the so called representative volume element is (Chiu, Stoyan, Kendall, and Mecke, 2013, pp. 236–237). In our case, the biggest concern is whether we have enough data in the $z$-direction, since we use rather thin samples to ensure stationarity. The range of the correlation in the $z$-direction is large compared to the size of the samples—as can be seen by comparing the correlation for the fitted model in the $z$-direction, with the size of the CLSM image in the $z$-direction, which is $20$ voxels. To see if the parameter estimates are sensitive to the size of the CLSM image sample in the $z$-direction, we fitted the model to the combination of two in the $z$-direction consecutive samples, both of size $20$ voxels in the $z$-direction, where the two samples have different pore sizes. With this larger combined sample, the point-estimates of all the parameters—and particularly the range parameters $\tau_{{\boldsymbol{\mathrm{s}}}}$ and $\tau_{z}$ which control the size of the pores in the ${\boldsymbol{\mathrm{s}}}$-plane and $z$-direction respectively—turned out to lie between the point-estimates for the model fitted to the two samples individually. Thus the estimate of the size of the pores in the combined sample lies somewhere between the estimates of the individual samples, which indicates that the observation windows in the individual samples are large enough for us to obtain valid estimates of the range parameters.

Despite the different techniques that were used to decrease the computational effort, the model fitting method is still computationally demanding. Other less demanding methods, such as the minimum contrast method, where the parameters are chosen to minimize the distance between a summarizing function and the estimate of the function computed from the data, and the pairwise likelihood have previously been used for this type of problems, see Wilson and Nott (2001). An advantage with the likelihood-based MCMC algorithm is that it uses more of the information in the data than the less computationally demanding methods. It would be interesting to perform a comparison of the parameter estimates from the MCMC algorithm with estimates from other model fitting methods, similar to what was done in Sidén, Eklund, Bolin, and Villani (2017).

Following Lindgren, Rue, and Lindström (2011), the model can be made non-stationary by letting the parameters of (3) be slowly varying functions, which is useful for capturing non-stationarities in the EC/HPC films. The model can also be extended using nested SPDEs (Bolin and Lindgren, 2011) to capture other types of anisotropies in the data. These extensions are interesting for future work where the model will be used to investigate how the mass transport depends on the microstructure of the material.

Acknowledgements

This work is funded by the Swedish Foundation for Strategic Research (SSF grant AM13-0066), the Knut and Alice Wallenberg foundation (KAW grant 20012.0067), and the Swedish Research Council (grant 2016-04187). We are grateful to Holger Rootzén for modelling ideas and valuable suggestions. We would also like to thank Henrike Häbel, Mariagrazia Marucci and Christine Boissier for supplying the data, and Jonas Wallin for sharing code that was used in the estimation method. Finally, we thank Aila Särkkä, Christian von Corswant, and everyone involved in the SSF project for valuable discussions.

Appendix A Details of the model and the model fitting algorithm

The triangulation used in the GMRF approximation

The triangulations of the domains $\Omega_{{\boldsymbol{\mathrm{s}}}}$ and $\Omega_{z}$ have a node in each voxel of the CLSM image. The domains are larger than the domain of the CLSM voxels, to reduce the effect of the boundary conditions imposed on the SPDE (see Lindgren, Rue, and Lindström, 2011). For an example, see Figure 12. The interior nodes in the dense subset of the triangulation of $\Omega_{{\boldsymbol{\mathrm{s}}}}$ are the microscopy data voxels. The exterior nodes are less densely spaced for computational efficiency. Because $\Omega_{z}$ is an interval, the $z$-triangulation is simply a division of the domain into subintervals. However, instead of using the model on the full domain $\Omega_{z}$, we consider the marginal distribution of only the interior nodes. This allows us to reduce the effect of the boundary conditions at a lower computational cost. In practice, this means that we use the precision matrix ${\boldsymbol{\mathrm{Q}}}_{z,\text{marg}}={\boldsymbol{\mathrm{Q}}}_{z,II}-{\boldsymbol{\mathrm{Q}}}_{z,IE}{\boldsymbol{\mathrm{Q}}}_{z,EE}^{-1}{\boldsymbol{\mathrm{Q}}}_{z,EI}$, instead of ${\boldsymbol{\mathrm{Q}}}_{z}$. Here,

[TABLE]

where $I$ denotes the indices that correspond to the interior nodes, and $E$ denotes the indices that correspond to the exterior nodes. Finding the marginal precision matrix ${\boldsymbol{\mathrm{Q}}}_{z,\text{marg}}$ can be done efficiently since ${\boldsymbol{\mathrm{Q}}}_{z,EE}$ is low-dimensional, but is too computationally demanding for ${\boldsymbol{\mathrm{Q}}}_{{\boldsymbol{\mathrm{s}}}}$.

The MCMC algorithm

The MCMC algorithm used to to estimate the model parameters is explained in more detail here.

To simplify the sampling of $\mathbf{w}$, we introduce the auxiliary variables $\mathbf{s}={\boldsymbol{\mathrm{A}}}\mathbf{w}+\boldsymbol{\epsilon}$. Then the distribution of $\mathbf{w}\mid\mathbf{s},\boldsymbol{\gamma},u,\mathbf{y}$ is multivariate Gaussian, which is easy to simulate from. The MCMC algorithm uses a Metropolis-within-Gibbs sampler, from which we obtain a Markov chain $\{(\mathbf{w}^{(i)},\mathbf{s}^{(i)},\boldsymbol{\gamma}^{(i)},u^{(i)})\}$ which has the desired posterior density as stationary distribution. The variables and parameters are updated in three separate blocks, which makes simulation of the vectors $\mathbf{w}$ and $\mathbf{s}$ easier while maintaining a low autocorrelation of the Markov Chain. A similar Gibbs sampler for binary data can be found in Albert and Chib (1993).

We start by selecting starting values $\{(\mathbf{w}^{(0)},\mathbf{s}^{(0)},\boldsymbol{\gamma}^{(0)},u^{(0)})\}$, and then repeat the following three steps for $i=1,2,\ldots$.

Sample $\mathbf{w}^{(i+1)}\sim\pi(\mathbf{w}\mid\mathbf{s}^{(i)},\boldsymbol{\gamma}^{(i)},u^{(i)},\mathbf{y})$ for $\mathbf{w}\mid\mathbf{s},\boldsymbol{\gamma},u,\mathbf{y}\sim\mathsf{N}\left(\hat{{\boldsymbol{\mathrm{Q}}}}^{-1}{\boldsymbol{\mathrm{A}}}^{\top}\mathbf{s}/\sigma^{2},\hat{{\boldsymbol{\mathrm{Q}}}}^{-1}\right)$ and $\hat{{\boldsymbol{\mathrm{Q}}}}={\boldsymbol{\mathrm{Q}}}(\boldsymbol{\gamma})+{\boldsymbol{\mathrm{A}}}^{\top}{\boldsymbol{\mathrm{A}}}/\sigma^{2}$. This step requires sampling a multivariate GMRF, which is the most computationally expensive step of the estimation procedure. How this is done is explained in the next section. 2. 2.

Sample $\mathbf{s}^{(i+1)},u^{(i+1)}\sim\pi(\mathbf{s},u\mid\mathbf{w}^{(i+1)}\boldsymbol{\gamma}^{(i)},\mathbf{y})$ using a Metropolis Hastings step. For this, a new threshold $u^{\prime}$ is proposed using a random walk proposal, and given the proposed threshold the proposal density for $\mathbf{s}$ is

[TABLE]

which is a product of independent truncated Gaussians. Here ${1}_{I}(s_{i})=1$ if $s_{i}\in I$, and ${1}_{I}(s_{i})=0$ otherwise. Finally ${\boldsymbol{\mathrm{s}}}$ and $u^{\prime}$ are accepted or rejected jointly. 3. 3.

Sample $\boldsymbol{\gamma}^{(i+1)}\sim\pi(\boldsymbol{\gamma}\mid\mathbf{w}^{(i+1)},\mathbf{s}^{(i+1)},u^{(i+1)},\mathbf{y})$ using a Metropolis Hastings step where each parameter is proposed separately. Constrained random walk proposals are used for the parameters $\theta_{{\boldsymbol{\mathrm{s}}}}$ and $\theta_{z}$, and log-normal proposals for the parameters $\kappa_{{\boldsymbol{\mathrm{s}}}}$ and $\kappa_{z}$. With the proposed parameters $(\theta_{{\boldsymbol{\mathrm{s}}}}^{\prime},\theta_{z}^{\prime},\kappa_{{\boldsymbol{\mathrm{s}}}}^{\prime},\kappa_{z}^{\prime})$, the gamma conditional density

[TABLE]

is then used as the proposal density for $\tau$. The parameters are accepted or rejected jointly.

Reducing the computational complexity of the MCMC algorithm

The most computationally demanding part of the MCMC algorithm is to sample from $\mathbf{w}\mid\mathbf{s},\boldsymbol{\gamma},u,\mathbf{y}$. We use the method introduced in Papandreou and Yuille (2010), and used in Sidén, Eklund, Bolin, and Villani (2017) to sample large GMRFs without having to compute the full Cholesky factor of $\hat{{\boldsymbol{\mathrm{Q}}}}$, with the following two steps:

Generate $\boldsymbol{\xi}=\tau({\boldsymbol{\mathrm{R}}}_{z}\otimes{\boldsymbol{\mathrm{R}}}_{{\boldsymbol{\mathrm{s}}}})^{\top}\mathbf{z}_{1}+{\boldsymbol{\mathrm{A}}}^{\top}\mathbf{z}_{2}/\sigma+{\boldsymbol{\mathrm{A}}}^{\top}\mathbf{s}/\sigma^{2}$, where $\mathbf{z}_{1},\mathbf{z}_{2}\sim\mathsf{N}(\mathbf{0},{\boldsymbol{\mathrm{I}}})$ are independent, and ${\boldsymbol{\mathrm{R}}}_{z}$ and ${\boldsymbol{\mathrm{R}}}_{{\boldsymbol{\mathrm{s}}}}$ are the upper triangular Cholesky factors of ${\boldsymbol{\mathrm{Q}}}_{z}$ and ${\boldsymbol{\mathrm{Q}}}_{{\boldsymbol{\mathrm{s}}}}$ respectively. By the properties of the Kronecker product, $\tau{\boldsymbol{\mathrm{R}}}_{z}\otimes{\boldsymbol{\mathrm{R}}}_{{\boldsymbol{\mathrm{s}}}}$ is the Cholesky factor of ${\boldsymbol{\mathrm{Q}}}(\boldsymbol{\gamma})$. All three matrix-vector products can be computed without evaluating the full matrices. The first product is computed as in (4), and the other two are computed using that ${\boldsymbol{\mathrm{A}}}$ is a binary matrix with one non-zero element per row. 2. 2.

Compute $\mathbf{w}$ by solving $\hat{{\boldsymbol{\mathrm{Q}}}}\mathbf{w}=\boldsymbol{\xi}$ approximately using a preconditioned conjugate gradient (PCG) method (Demmel, 1997).

The PCG-method in the second step solves the system ${\boldsymbol{\mathrm{P}}}^{-1}\hat{{\boldsymbol{\mathrm{Q}}}}{\boldsymbol{\mathrm{w}}}={\boldsymbol{\mathrm{P}}}^{-1}\boldsymbol{\xi}$ approximately where ${\boldsymbol{\mathrm{P}}}={\boldsymbol{\mathrm{R}}}_{{\boldsymbol{\mathrm{P}}}}{\boldsymbol{\mathrm{R}}}_{{\boldsymbol{\mathrm{P}}}}^{T}$ is a positive definite and symmetric preconditioner. A good preconditioner approximates $\hat{{\boldsymbol{\mathrm{Q}}}}$ in such a way that ${\boldsymbol{\mathrm{P}}}^{-1}\hat{{\boldsymbol{\mathrm{Q}}}}$ is well-conditioned and ${\boldsymbol{\mathrm{P}}}^{-1}{\boldsymbol{\mathrm{v}}}$ can be computed efficiently. Sidén, Eklund, Bolin, and Villani (2017) used the preconditioner with ${\boldsymbol{\mathrm{R}}}_{{\boldsymbol{\mathrm{P}}}}=ichol(\hat{{\boldsymbol{\mathrm{Q}}}})$, where $ichol$ denotes an incomplete Cholesky factorization. We instead use the preconditioner with Cholesky factor $ichol({\boldsymbol{\mathrm{Q}}}_{z})\otimes ichol({\boldsymbol{\mathrm{Q}}}_{{\boldsymbol{\mathrm{s}}}})$, and use the structure of the preconditioner and of $\hat{{\boldsymbol{\mathrm{Q}}}}$ to compute ${\boldsymbol{\mathrm{P}}}^{-1}(\hat{{\boldsymbol{\mathrm{Q}}}}{\boldsymbol{\mathrm{w}}})$ using (4) and (5) without evaluating the full matrices. Table 1 compares these two choices with some other natural choices of preconditioners for sampling $\mathbf{w}$, and shows that using $ichol({\boldsymbol{\mathrm{Q}}}_{z})$ and $ichol({\boldsymbol{\mathrm{Q}}}_{{\boldsymbol{\mathrm{s}}}})$ is better than using ${\boldsymbol{\mathrm{R}}}_{z}$ and ${\boldsymbol{\mathrm{R}}}_{{\boldsymbol{\mathrm{s}}}}$ in this case, since it leads to a more sparse preconditioner. From the results in the table, it is also clear that it is better to avoid evaluating the full Cholesky factor of the preconditioner by using its Kronecker structure, if possible.

Using (4) and (5) reduces the amount of elements we need to store in the MCMC algorithm from $\mathcal{O}((n_{z}n_{{\boldsymbol{\mathrm{s}}}})^{2})$ to $\mathcal{O}(n_{z}^{2}+n_{{\boldsymbol{\mathrm{s}}}}^{2})$, where ${\boldsymbol{\mathrm{Q}}}_{z}\in\mathbb{R}^{n_{z}\times n_{z}}$ and ${\boldsymbol{\mathrm{Q}}}_{{\boldsymbol{\mathrm{s}}}}\in\mathbb{R}^{n_{{\boldsymbol{\mathrm{s}}}}\times n_{{\boldsymbol{\mathrm{s}}}}}$—not taking the sparsity of ${\boldsymbol{\mathrm{Q}}}_{z}$ and ${\boldsymbol{\mathrm{Q}}}_{{\boldsymbol{\mathrm{s}}}}$ into account—since we never store full Kronecker products. As can be seen in Table 1, the gain from using (5) can be substantial for the PCG step, which is the most computationally demanding part of the algorithm. Another potentially demanding part is the computation of the Cholesky factor of ${\boldsymbol{\mathrm{Q}}}(\boldsymbol{\gamma})$, which is done to evaluate $\varphi_{(\mathbf{0},{\boldsymbol{\mathrm{Q}}}^{-1}(\boldsymbol{\gamma}))}(\mathbf{w})$ when sampling the parameters $\boldsymbol{\gamma}$. The complexity of the Cholesky factorization of the sparse ${\boldsymbol{\mathrm{Q}}}_{z}$ is $\mathcal{O}(n_{z})$, and the complexity of the factorization of the sparse ${\boldsymbol{\mathrm{Q}}}_{{\boldsymbol{\mathrm{s}}}}$ with an optimal reordering is $\mathcal{O}(n_{{\boldsymbol{\mathrm{s}}}}^{3/2})$ (Rue and Held, 2005, Ch. 2.4). Compare this with the cost of factorizing a dense $n\times n$ matrix, which is $\mathcal{O}(n^{3})$. Performing the Cholesky factorization on the full matrix ${\boldsymbol{\mathrm{Q}}}(\boldsymbol{\gamma})$ would therefore take $\mathcal{O}((n_{z}n_{{\boldsymbol{\mathrm{s}}}})^{3})$ operations if ${\boldsymbol{\mathrm{Q}}}(\boldsymbol{\gamma})$ were dense, whereas the lower bound using the Kronecker structure of its Cholesky factor $\tau{\boldsymbol{\mathrm{R}}}_{z}\otimes{\boldsymbol{\mathrm{R}}}_{{\boldsymbol{\mathrm{s}}}}$ is $\mathcal{O}(n_{z}+n_{{\boldsymbol{\mathrm{s}}}}^{3/2})$ operations, taking into account the sparsity.

Thus it is clear that the algorithmic complexity of the MCMC algorithm is greatly reduced by using the sparsity and Kronecker product structure of the precision matrix ${\boldsymbol{\mathrm{Q}}}(\boldsymbol{\gamma})$ and its Cholesky factor. This allows us to fit the model to much larger data sets. Using (5) we can also generate samples from the GMRF ${\boldsymbol{\mathrm{w}}}$ efficiently.

Appendix B Estimation of the size distributions

The local size in a point in a (random) set $\Xi$ with respect to a structuring element $B$ can be defined using the morphological operation opening, denoted $\circ$. The opening of $\Xi$ by $B$ consists of an erosion followed by a Minkowski addition with $B$. The opened set can be written

[TABLE]

The local size $h$ is defined as $h(x,\Xi,B)=\sup\{\lambda:x\in\Xi\circ\lambda B\}$, where $\lambda B=\{\lambda x:x\in B\}$. It follows that, if the supremum can be attained, the local size in a point $x\in\Xi$ is the largest value $\lambda$ for which the rescaled structuring element $\lambda B$ can be translated so that it fits within the set $\Xi$ and covers the point $x$.

We define the size distribution for a random set $\Xi$ with respect to a convex structuring element $B$ as

[TABLE]

i.e. the size distribution $S(\lambda;x,\Xi,B)$ gives the probability that the local size with respect to $B$ in a point $x$ in $\Xi$ is greater than or equal to $\lambda$. Since the structuring element is convex, $h(x,\Xi,B)\geq\lambda$ is equivalent to $x\in\Xi\circ\lambda B$, and so the size distribution can be written as $S(\lambda;x,\Xi,B)=P(x\in\Xi\circ\lambda B\mid x\in\Xi)=P(x\in\Xi\circ\lambda B)/P(x\in\Xi)$ (Matheron, 1975; Delfiner, 1972).

When estimating the size distribution of a stationary random set $\Xi$ defined on $\Omega$, which is only observed in a bounded window $W\subset\Omega$, we have to consider the boundary effects to get an unbiased estimate. If we do not have any information about the set outside of $W$, we can in general only observe the opened set $\Xi\circ\lambda B$ in the smaller window $\tilde{W}(\lambda B)=(W\ominus\lambda B)\ominus-\lambda B\subset\lambda B$ (Ohser and Schladitz, 2009). This means that we cannot determine the local size in a point close to the boundary of the observation window $W$, since we need information about the surrounding of the point. To account for this, we can use a minus-sampling estimator of the size distribution based on the unbiased estimate $v_{d}(\Xi\circ\lambda B\cap\tilde{W}(\lambda B))/v_{d}(\tilde{W}(\lambda B))$ of $P(x\in\Xi\circ\lambda B)$, where $v_{d}(\cdot)$ denotes the volume of the set. Properties of minus-sampling estimators for size distributions for random closed sets can be found in Moore and Archambault (1991), and properties for discretely indexed random sets in Sivakumar and Goutsias (1997).

Because we select microscopy data samples that have homogeneous pore sizes, our observation windows are small in the $z$-direction. Therefore, we do not use the unbiased minus-sampling estimator since that would reduce the observation window of the sample further. Instead we approximate $\Xi\circ\lambda B\cap W$ with the opened set $(\Xi\cap W)\circ\lambda B$, which is equivalent to assuming that there are no pores outside of the observation window $W$. This gives us a biased estimate $v_{d}((\Xi\cap W)\circ\lambda B)/v_{d}(W)$ of $P(x\in\Xi\circ\lambda B)$. However, the bias has the same effect for the filtered CLSM pore structure and the stochastic simulations from the model, since the observation windows are of the same size, and hence will not affect comparisons.

The size distribution for a structuring element $B$ is thus estimated as

[TABLE]

for $\lambda=1,2,\dots,N$, where the maximal local size $N$ is determined by the observation window $W$ in which the set $\Xi$ is observed. The size density is estimated as $\hat{s}(\lambda)=\hat{S}(\lambda)-\hat{S}(\lambda+1)$.

Appendix C Results for all CLSM samples

Here we present results for the model fitted to all CLSM samples. More results for sample $HPC30_{1}$ can be found in Section “Results”. For some samples, only every second or third voxel in the CLSM image was kept, depending on the average size of the pores. There is also a resolution difference between the $30\%$ and $40\%$ CLSM images. Thus, the distances between voxels are not the same in all samples (see Table 2). However, we chose the size of the domains $\Omega$ so that the fitted parameters are comparable between the samples. The estimated oscillation parameters are fairly similar for all model fits (Table 3).

A larger range parameter for the model corresponds to a smaller covariance range and hence smaller pores. The samples for the two CLSM images were ordered so that the first two samples had smaller pore sizes than the last two samples for each image, which can be seen in the estimated range parameters. One can also note that the estimated range parameters are larger for the samples with $40\%$ HPC than for the samples with $30\%$ HPC.

The differences between the sample and fitted model volume fractions are no more than two percent, and the differences between the sample and fitted model surface areas are $3$–$13\%$, where the surface areas were calculated using the algorithm provided by Legland, Kiên, and Devaux (2011). The volume fractions of the samples differ from the theoretical HPC weight ratios. The depth dependent threshold was always determined using a larger part of the CLSM image than the sample itself, and there may e.g. be edge effects or dust particles outside of the sample contributing to the volume fraction discrepancies.

The covariance functions and size distributions (Figures 13–19) for the remaining samples show that the model provides a relatively good fit. It is clear that the signal to noise ratio of the covariance functions corresponding to the noise model (1) is underestimated (Panels (a)–(b)). However, as was the case for sample $HPC30_{1}$, this underestimation does not have a big effect on the pore model.

The marginal covariance functions in the ${\boldsymbol{\mathrm{s}}}$-plane, $C_{{\boldsymbol{\mathrm{s}}}}$, for both the noise and pore model, fit better at smaller distances (Panels (a) and (c)). For half of the samples, the ranges of the estimated covariances $\hat{C}_{{\boldsymbol{\mathrm{s}}}}$ are slightly underestimated, although the shapes of $\hat{C}_{{\boldsymbol{\mathrm{s}}}}$ are consistent with the model envelopes (Figures 13, 15–17, Panels (a) and (c)). A possible explanation comes from the fact that even though the samples were kept small in the $z$-direction, the pore sizes decrease slightly along the $z$-axis from the bottom to the top of each sample. The estimated range parameters could then have been chosen to better fit the part of the sample with smaller pore sizes. For samples with high estimated oscillation parameters $\theta_{{\boldsymbol{\mathrm{s}}}}$, the smaller peaks and troughs of the covariances $\hat{C}_{{\boldsymbol{\mathrm{s}}}}$ do not seem to be as pronounced as those of the model envelopes (Figures 15–18, Panels (a) and (c)). This might be because the peaks and troughs shift and therefore cancel out due to the decreasing pore size along the $z$-axis. The size distributions lie within their envelopes to a higher extent than the covariance functions (Panels (e)–(l)). Thus the pore shapes of stochastic simulations from the fitted pore models correspond well to the pore shapes of the pore structures obtained from the CLSM samples.