Pore space extraction and characterization of sack paper using μ‐CT

We show that attenuation X‐ray microcomputed tomography (μ‐CT) offers a route to extract the three‐dimensional pore space of paper reliably enough to distinguish samples of the same kind of paper. Here, we consider two sack kraft papers for cement bags with different basis weights and thicknesses. Sample areas of approximately 5 mm2 with a resolution of 1.5 μm are considered, i.e. sizes that exceed sample areas of 2 mm2 for which the pore structure was previously studied in the literature. The image segmentation is based on indicator kriging as a local method that removes ambiguities in assigning voxels as pore or as fibre. The microstructures of the two samples are statistically compared in terms of descriptors such as sheet thickness, porosity, fractions of externally accessible pores and mean geodesic tortuosity. We demonstrate that a quantitative comparison of samples in terms of porosity and thickness requires a common definition of the sheet surfaces. Finally, the statistical pore space analysis based on the μ‐CT scans reliably reveals structural differences between the two paper samples, but only when several descriptors are used.


Introduction
There is tremendous interest to comprehend the pore structure of paper, because the porosity of paper is a key material property that enables or crucially supports a wide range of applications. Amongst the methods to experimentally determine the three-dimensional (3D) microstructure of paper (see Chinga-Carrasco, 2009 for a comprehensive list of references), X-ray microcomputed tomography (μ-CT) is a well-established nondestructive method, in particular to study the geometry of the * Equally contributing authors.
However, extracting and statistically analysing the pore space of paper still poses unique, paper-born challenges. These challenges root in two factors: On the one hand, sample sizes are required to be as large as possible, because paper fibres are highly nonuniform in shape. On the other hand, paper pores span sizes from more than 50 μm down to a few nanometres; while the larger pores are predominantly formed between fibres, the nanometre-size pores are typically nested in the fibre walls. The latter ones are closed in dry paper and open only when water penetrates the fibre (Stone et al., 1966).
To resolve the pores, the microstructure acquisition method ought to resolve features as small as possible. The resolution of a given method is defined as the smallest distance at which two objects can be distinguished. In μ-CT, the resolution is, at best, as small as the voxel edge length. In fact, resolutions as low as 0.7 μm were reported for 3D imaging of paper (Chinga-Carrasco et al., 2008;Axelsson & Svensson, 2010); however, at the cost of rather small sample areas ranging between 0.4 and 2 mm 2 (Axelsson et al., 2006;Rolland du Roscoat et al., 2007;Chinga-Carrasco et al., 2008;Axelsson & Svensson, 2010;Sharma et al., 2015).
In the present work, we pursue the question whether it is possible to achieve a reasonable compromise between sample size and resolution when employing μ-CT. Holmstad et al. (2006) provided evidence that a mediocre resolution of 1 μm in phase-contrast μ-CT preserves the topology of the pore space. In going to an even lower resolution, we aim at characterizing paper samples of approximately 5 mm 2 in area with attenuation μ-CT. Attenuation μ-CT yields absorption contrast images (Rolland du Roscoat & Bloch, 2005) and can be carried out for such sample sizes with a laboratory desktop setup. We particularly focus on (i) whether we can reliably extract the pore space from the scanned data, (ii) whether the statistical analysis of the acquired pore space allows us to distinguish two different paper sheets and (iii) to which extent a statistical pore space analysis corroborates the ranking of the paper sheets obtained from standard paper characterizations. Holmstad et al. (2006) and Chinga-Carrasco et al. (2008) demonstrated that the analysis of the pore structure acquired with μ-CT allows to distinguish paper types as different as filter paper, newsprint and handsheets. To challenge our acquisition and analysis methods, we intentionally investigate two samples of the same kind of paper, i.e. samples that share the same fibre type, but differ markedly in grammage (i.e. basis weight) and thickness. Specifically, we consider two types of sack kraft paper that are employed to produce cement bags. They combine a high porosity with a superior mechanical strength.
Last but not least, we show that pores and solid materials can be distinguished in the 3D scan of mediocre resolution by applying a local thresholding method. The corresponding segmentation method relies on the determination of a pair of thresholds followed by indicator kriging (Oh & Lindquist, 1999). In a subsequent analysis based on tools from spatial statistics (Chiu et al., 2013) and mathematical morphology (Matheron, 1975), we quantify the pore structure in terms of sheet thickness, porosity, the spherical contact distribution function, and the pore-network characterizing quantities such as fraction of accessible pores and tortuosity. In particular, we explore the question to which extent ambiguities in defining paper sheet surfaces affect the paper thicknesses and subsequent geometrical pore-related properties.

