2D & 3D particle size analysis of micro-CT images G. van Dalen, M.W. Koster Unilever Research & Development, Advanced Measurement & Data Modelling, Olivier van Noortlaan 120, NL-3133AT Vlaardingen, [email protected]

Introduction Particles are everywhere in our everyday life. Particles can be found in nearly all products. For example bubbles in aerated products 1, water or oil droplets in emulsions, grains in cereals, salt crystals in bouillon cubes and granules in washing powders. Particles may be dispersed in a continuous phase, which in consumer products is mostly air or water. The size of these particles may influence the physical, chemical or physiological behaviour of products or processes. A wide variety of different techniques are used to measure the particle size distribution. These include physical, electrical, imaging and optical techniques. Imaging techniques with image analysis offer the possibility to analyse the size and shape of each individual particle present in a sample. Especially X-ray microtomography (μCT) provides a powerful tool for non-destructive analysis in 3D. Various image analysis methods are used to obtain quantitative information about the size and shape of particles. In this paper 2D and 3D discrete image analysis methods are compared. For each particle a variety of different specific measurements can be made (e.g. location, size and shape)2. In these methods touching particles are separated, so that each individual particle can be identified, counted and measured. In the SkyScan CT-Analyser (CTan) software the separation possibilities are limited. Therefore, it was investigated if the determination of the 3D structure thickness can be used as an indirect method for the estimation of the particle size distribution. These methods were tested on μCT images of aerated emulsions imaged at different storage times3. The bubble size of these samples increases in time (coarsening). This will not only change the mean bubble size, but also the shape of the size distribution curve.

Methods Images were obtained using a SkyScan 1172 desktop µCT system. Details about the imaging acquisition and tomographic reconstruction can be found in a separate paper3. The image analysis methods were tested on identical binary images to exclude influences of different thresholding procedures. For 2D and 3D image analysis of the individual particles, the image analysis toolbox DIPlib (vers. 2.3) from the Delft University of Technology, running under MATlab (vers. 2009a) from MathWorks was used. Also Avizo Fire (vers. 7.0) from the Visualization Sciences Group was used for 3D image analysis of the indivual particles. Avizo Fire combines 3D visualisation with 3D quantitative measurement. For visualisation in 3D space, isosurface rendering was used. Images of time series were aligned in Avizo Fire using affine registration. Before registration a coarse manual alignment was performed using marks on the sample holder wall. The 3D structure thickness was analysed using SkyScan CT-Analyser (CTan) version 1.11.10.0.

2D image analysis A 2D image will not represent the true sizes of objects. Also the number of objects observed in a 2D image does not correspond to the number per unit volume. A 2D image is the intersection of a plane with the three-dimensional space. The observed apparent size distribution is generally different from the true distribution because: - The apparent diameters are as a rule equal or smaller than the true diameters (with a random section the plane may pass through the sphere at the equator, near one of the poles, or anywhere in between) - A random section will contain a relatively greater number of actually large droplets than of small ones (compared to the true distribution), because the former will more frequently be intersected by the section plane. The two causes will, however work in opposite directions, thus having a chance to balance each other. Methods for calculating the true from apparent size distributions are provided by the science of stereology. These methods are statistical in nature, meaning that the size of an individual particle is not determined exactly, but the distribution of many particles can be described. For the relative simple case of samples with disjoint spherical particles, intersected by a plane with zero-thickness a numerical method was developed by Wicksell4 in 1925. This approach is still used as a basic part of more sophisticated methods. The Wicksell transform is based on the distribution of planar intercepts (circle) through a sphere. The more distant the centre of the spheres from the plane the smaller will be the image diameter, in comparison to the actual diameter. The shape of the circle size distribution can be calculated as shown in Figure 6.

Figure 6

Random sectioning of a sphere produces a distribution of circle sizes, which can be calculated from analytical geometry2,5.

For different sizes of spheres the measured distribution will show the superimposition of data from different sizes in proportion to their relative abundance and to their size (as mentioned earlier, it is more likely for the plane to strike a large sphere than a small one). The distribution might be unfolded sequentially since the largest circles could only come from the largest spheres, while all of the spheres may contribute to the smaller circles. For the case where the planar intercept or viewing section is of finite thickness, the Wicksell method has to be extended 6. In many papers unfolding methods are described which were developed for analysing microtome slices by light or electron microscopy. An overview is given by CruzOrive7. These methods are also used for optical sections. For these methods it is assumed that one particle is not overlapped by another particle deeper in the slice. With a low particle

density and thin viewing slice as observed in μCT, this effect will be small. The apparent diameter will be equal to true diameter when the centre of the sphere is within the cross section. The numerical methods for a zero thickness slice4 can easily be programmed by using a matrix of coefficients (i,j), published in the corresponding literature: Si =  i,j . Cj In which : Si = number of spheres in size class i Cj = number of circles in size class j A disadvantage of these methods is that in certain circumstances frequencies of a derived true size distribution become negative for certain droplet diameters. Negative frequencyregions seem to be predominantly found at small diameters and in diameter ranges where the apparent distribution function approaches zero or undergoes a rapid or even unsteady change. These negative frequencies occur because the measured apparent frequencies conflict with one or more theoretically required conditions (e.g. by scattering of the frequencies around 0). The negative frequency effect can be avoided by Wicksell transformation of true log-normal distributions (see section about size distributions). The observed apparent size distribution can now be fitted by a linear combination of the apparent size distributions derived analytically from the log-normal distributions, in order to estimate the parameters of the true size distribution (iterative process from initial to final parameter estimates). This procedure has been implemented in the computer program SGYNWAER using a non-linear least-squares method for estimation of the parameters of a mixture of up to four component distributions. The non-linear least-squares method is Marquardt's8 procedure for least-squares estimation of non-linear parameters. The program has been derived from ROKE, a computer program for decomposition of mixtures of distributions9. The estimated true distribution parameters are then used to represent the expected apparent and the true size distributions and to compute the mean diameters D(p,q) for the true distribution. Input data consists of histogram values, a set of initial estimates for the distribution parameters and the thickness of the section. The name SGYNWAER is a combination of the Dutch words for "Apparent" and "True", written in what seems an old-Dutch spelling10. An example of the calculation of the estimated true distribution from a 2D apparent distribution is shown in Figure 7. 4.0

relative volume, %

3.5

true apparent

3.0 2.5 2.0 1.5 1.0 0.5 0.0 0

1000

2000

3000

4000

equivalent diameter, µm Figure 7

Comparison of the apparent size distribution (2D distribution of circles) and the estimated true distribution (3D distribution of spheres). Bubbles in fresh aerated products are reasonably spherical. However, after storage, also irregular-shaped bubbles are formed. The Wicksell transformation is developed for the distribution of uniform shaped spherical objects and will introduce a bias for other shapes. Other unfolding methods were developed (e.g. for spheroids Wicksel11 and Cruz-Orive7), but

these models are only useful when the three-dimensional shape of the objects is known and is the same for all objects present. Unless the inherent 3D nature of the microstructure is recognised and carefully accounted for the use of 2D image analysis methods can lead to significant errors (bias).

3D image analysis Details about the quantitative analysis of the bubble size distribution are reported in a separate paper12. These methods use the following steps: a) separation of touching particles, b) removing particles touching the lower and upper edge, c) correction for over-segmentation (only in 3D for DIPlib)12, d) measurement of particles and e) generation of a size distribution and calculation of mean particle diameters. For segmentation a watershed transform of the Euclidean distance map was used. In Avizo Fire a build in function was used combining both procedures. The image processing steps in Avizo Fire are shown in Figure 8. After particle size measurement a label image can be generated showing the size classification of each particle (classifying measures with sieves).

Figure 8

3D binary image analysis using Avizo Fire showing images A: used as input, B: after separation of bubbles, C: after labelling (showing different colours for individual separated particles), D-G: after measuring and assigning each bubbles to a predefined size class by using different colours (box size = 11.5mm*11.5mm*21.6mm, pixel size = 24.0 μm, after 4*resampling).