Samples
Two paper samples of 3 × 3 mm 2 are considered in the present investigation based on 3D imaging, denoted by sample I and sample II. The paper types of both samples are made of virgin unbleached kraft pulp, and possess different basis weights, i.e. according to the supplier, nominally 70 g/m 2 for sample I and 110 g/m 2 for sample II. In particular, sample II is taken from a just-filled valve cement sack kraft paper. Since the cellulose is known to yield a low contrast in attenuation μ-CT, the presence of impurities, such as cement grains, containing elements heavier than carbon and oxygen is expected to further lower the contrast between the cellulose fibres and the air-filled pore space. This will aggravate the expected issues to distinguish pore space and cellulose.

Acquisition, pre-and postprocessing of μ-CT data
The μ-CT scans were performed at the company MITOS GmbH (MITOS, https://www.mitos.de) using the Xradia 500 Versa 3D X-ray microscope (Zeiss, Germany). In this work, an X-ray tube voltage of 60 kV was used without an additional X-ray filter. For each sample, 2401 projections over 360 • of rotation were recorded. This system provides a two-stage magnification. In the first one, the projection image is enlarged through geometric magnification. For this case, the source-to-sample and sample-to-scintillator distances were 17 and 20 mm, respectively. In the second stage, the scintillator converts X-rays to visible light, which is optically magnified and then captured by a CCD camera, with 2048 × 2048 pixels. Specifically a 4× objective lens was used. The data were reconstructed from the projections using the standard Feldkamp-Davis-Kreiss algorithm (Feldkamp et al., 1984). No additional efforts were directed to enhance the resolution by exploiting phase contrasts or to suppress phase artefacts. The Versa system is designed so that even at high resolutions, the spatial resolution corresponds to the voxel size. In the present case, the voxel edge length is 1.5 μm. After reconstruction, two sets of 1983 (1984) raw images (16-bit-TIFF) are obtained for samples I and II, respectively.
As mentioned before, phase artefacts were not removed during reconstruction. Phase artefacts lead to an artificial edge enhancement (Wernersson et al., 2013). Correspondingly, they can result in artificially increased grey values of the cellulose and decreased values related to pores. This impacts a subsequent porosity analysis based on the image histograms. We reduce the phase artefacts from the reconstructed images by applying a spatial phase-removal 3D filter, as proposed by Wernersson et al. (2013). This filter yields results comparable to common single-image in-line phase-retrieval and phase-correction algorithms. Moreover, it does not require the projection data. The filter requires two parameters, σ and c, which relate the amplitude and width of the fringes in images to the mix of absorption and phase contrast. In the present work, we choose σ = 1.975 and c = −29.3 for both samples, for a more detailed description about how to fit these parameters, we refer to the paper of Wernersson et al. (2013). The such phase-corrected image data are then filtered by applying a Gaussian low-pass filter with standard deviation 0.75.
During the scan, the sample is placed in a holder, whose shape results in a slight bend of the paper sheet. This bend can be observed in two raw images shown in Figure S1 in the Supporting Information (SI). We remove this slight bend of the paper with a geometric transformation to facilitate the identification of the sheet surface. Details of this algorithm are included in Figure S2. This transformation does not affect the porosity; we find a vanishingly small difference in the calculated porosity for preliminarily segmented images prior and after unbending. In a last step, each image is cropped to the area containing the paper sample, thus discarding regions solely containing the sample holder. To this end, the exact position of the edges of the sample holder is determined using an edge-detection algorithm (Canny, 1986) on the corresponding enhanced image. Details of the separation of the region of interest from the holder and background are given in Figure S2 in the SI.
The upper and bottom parts of the sample are not always in the field of view of the cone beam during sample rotation. To unambiguously remove images affected by an incomplete field of view, we inspected the global thresholds in the subsequent segmentation step (cf. 'Determination of thresholds' section). In brief, the segmentation method identifies two global thresholds, i.e. specific grey values, considering all images in the 3D data set. When we repeat the threshold determination for each image individually, we find that these image-related thresholds correspond to the globally obtained thresholds for the bulk of the images. However, images located at the border of the reconstructed volume, i.e., images likely to be affected an incomplete projection, possess an upper threshold that markedly differs from the global threshold. Such images were discarded. Hence, the considered final sample areas are 5.6 mm 2 (2.001 × 2.802) for sample I and 4.7 mm 2 (1.666 × 2.812) for sample II .

Experimental paper characterization
The thickness of the paper sheets is determined according to the standard DIN EN ISO 534 (Paper and board: Determination of thickness, density and specific volume), whereas the grammage (i.e. basis weight) is determined according to DIN EN ISO 536 (Paper and board: Determination of grammage). Specifically, the paper thickness is measured with the mechanical contact method thickness tester Lehmann LDAL -03. In this measurement, the paper sheet is placed between two metal plates (Schultz-Eklund et al., 1992). The metal plates are pressed together until a certain, standardized pressure is reached. Then, the distance between the plates is a descriptor for the apparent thickness of the paper. Note that it is possible that the paper sheet is slightly compressed during this measurement. The resulting values for the thicknesses are 106±4 μm in sample I and 151±3 μm in sample II. In the following, we will refer to an apparent thickness whenever a thickness has been obtained with such a caliper-based measurement. The basis weights amount to 70±1.5 g/m 2 in sample I and 110±3.3 g/m 2 for sample II. We estimate a common, maximum relative error of the basis weight of ca. 3%, that consists of the largest assessable relative errors from the area measurement (≤2%) and the weight determination (≤1%). The ratio between the basis weight and the thickness yields the mass density of the paper sheets, i.e. 0.64 gcm −3 for sample I and 0.73 gcm −3 for sample II. Correspondingly, the maximum relative error of the mass density is 13%.

Segmentation
The data obtained from μ-CT measurements correspond to a 3D greyscale image that consists of stacked 2D images.
A key challenge to analyse greyscale images is to precisely distinguish between different materials within the samples (Rolland du Roscoat et al., 2012). In the segmentation step, each voxel with a given grey value is assigned to one of the present materials, in our case, a cellulose-based fibre, an airfilled pore or an impurity within the paper sheet. The success of any subsequent analysis and its interpretation strongly relies on the accuracy of this assignment step, because an incorrect segmentation distorts the size and the shape of the pore space.
Hence, the practical challenge is to find suitable criteria for such assignments. The definition of optimal criteria for the segmentation of greyscale images is still an active field of research. There are several methods available, for a review see, e.g. Wildenschild & Sheppard (2013). However, their individual strengths and weaknesses can be compared only for synthetic images, i.e. images in which each voxel is affiliated to a material in advance (Oh & Lindquist, 1999;Schlüter et al., 2014). Hence, the choice of a segmentation method is mainly determined by the type of material under study as well as by the quality and resolution of the image data at hand.
Below, we illustrate the main features of our chosen method as well as its specific application to our image data that represents the microstructure of paper sheets. For the sake of clarity, we illustrate the different steps with the help of an excerpt of only one two-dimensional cross-section image from the stack of sample I. The chosen image is shown in Figure 1(A). However, it is important to highlight that the segmentation of both samples is done in a 3D fashion. This implies that the analysis is made in terms of voxels rather than pixels. Figure 1(B) is a magnified excerpt of the image in Figure 1(A), in which the underlying heterogeneous structure is even more striking. The histogram of grey values for the complete stack of sample I is presented in Figure 1(C). In the same figure, we superimpose the corresponding entropy function as described in Kapur et al. (1985).
A first visual inspection of both, images and histogram, strongly suggests the presence of three distinct materials that we assign to cellulose-based fibres, pores and impurities. Each material contributes distinctly to the histogram of greyscale values in Figure 1(C). Pores are responsible for the cluster of greyscale values that peak at approximately 16 050 in Figure 1(C). A solid, mostly cellulose-based material causes a second cluster of grey values with a peak at approximately 17 700. These two clusters possess a large degree of overlap. This pattern of two strongly overlapping peaks is found in the histograms of both samples and seems to be a common feature of paper μ-CT scans (Rolland du Roscoat & Bloch, 2005). Voxels, whose grey values stem from the overlapping region, will be particularly difficult to assign to one of the materials, since these voxels are predominantly located at the fringes between pores and solid domains. The histogram in Figure 1(C) contains further, nonnegligible contributions at grey values that largely exceed those of the cellulose-based peak. We find such contributions in both samples I and II and relate them to impurities. Since sample II was taken from a just-filled cement bag, particularly plausible candidates for such impurities in that sample are cement particles.

Segmentation using indicator kriging
To segment the scale images, we perform an indicator kriging in the spirit of Oh & Lindquist (1999). Indicator kriging is a widely used and successful method for the segmentation of greyscale images (Ohser & Schladitz, 2009). However, to the best of our knowledge, this is the first time that indicator kriging is applied to paper samples. This thresholding method does not only rely on the analysis of the histogram of greyscale values, but also on information about spatial dependencies. It is a procedure that (i) assigns a population to the set of voxels which can be clearly identified by global thresholding and (ii) enhances the assignment accuracy for the remaining voxels while keeping the surfaces between regions of different materials smooth.
In a nutshell, the two steps of indicator kriging can be explained as follows: (i) The initial thresholding step, which is based on the analysis of the global histogram of the entire greyscale image. By the aid of this analysis, two thresholds T R , with R ∈ {pore, solid}, are determined. In absorption contrast mode, the greyscale values are proportional to X-ray attenuation, which is a function of the density and type of atoms. Hence, in a first step, a given threshold can be used to reliably label voxels. Thus, all those voxels with a grey value z with z ∈ [mi ni mum, T pore ] are assigned to pores, while those with z ∈ [T solid , ma xi mum] are considered as solid. All voxels with z ∈ (T pore , T solid ) remain unassigned. (ii) The second step comprises the segmentation of unassigned voxels. For each unassigned voxel at x o , the method estimates the probability that x o belongs to one of the two contributing materials. This estimation is based on the grey values of voxels located in the 3D neighbourhood of x o .

Determination of thresholds.
This section provides the initial thresholds T pore and T solid with the so-called entropy method. The determination of these initial thresholds has a decisive impact on the extent of voxel misidentification. Oh & Lindquist (1999) propose two methods of choosing the thresholds, the binormal mixture and the entropy method. The former is applicable when the distribution of the greyscale values in the histogram is due to two distinct contributions of approximately univariate normal shape. Unfortunately, however, each of the peaks in the histogram of Figure 1(C) cannot be described by a normal distribution. Specifically, by using an expectation maximization algorithm (Dempster et al., 1977;Benaglia et al., 2009) to only the grey values around both peaks, it is not possible to find a Gaussian-mixture model with only two components.
The entropy method proposed by Kapur et al. (1985) itself is a segmentation method based on the Shannon information theory (Shannon & Weaver, 1949) that entirely relies on the histogram information. Its goal is to find a single optimal threshold z * , which is defined as the grey level which maximizes the entropy function. If we look at Figure 1(C), the maximum of the entropy function stands out clearly at a grey value higher than those related with the second cluster in the histogram, i.e. the cellulose-based one. This becomes clear when having in mind that the Shannon formulation is based on the concept that the gain in information from an event is inversely related to its probability of occurrence (Lombardi et al., 2016). Thus, due to the rare-event character of the impurities on both samples, the maximum of the entropy function is related to them. However, on closer inspection, it emerges that the entropy function has a local maximum in the region between the peaks of the histogram related to pore and solid populations, see Figure 2(A). This feature is independent of the sample.
Based on this locally entropy maximizing grey value z * , we compute the thresholds T R , with R ∈ {pore, solid}, by where r e ∈ (0, 1) and ψ the entropy function (Oh & Lindquist, 1999). In the present work, a value of r e = 0.015 is used for both samples. This value keeps the thresholds in the region of the two peaks where the overlapping occurs. The calculated thresholds T pore and T solid are indicated by vertical lines in the histogram shown in Figure 2(A). All voxels with greyscale values below and above these thresholds are now assigned either to the pore (dark blue) or the solid population (lime green) as indicated in Figure 2(B). A necessary step before indicator kriging can be performed is to calculate the corresponding standard deviation for each of the thresholded populations, i.e. σ pore and σ solid . All unassigned voxels in Figure 2(B), which remain in the white region z ∈ (T pore , T solid ) of Figure 2(A), will be now treated by the indicator kriging.

Segmentation of unassigned voxels.
To give a clear idea, how we utilize the two thresholds T pore and T solid to label the remaining, unassigned voxels, we briefly explain some details of the indicator kriging method. For a more detailed description, we refer to the paper of Oh & Lindquist (1999). For an arbitrary unassigned voxel at x o we define a 3D 'spherical' neighbourhood N(x o ) = {x 1 , . . . , x 256 } of x o with constant radius 4, i.e. x ∈ N(x o ) if the Euclidean distance between x and x o is smaller or equal to 4, see Figure 3(A). The neighbourhood N(x o ), called the kriging window, encloses the set of neighbouring voxels, see Figure 3(B). The voxel x o will be labelled as pore or solid depending on max(P (T pore ; x o |n), 1 − P (T solid ; x o |n)), where is the likelihood that x o belongs to the population R ∈ {pore, solid}. The smoothed indicator kriging functionsî (T R ; x), see Figure 3(C), evaluate the proportion of voxels in the kriging window valued below/above a given threshold, for each x ∈ N(x o ) and each population R. Here, z(x) denotes the grey value of a voxel located at x, F (z) the cumulative distribution function estimated from the original greyscale values. Expressions s v R and s u R in Eq.
(3) correspond to s v pore = s u solid = 0 and s u pore = s v solid = (σ pore T solid + σ solid T pore )/(σ pore + σ solid ). The original greyscale values in the kriging window are transformed by the smooth indicator kriging into two sets with values between zero and one as shown in Figure 3(D). It is important to emphasize that the evaluation ofî (T pore ; x) and i (T solid ; x) requires the original greyscale values of the voxels inside the kriging window, no matter whether they have been already assigned during the inspection of a neighbouring voxel or in the initial thresholding run.
The weights λ x (T R ; x o ) are determined by the so-called ordinary kriging system as described in Oh & Lindquist (1999), i.e. the weights are determined based on the spatial correlation of the greyscale values. Note that the grey value z(x o ) does not play any role in the assignment of x o itself, but rather determines the identification of the neighbouring unassigned voxels. Once the voxel x o is assigned to one population, the kriging window is moved to the next unassigned voxel and so on. The process is repeated until each voxel of the entire stack is assigned to one of the two populations, as shown in Figure 3(E).

Impurity distinction
Once each image is segmented in pore and solid contributions, see Figure 4(A), we turn to discriminate between the cellulose and the high grey values that may appear in the solid regions. Due to the heterogeneous composition and to the rare-event character of the impurities, indicator kriging is neither suitable for separation of the solid-phase voxels into two new populations, i.e. cellulose and impurities, nor for the assignment of the boundaries between them. On the other hand, this argument validates the use of the global entropy-maximizing grey value z * as a single threshold between paper and impurities. Thus, each voxel with a grey value z ≥ z * is labelled as impurity.
However the spatial phase-removal 3D filter undesirably alters the neighbourhood of the impurities. This is reflected in an artificial increase of the grey values of voxels in the close vicinity of the impurities. If left untreated, this effect can suppress the porosity in these regions during the segmentation. Therefore, prior to the segmentation step, a local correction has been applied to the voxels in these regions. For this purpose, we make use of a growing-type algorithm, the so-called fast marching method (FMM) (Sethian, 1996) as implemented  in Matlab (MathWorks, Natick, MA, USA), to identify the set of voxels that are associated with these regions. The seed positions for the propagating fronts are the voxels already labelled as impurities. Having defined these regions, the grey values are corrected with the help of a transfer function. Further details about this correction are provided in Figure S5. Figure 4(B) shows the final segmentation of an image from sample II where the impurities are represented in yellow. For a complete visualization of the impurities in the volume, we refer the reader to Figure S6.

Pore space characterization
Having achieved a complete segmentation of the μ-CT data, we turn to the determination of the (inner) pore space and its characterization. After having separated intrinsic (inner) pore space from outer void space, we particularly pay attention to (i) whether we can discriminate samples I and II and (ii) whether the statistical analysis of the pore space corroborates the findings from the standard experimental paper characterization.

Determination of paper volume and pore space fractions
A major difficulty to determine the porosity of paper, in particular when starting from its microstructure, is the lack of a well-defined geometry of the paper sheets, which consist of a network of fibres and possess corrugated surfaces as well as significant variations in the thickness. These variations can be clearly seen in the volume reconstruction of the two paper samples, as illustrated in Figure 5. Without having clearly defined paper surfaces, neither the thickness, the total volume, nor the porosity of the sample can be unambiguously determined. In turn, one may expect that a given recipe to determine the surface is going to impact the porosity values and the qualitative trends when comparing both samples.

Determination of the paper thickness.
Due to the complex geometry of the paper surface, it is straightforward to accept that the thickness of a paper sheet is not a constant, but rather corresponds to a distribution of local thickness values. At each position of the paper sheet, the local thickness is said to be the distance between the top and bottom surface of the paper sheet. Hence, it is essential to define the paper surface within the segmented 3D image data. In the absence of a standard definition, we explore three plausible definitions for paper surfaces. Each of these scenarios mimics an experimental method to determine the thickness.

Rolling ball approach.
To take the complex details of the top and bottom surfaces into account, the rolling ball approach is a reasonable and intuitive method. This method has previously been used to determine the paper surface based on segmented μ-CT images (Holmstad et al., 2006;Chinga-Carrasco et al., 2008;Axelsson & Svensson, 2010). Moreover, the rolling ball method has been used to define the outer boundary of a certain type of coated paper based on scattering electron micrographs in Chinga & Helle (2002). This approach mimics the key concept of an experimental local thickness probe. Such a  setup consists of two metal spheres with the same diameter that, prior to a measurement, touch each other with a well-defined force. During the measurement, the paper is fed through these two spheres with a constant speed. The two spheres can be displaced along their connecting axis, however, with the initial force held constant. Then, the displacement of the spheres is a measure for the local thickness of the sheet (Schultz-Eklund et al., 1992;SCAN-P 88:01, 2001). In transferring this idea to our microstructures, we let a ball of radius r = 15 μm freely role within the void space in our segmented images as in Chinga-Carrasco et al. (2008). Formally, we perform a morphological closing (Soille, 2003) with radius r on the set of voxels being assigned either as cellulose material or as impurities. The resulting set of voxels is denoted by C . Each voxel belonging to a certain (discretized) box W is considered to be a part of the pore space if it is neither labelled as cellulose material nor as impurity and if one of the following conditions is fulfilled: 1. the voxel is contained in C ⊂ W and 2. the voxel cannot be connected to the upper or lower boundary of the image through the complement W \ C of C . Whenever connected clusters have to be computed from image data, the algorithm introduced in Hoshen & Kopelman (1976) is used. Roughly speaking, each void voxel is defined as a pore if it cannot be reached by the rolling ball from outside. The resulting surfaces are visualized in Figure 6(A). Due to the generally geometrically complex corrugation on both surfaces of the sheet, the thickness of the paper sheets becomes notably positiondependent (dark grey in Figure 6A).

Visibility depth approach.
In a second case, we take the rolling ball idea to its extreme. Here, we classify a voxel as surface voxel if we cannot find further solid voxels above (top surface) or below (bottom surface) along a vertical voxel column ( Figure 6B). This approach probes the surfaces in spirit of a sharp, idealized tip that is vertically approaching the paper sheet either from above or below. This translates into the following formal definition: Let x = (x 1 , x 2 , x 3 ) be an arbitrary voxel where x 1 , x 2 , x 3 are non-negative integer numbers. Moreover, a (discretized) box W (i.e. the discretization of a rectangular cuboid) can be described as a set of voxels with W = [0, w 1 ] × [0, w 2 ] × [0, w 3 ] ∩ Z 3 for some non-negative integers w 1 , w 2 , w 3 . Each void voxel x = (x 1 , x 2 , x 3 ) ∈ W is said to be a pore voxel of the sheet if the following condition denote the x 2 -value of the highest/lowest solid voxel for given values of x 1 and x 3 .
Box approach. Finally, let us assume the opposite extreme case of perfectly flat, parallel surfaces ( Figure 6C). That is, the solid material is fully contained in a box the volume of which is given by a constant thickness, d , and the nominal sample area. This assumption is similar to the standard acquisition of the so-called apparent thickness; in this standard measurement, probing metal plates touch the paper sheet only at the most protruding features. To determine the position of the smallest box containing the solid material of the entire sample, we can simply proceed with the local information collected in the previous case: Here, m global is the 5% quantile of the full set of encountered local minima m loc (0, 0), m loc (0, 1), . . . , m loc (w 1 , w 3 ). Analogously, M global is the 95% quantile of the full set of encountered local maxima M loc (0, 0), M loc (0, 1), . . . , M loc (w 1 , w 3 ). Note that using quantiles rather than the actual global minimum and maximum allows us to account for spurious outliers due to, e.g. fibres that markedly stick out of the paper sheet. Consequently, all voids encountered in the box are going to be considered as pores of the paper sheet.
Starting out from each of the above definitions for the paper surfaces, we determine the local thickness of the microstructure of samples I and II. In Figure 7, the distribution of the locally encountered thicknesses is shown together with the experimental value (black dashed line) for both samples. The thicknesses derived from the box approach are indicated with the green lines. Note that, due to the local outliers, the paper thickness obtained by the box approach is also a distribution of local thicknesses; local thicknesses larger than d form a narrow tail towards larger thickness values. Strikingly, the box-derived thickness d (green line in Figure 7) is in excellent agreement with the experimental value for sample I. The value d of sample II is smaller than the experimental thickness. We conjecture that this underestimation occurs due to an insufficient sample size. The inspection of the images discarded for the sake of edge correction during the segmentation step revealed the protruding features that corroborate the experimental thickness.
The local thicknesses obtained with the rolling ball and the visibility depth method exhibit broad distributions whose mean values are several tens of microns below the experimental value. This deviation is approximately 40 μm for both samples, i.e. around 30% lower with respect to the apparent thickness. Such a marked discrepancy between apparent and mean local thicknesses is commonly encountered, independent whether the local thicknesses are determined by image analysis of microstructures (Holmstad et al., 2006;Chinga-Carrasco et al., 2008) or local probing techniques (Schultz-Eklund et al., 1992;Dodson et al., 2001;Sung & Keller, 2008). Our analysis suggests that this notable discrepancy essentially originates from an inherent overestimation of sheet thicknesses in caliper-based setups. The latter method inherently records the largest vertical extension of the entire sheet. This is fully in line with the operation of our box approach. When we inspect the distributions of local thicknesses more closely (Figure 7, blue and orange lines), we find that, in fact, thicknesses as large as d are encountered with small, yet nonnegligible contributions at the high-thickness tail of the distributions. The extent of the deviation of the mean local thickness from the maximum thickness is hence expected to enlarge with increasing corrugation of the surfaces.
A further comparison of the rolling ball and the visibility depth method reveals that the latter (i) tends to predict systematically lower thickness values and (ii) encounters markedly lower minimum thicknesses than the former. The reason for this deviation is readily spotted in Figures 6(A) and 6(B). The idealized tip employed in the visibility depth method occasionally penetrates deeply into the paper sheet, yielding a surface with deep pits. At such positions, unrealistically small thicknesses are obtained.

Porosity.
We proceed to determine the porosity, i.e. the ratio between the volume of pore space and the total volume of the considered region of interest. Note that we consider all voxels encountered in between the top and bottom surfaces of the sheet, thus tracking the impact of our surface definition on the porosity prediction. The latter cannot be straightforwardly anticipated, because both the available pore volume and the total volume (via the thickness) will be affected.
The porosity values of samples I and II obtained for each surface definition are displayed in Figure 8(A). Remarkably, the porosity is consistently predicted to be the same for both samples regardless of the applied method. Compared to the rolling ball case, where the porosity is 35.1% for sample I and 38.3% for sample II, the box approach yields large porosity values of 51.5% for sample I and 50.4% for sample II; the boxes incorporate a much larger pore volume due to the lack of discriminating geometric details at the surfaces. In contrast, the visibility depth approach artificially excludes pores in the regions where the pits are formed. Thus, this approach yields markedly lower porosities, i.e. 31.7% for sample I and 35.3% for sample II. On the basis of this comparison, the rolling ball approach emerges as the most reasonable method for determining the paper surface. Thus, the pore space defined by the rolling ball approach will be the basis for the further characterization of the microstructure of samples I and II.
It is intriguing to check whether the similarity in the porosity values is reflected in the mass density of the two samples. The mass density is given by the grammage divided by the thickness. In a first step, we established for each sample that the porosity variation is consistent with the variation in mass density due to differently obtained thickness values. Note that this consistency does not imply that the mass density of the samples can be unambiguously determined. More details regarding the validation of the porosity using mass densities are given in Figure S3. In a second step, the two samples were compared. As is shown in Figure 8(B), the thicker sample II possesses with ca. 110 g/m 2 a profoundly larger grammage than sample I (ca. 70 g/m 2 ); note that the measured grammage is in excellent agreement with the nominal grammage specified  for the paper by the supplier. The values of the corresponding mass density for samples I and II are shown in Figure 8(C). Even though the mass density of sample II surpasses that of sample I, both samples have indeed a comparable density within the margins of experimental uncertainty.

Spherical contact distribution function.
In order to get further information about the geometry of the pore space, we consider the spherical contact distribution function, which is a commonly used characteristic in statistical microstructure analysis (Chiu et al., 2013). The value of the spherical contact distribution function at position r > 0 is defined as the probability that a randomly chosen point of the pore space has at most distance r to the closest voxel of the solid phase. The distance from each pore voxel to the closest voxel of the cellulose material is obtained by computing the Euclidean distance transform according to the algorithm described in Felzenszwalb & Huttenlocher (2012). By virtue of its definition, the spherical contact distribution function accounts for the pore shapes and the pore size distribution (Bezrukov et al., 2002). Hence, it is highly remarkable that the spherical contact distribution functions of the two samples are nearly identical, see Figure 9(A).

Pore network properties
As the two samples are indistinguishable in terms of porosity and spherical contact distribution function, we turn to a more detailed analysis of the network spanned by the pores. Particularly promising is to explore (i) the fraction of pores that can potentially contribute to transport through the paper sheets and (ii) the geometry of correspondingly formed pathways. Basically leaned on the concept of the pore size distribution from simulation of mercury intrusion porosimetry, we geometrically determine the part of the pores that can be reached by an intrusion (from the bottom to the top surface of the sheet) of spheres with a certain radius r (Holzer et al., 2013). In this way, all pores are discarded, which are either not connected to the bottom surface of the paper sheet or which can only be reached by paths going through extremely narrow bottlenecks. That is, pores are discarded which are less relevant for transport from the bottom to the top surface of the paper sheet.
The fractions of pores in samples I and II that can be reached by an intrusion of spheres with different radii is presented in Figure 9(B); two-dimensional visualizations for correspondingly explored pore space are given in Figure S4 in the SI. For r = 0, we obtain those pores, which are connected to the bottom surface of the paper sheet. The next indicated radius 1.5 μm is the size of one voxel, and thus it is the lowest reasonable value for the radius r . For small radii, the predicted accessible pore fraction is nearly identical for both samples. However, when we restrict ourselves to pore radii between 4.5 and 6 μm, less pore space can be reached in sample II compared to sample I.
In terms of paths through the pore network of the samples, the mean geodesic tortuosityτ geod is an illustrative measure for the distribution of encountered paths through accessible pores or, more precisely, their sinuousness with respect to the sheet thickness (Adler, 1992;Axelsson & Svensson, 2010). For an arbitrary pore voxel located at the bottom surface of the paper sheet, its geodesic tortuosity τ geod is defined as the ratio between (i) the length of the shortest pathway from this voxel to the top surface and (ii) the local thickness of the paper sheet at the position of the considered voxel. Then,τ geod is the average of τ geod over all pore voxels at the bottom surface. Values ofτ geod close to unity indicate a dominance of pathways straight through the sheet. For a more detailed description of mean geodesic tortuosity and its quantitative influence on transport properties, we refer to Stenzel et al. (2016).
While there are many coexisting definitions of tortuosity (for an overview, see, e.g. Clennell, 1997), the geodesic tortuosity is of purely geometrical nature. Hence, its value solely refers to the microstructure without being additionally tied to a certain transport process, as would be the case, e.g. for hydraulic tortuosities (Holmstad et al., 2006). To determine the shortest paths from the bottom to the top surface of the sheet, we apply Dijkstra's algorithm on the voxel grid, see, e.g. Thulasiraman & Swamy (1992) and Axelsson & Svensson (2010). Specifically, we determine the shortest paths within the pore space connected to the bottom surface by pores, which can be filled by spheres of radius 1.5 and 3 μm. Furthermore, we consider the extremal intrusion radius r = 0, for which the entire pore space can be considered as accessible pore space. The probability density functions of the encountered τ geod are estimated by a kernel density estimation via diffusion (Botev et al., 2010).
For both samples, the resulting probability density functions of τ geod for intruding spheres of varying radii are shown in Figures 10(A) and (B). Note that the distributions of the geodesic tortuosity do not feature any sharp peak exactly at unity (indicated with a dashed line). Hence, only a vanishingly small fraction of the paths go straight through the pore network. This absence of straight through holes is quite remarkable. Due to the large variation of pore sizes, one would rather expect a superposition of predominantly straight-through pathways through large pores (several μm wide) with longer, more broadly distributed path lengths.
For ideally small intruders (r = 0 in Figs. 10 A and B), both distributions reveal a peak at about 1.2. In other words, in both samples, a large number of shortest path lengths are 20% longer than the thickness of the material. However, the shape of the distributions differs profoundly between both samples. The different shapes are reflected in the values ofτ geod (see Fig. 10 C). At r = 0,τ geod for sample II is around 25% larger than unity, whereas for sample I exceeds unity by around 42%. Remarkably, there are, on average, shorter transport paths in sample II for each of the considered radii. Moreover, the distribution of shortest path lengths is slightly narrower for sample II than for sample I (compare Figs. 10 A and B). It is important to note the perceived differencesτ geod are not necessarily exclusively due to subtle differences in the papers themselves. Rather, the difference observed in sample II might be, at least partially, caused by the additional filling of the cement bag. Potential explanations include aspects such as (i) structural changes in the bag paper, (ii) changes in surface roughness or (iii) the blocking of larger pores due to impurities such as cement grains. Since these aspects may occur in combination, we cannot unambiguously identify the reason for the difference inτ geod from the available amount of data.

Conclusion
We investigated the pore structure of two sack kraft papers, i.e. of the same paper kind, with distinct basis weights and thicknesses. To reconstruct the microstructure, we utilized attenuation μ-CT for a large scanned sample area of 3 × 3 mm 2 with a voxel length edge of 1.5 μm. We employed indicator kriging as a local thresholding method to distinguish between pore and solid regions in the image data. The microscopic pore structure of both samples can only be distinguished when considering a set of multiple descriptors of the pore space and the pore network obtained by methods from spatial statistics and mathematical morphology. The porosities and the spherical contact distribution functions are nearly identical for the two samples; the comparable mass densities of the two samples corroborate this finding. A distinction of the two pore spaces is only possible by means of a network analysis, i.e. the determination of accessible pores and tortuosity. We found that the thickness and porosity values extracted from the microstructure strongly depend on the actual definition of the paper surface. Hence, such values are meaningful only if they are quoted together with the surface definition. For the comparison with experiments, this implies that the microstructure analysis can be adapted to mimic the 'perception' of the paper surfaces by the underlying measuring technique.

Supporting Information
Additional supporting information may be found online in the Supporting Information section at the end of the article. Figure S1. Two cross-sectional images after reconstruction and their corresponding histograms of grey values for sample I and sample II. Grey values are represented with the viridis colour-map. Figure S2. Cut-out of a cross section before (A) and after (B) unbending transformation. Figure S3. Density of sample I. The experimental density based on the caliper is shown in blue. Figure S4. Visualization of the pore space (white), which can be filled by an intrusion of spheres with a certain radius r from the bottom to the top boundary of the sheet. Figure S5. Grey-level correction around impurities. (A) Greyscale image around an impurity (yellow). Figure S6. (A) 3D volume reconstruction of the full volume of sample I with area 2.001 × 2.802 mm 2 and (B) sample II with area 1.666 × 2.812 mm 2 .