Size distributions and distribution parameters A histogram is one of the simplest ways to display a particle size distribution. It is a particle frequency distribution that shows the percentage of particles found in each size range (bin). The frequency can be plotted on the Y-axis by the number count, surface area, or mass. For many applications a volume-weighted size distribution is preferred. Image analysis of individual particles resulted primarily in a number-weighted size distribution (counting the number of particles/bin) which can be converted easily to an area – or volume weighted distribution. The particle distribution can be presented as a bar or line chart (histogram or curve). An example of a volume-weighted size distribution is given in Figure 9. In this example colours were assigned to the size classes or to a range of size classes. These colours were used to label each particle in the 3D image (Figure 9 & Figure 8). In Figure 10 volume - and number based size distributions are compared. The volume-weighted distribution shows clearly the effect of bubble coarsening whereas the number distribution shows only a small shift to larger sizes. The volume-weighted distribution show 2 different peaks (local maxima). This so called bimodal distribution arises from a mixture of two distributions. For the fresh sample these peaks are not fully separated. After storage two distinct peaks can be observed: relative small bubbles beside large bubbles.

Figure 9

A: 3D visualisation of particles (clipped surface rendering) with B: corresponding volume-weighted particle size distribution. Colours representing a size class or to a range of size classes.

In this paper it is shown that the bubble size distribution of emulsion systems can be adequately described using log-normal distributions. The log-normal distribution plays an important role in particle-size studies13,14,15,16. Massey17 also showed that the bubble size distribution in cake batters under any given set of conditions could be described by lognormal distributions. Jakubczyk & Niranjan18 observed that the lognormal distribution of bubbles in dairy cream changes from unimodal to bimodal during prolonged whipping. The validity of the log-normal distribution has also been derived from theoretical considerations19,20. The process underlying change or growth (forming of bubbles or droplets during processing) is multiplicative rather than additive. The resulting distribution of bubbles in food emulsions can always be described as a sum of lognormal distributions, i.e. in terms of a multimodal distribution. In practice the first two or three lognormal components suffice for the description of the complete distribution. The bubble size D is lognormal distributed if the function f = ln D obeys the normal law (in which f = the frequency of droplets with size D).

Figure 10 Comparison between a volume and a number based size distribution for an aerated sample after preparation (top) and after 1 month storage (bottom). Generally the particle size is reported as an equivalent spherical diameter. This is the diameter of a sphere having the same volume as the feature. The volume is the number of pixels within the particle, which is straightforwardly determined by counting. The equivalent spherical diameter is not only used for particles approaching a spherical shape but also for irregular particles. Only for elongated of fibrous particles the length is used. Usually, an average diameter is used to characterize the size of the distribution. For aerated products often the arithmetic mean diameter (D̅̅ 1,0), the surface-weighted mean diameter (D̅̅ 3,2) and the volume-weighted mean diameter are used (D̅̅ 4,3). The definition and nomenclature of the mean particle diameters D̅̅ p,q have also been reported by Alderliesten21,22,23. A log-normal distribution can further be characterised by the standard deviation () of the logarithm of the bubble diameter (or by exp()). A single mean diameter cannot describe the size distribution of a sample. A better approach is to report a mean diameter and standard deviation or several mean diameters (e.g. D̅̅ 1,0, D̅̅ 3,2 and D̅̅ 4,3).

Time series The image analysis methods were tested on images of aerated products which were acquired over a period of 1-2 months12. For this purpose exactly the same sub sample was imaged. Three samples were selected, showing a different coarsening over time (Figure 11 - Figure 13).

Figure 11 Size labelled μCT images with size distribution of bubbles in aerated sample 1 (storage time in days, box size = 14.3mm*14.3mm*14.4mm, pixel size = 16.0 μm, after 2*resampling).

Figure 12 Size labelled 3D image images with size distribution of bubbles in aerated sample 2 (storage time in days, box size = 14.3mm*14.3mm*14.4mm, pixel size = 16.0 μm, after 2*resampling).

Figure 13 Size labelled 3D image images with size distribution of bubbles in aerated sample 3 (storage time in days, box size = 11.5mm*11.5mm*21.6mm, pixel size = 24.0 μm, after 2*resampling). For all samples large bubbles were formed in time, at the cost of small bubbles. The relative volumes of small bubbles (