A new method to measure pore radius distribution of powders

Abstract

The pore radius distribution plays a significant role in characterizing the porous structure of powder particles, yet quantifying radii spanning 100 nm–10 µm—including both open and closed pores—remains challenging. Here, we propose a new method to measure the pore radius distribution. The method comprises two parts. First, the cross-section radius distribution is measured by an automated routine combining scanning-electron-microscopy (SEM) and deep-learning models. Second, a novel algorithm was developed to convert the cross-section radius distribution into a pore radius distribution. This requires a numerical solution to Wicksell's corpuscle problem, where the new algorithm outperforms the commonly used Saltikov-GCO method. We apply the proposed method to powder samples and compare the result with data from synchrotron X-ray tomography. Our approach provides more information on the distribution and agrees with the result of synchrotron X-ray tomography at a larger scale. As a secondary outcome, the algorithm can also be applied to geology and metallurgy when 3D grain size distribution is calculated from 2D grain size distributions.

---

1 Introduction

The internal porous structure of a powder is a critical determinant of its macroscopic functional properties. The pore radius distribution (PRD), in particular, directly influences key performance characteristics such as rehydration and dissolution kinetics, bulk density, powder flowability, and mechanical strength.^1,2 For instance, in food and pharmaceutical applications, a well-controlled porous network is essential for rapid reconstitution of dehydrated products.³ Consequently, accurate and reliable measurement of the PRD, encompassing both open and closed pores across a wide range of radii (from nanometers to micrometers), is crucial for product development and quality control.

Despite its importance, quantifying the PRD of powders, especially in the 100 nm to 10 µm range, remains a significant challenge. Current characterization methods can be broadly categorized into indirect macroscopic techniques and direct microscopic imaging.

Indirect methods infer pore characteristics from the interaction of the material with a fluid. Mercury Intrusion Porosimetry (MIP) is a widely used technique that calculates pore radii from the pressure required to force mercury into the porous network, based on the Washburn equation.⁴ While it covers a broad measurement range, MIP has notable drawbacks. The high pressures required can compress and damage the sample's delicate structure, leading to inaccurate results.^5,6 Furthermore, the method is unable to characterize closed porosity, poses environmental and health risks, and can misinterpret pore systems with narrow entrances leading to larger cavities (so-called “ink-bottle” pores).^7,8

An alternative is gas physisorption, where the pore size is derived from the volume of an adsorbed gas (typically nitrogen) at various pressures. The Barrett–Joyner–Halenda (BJH) model, which is based on the Kelvin equation, is commonly used for this calculation.⁹ However, as outlined in the IUPAC recommendations, this method and its underlying assumptions are best suited for nanoporous materials and are generally not applicable to pores with radii larger than approximately 50 nm, as capillary condensation effects become negligible.¹⁰ Applying these models outside their valid range can lead to significant errors and misinterpretation of the true pore structure.¹¹ Like MIP, gas physisorption is also limited to characterizing open pores accessible from the particle's exterior.

Direct imaging techniques offer a way to visualize the pore structure without the geometric assumptions inherent in indirect methods. X-ray micro-computed tomography (µCT) can provide a full 3D reconstruction of a particle's internal and external structure, allowing for direct measurement of both open and closed pores.¹² However, a trade-off exists: the spatial resolution of common laboratory-based µCT systems is often limited to several micrometers, while high-resolution synchrotron-based systems, though capable of sub-micron analysis, are constrained by high cost and limited accessibility.^13,14

Scanning electron microscopy (SEM), by contrast, is widely available and offers excellent resolution down to the nanometer scale. However, it presents two significant challenges in interpreting the data. First, classical image processing techniques are unreliable for the present conditions. Classical image processing techniques have been applied to segment images based on well-defined features such as intensity. They have been reported to yield fairly accurate results on SEM images of resin-impregnated and polished geological samples.¹⁵ However, for porous powder samples, resin impregnation can induce internal structural damage and cannot access closed pores. To avoid these issues, the powder is fractured with a blade to expose the internal structure. SEM images are complicated by artifacts such as irregular fracture planes, debris, uneven intensity caused by variable surface direction, and non-circular pore edges. This necessitates the use of more advanced image processing techniques, such as deep-learning models, to robustly identify and measure the pores. Second, the observed circles are cross-sections of the 3D spherical pores, and their radii are necessarily smaller than or equal to the true pore radii. The task of inferring the 3D radius distribution of spheres from the distribution of their 2D cross-section radii is a classic stereological challenge known as Wicksell's corpuscle problem.^16–18 An introduction to numerical methods for solving Wicksell's corpuscle problem can be found in Section 3.1.

This work introduces a novel method to overcome these limitations. We aim to establish an accessible and accurate technique for measuring the pore radius distribution of powders. The approach detailed herein is twofold: (i) developing an automated high-throughput workflow combining SEM with deep-learning-based image analysis to rapidly measure the 2D cross-section radius distribution, and (ii) developing and validating a robust numerical algorithm that accurately solves the inverse Wicksell's corpuscle problem, transforming the measured 2D cross-section data into a true 3D pore radius distribution.

2. Imaging

2.1. Sample preparation for imaging

We prepared samples for SEM imaging as shown in Fig. 1. For each SEM stub, we adhered powders to the top of the stub with conductive tape and removed loose powder by shock and vibration. To expose the internal pore structure upward, we then cut the powder with a razor blade at a low angle relative to the surface and with a uniform and dense spacing. Finally, a plasma gold coating was applied on top to avoid the charging effect during the SEM scan.

Fig. 1 Illustration on cutting particles on SEM stub using razor blade. The angle of the cut makes the cross-section of pores appear to be slightly elliptical when viewed from vertically above. For illustration purpose, the particles and pores are drawn extra large compared with the aluminum stub.

2.2. SEM configuration

A Thermo Fisher Quattro s SEM was used to perform the imaging task. A BSE detector was selected because the intensity is less influenced by the orientation of the sample in comparison with ETD. The electron voltage was set to 10 kV with a spot size of 3.5 and the working distance was 10 mm.

2.3. Automated cross-section radius measurement with SEM

We developed a fully automated routine to measure the cross-section radius distribution of the exposed porous structure with SEM. The routine consisted of 2 stages: Region of interest (ROI) detection at low magnification and cross-section radius detection at high magnification.

2.3.1. ROI detection at low magnification. Each SEM stub was scanned in a non-overlapping 4 × 4 grid of tiles. The magnification was determined visually and kept constant throughout the ROI detection process. Regions of Interest (ROIs) were defined as locations where the porous structure was exposed by cutting. A custom-trained instance segmentation model based on YOLO11-seg¹⁹ was employed to automatically detect ROIs from the scanned images (see SI). Among outputs of the model, only the bounding box outputs were utilized. The coordinates for each ROI were registered as the center of the bounding box, offset by the scanning coordinates. These ROI coordinates were then re-scanned for cross-section detection. Fig. 2 presents an example of the detected ROIs. For visualization purposes, each ROI is marked with a circle with a diameter equal to the average bounding box size.

Fig. 2 SEM image of particles under low magnification. The detected ROIs are marked by green circles. The ROIs are open pore areas resulting from cutting.

2.3.2. Cross-section radius measurement at high magnification. The SEM re-scanned the detected ROI coordinates at a higher fixed magnification, determined by optimizing the trade-off between resolution and coverage. Specifically, the magnification was set to maximize the field of view, subject to the constraint that the smallest pores of interest have size larger than 10 pixels. A custom-trained instance segmentation model based on YOLO V11 seg¹⁹ was then used to predict masks for each pore instance as shown in Fig. 3. To calculate the pore cross-section radii, we assumed the pores are spherical, implying their true 3D cross-sections are circular. However, when a circular cross-section is projected onto an angled 2D plane, it appears as an ellipse. Since the major axis of this ellipse corresponds to the true diameter of the cross-section, the radii were calculated using the minimal enclosing circle of the predicted mask, which effectively captures this maximum dimension. The enclosing circles are visualized on the original image in Fig. 3.

Fig. 3 SEM image of a particle under high magnification. The predicted pore segmentation are masked in blue and the detected pores are marked by green circles on the left hand side of the image. The right hand side of the image is unprocessed for comparison.

Following the automated routine, pore cross-section predictions were subjected to a structured visual inspection. While quantitative performance metrics indicated the model's ability to learn the training dataset (see SI), the limited size of the dataset meant that generalization to new samples was not guaranteed. Therefore, we validated the performance of the model on each new sample at the end of the routine.

To avoid creating full ground-truth annotations for every new sample, we used a human-guided protocol. For each new sample, three random SEM scans were overlaid with model predictions, as shown in Fig. 3. False negative and false positive counts were determined by visual inspection, while the positive count was obtained from the model. Detection results were accepted if (1) the recall was higher than 90%, and (2) the precision was higher than 90%. In cases where detections failed to meet these criteria, the model was retrained using annotated images from the current sample to improve performance. These failure cases are reported in the SI.

3. Stereological transformation

The automated imaging routine provides a 2D cross-section radius distribution. However, the true property of interest is the underlying 3D pore radius distribution. To connect the measured 2D data to the actual 3D structure, we must first find a mathematical transformation that represents the measuring process.

We abstract the measuring process as finding the cross-section radius distribution resulting from uniformly distributed parallel cutting planes intersecting a set of spherical pores with known radius distribution arbitrarily dispersed in a finite 3D space without overlapping, under the following assumptions:

• Pores are spherical

• Distribution of particles on the SEM stub is not biased by the adhesion process

• With respect to the process of cutting particles:

– Section planes resulting from cutting are planar

– Section planes are parallel

– Probability density of section planes is uniform on the SEM stub

• Open pores at the section plane are detected by the model

The explicit expression for the transformation of the measurement process is

where r is the cross-section radius, R is the pore radius in 3D, f_r(r) is the probability density function (PDF) of cross-section radii and f_R(R) is the PDF of pore radius. Eqn (1) has the same form as the original Wicksell corpuscle problem.¹⁶ The Wicksell corpuscle problem asks how to find the inversion of the transformation expressed by eqn (1). Wicksell derived the general solution to this problem:¹⁶

However, in practice, this general solution is rarely applied. This is because the cross-section radius distribution is not differentiable without regularization.²⁰

It is important to note that, unlike Wicksell's original formulation,¹⁶ our derivation does not require a homogeneous spatial distribution of pores (see SI). In our model, the pores are treated as a fixed collection of spheres within the sample volume. By employing cuts that are parallel, random, and uniformly distributed across the entire space, the probability of intersecting any given pore depends solely on its radius and the global sampling density. Consequently, the resulting cross-section radius distribution is independent of the pores' spatial arrangement.

3.1. Numerical methods for the Wicksell corpuscle problem

Here we give a brief review of available numerical methods. Current numerical methods that have been applied to the Wicksell problem can be grouped into two main categories.

The first category consists of distribution-free methods,^21–24 which calculate pore radius distributions without prior assumptions on distribution type. In 1967, Saltikov proposed a numerical technique that iteratively solves for the distribution, conventionally known as the Saltikov method.²¹ In this approach, the continuous pore radius range is discretized into classes. The algorithm employs an iterative procedure: starting from the largest size class, it calculates the contribution of these largest pores to the cross-sectional distribution and subtracts it from the total. This process is repeated sequentially down to the smallest class. In the original Saltikov formulation, the representative pore radius is assumed to be the maximum of the class range. Independently, Goldsmith (1967) and Cruz-Orive (1976) developed an algorithmically identical method, later termed the GCO method, which differs only by defining the representative radius as the midpoint of the class interval.²⁵ While the GCO refinement is reported to yield higher accuracy,²⁶ the mathematical framework of both approaches is indistinguishable. To resolve the nomenclature inconsistencies in the literature and acknowledge both the foundational algorithm and the midpoint refinement, this study classifies them as a single numerical approach referred to as the Saltikov-GCO method.

The Saltikov-GCO method remains popular due to its computational simplicity; however, several critical limitations have been reported. First, the method exhibits a systematic bias, often causing a shift of the distribution median toward larger values.²⁷ This occurs because the algorithm prioritizes fitting the largest radius classes first. Consequently, errors accumulate during the sequential subtraction process, frequently leading to an underestimation of the population of smaller pore radii. Second, practical experience indicates that the method requires typically exceeding 1000 cross-section measurements to yield statistically stable results.²⁷ This high sample size requirement renders the method labor-intensive and potentially expensive in the absence of automated data collection systems.

The second category assumes a distribution with parameters (log-normal, Rayleigh, Weibull, etc) before calculation.²⁸ Using formula like eqn (1), the cross-section radius distribution can be calculated based on a given set of distribution parameters. In this approach, the best fitting distribution parameters can be found based on residual minimization or likelihood maximization. The advantage of parametric methods like this is that they require less measurement data of cross-section radii to achieve similar accuracy provided that the underlying distribution type is known. Evidently, the disadvantage is that it only works if the distribution type is known. For this work, a distribution-free method is preferred. This is because the cost of cross-section radius measurement is made negligible by the automation of the measurement process and the powder pore radius distribution is usually unknown prior to the calculation.

The proposed distribution-free numerical method differs from the previous Saltikov-GCO method in two ways:

1. This method uses a uniform distribution of pore radii within each radius class instead of assigning a single radius to each radius class.

2. This method determines the pore radius distribution through residual minimization instead of iterative solving.

3.1.1. Discretised formulation of the problem. Consider a 3D pore radius distribution f_R(R) approximated as a piece-wise constant function with radius classes defined by a_j

f_R(R) = c_j for a_j−1 ≤ R < a_j, and j = 1, 2,…, M

and the c_j should satisfy normalization condition

Using (1) and the fact that f_R(R) = c_j in each radius class, following simplification can be made,

Based on following integral property

the expression for the forward conversion can be further simplified:

Eqn (3) is well behaved as the piece-wise integration removes the singularity in eqn (1). More importantly, eqn (3) shows that the f_r(r) can be formulated as a linear combination of Φ_i(r)/〈R〉 weighted by probability density c_i for each part of piece-wise function. Similar treatment has been done as preparation for parametric method.²⁸

3.1.2. Least-mean-square-error to estimate the pore size distribution. The residual is defined as a function of class weight c_i

It measures how well does a particular set of c_i can account for the measured cross-section radius distribution. By minimizing the residual, we can find the best estimation of c_i. At first sight, the minimization problem with c_i is not guaranteed to be convex because 〈R〉 is also dependent on c_i. However, by making a change of variable,

we have

Then the minimization problem becomes convex† with respect to the new variable α_i, so that the unique global minimum of Res is where

By collecting the terms, it can be expressed as a linear problem

with

It can be shown that Ψ_ij is a positive definite matrix. Based on the rules of linear algebra, there must therefore exist an inverse Ψ_ij⁻¹ which is also positive definite. However it does not guarantee a positiveness in α_i. Therefore, a non-negative least square solve was used to find the positive solution.²⁹ Based on relation

c_i = 〈R〉α_i

and the normalization condition

The probability density in each radius class, c_i is

3.2. Validation of the numerical method with a simulated experiment

We use a synthetic ground truth and simulated cross section radius measurement to validate the numerical method. The verification follows these steps:

1. Prepare a synthetic pore radius distribution f_R⁰(R) as the ground truth

2. Prepare a cross-section radius distribution f_r¹(r) by a direct simulation using f_R⁰(R)

3. Estimate f_R²(R) by the numerical method based on the f_r¹(r) and compare it with the ground truth f_R⁰(R)

A synthetic f_R⁰(R) is generated by taking 2000 values from a mix of two Rayleigh distributions. This distribution is chosen because a bimodal distribution of pores is common in porous powders and the properties of the Rayleigh distribution guarantees both a positive radius and a centered distribution.

Out of the 2000 randomly generated pores radii, only pores with radii within the range (0.1,10) are kept to reflect the realistic range of values. The range of the pore radius is partitioned into 20 classes with equal width on a logarithmic scale to approximate the probability density of pore radius. The probability density for each class is

where n_i is the number of pores within ith class, and w_i the width of the ith class. The synthetic pore radius distribution f_R⁰(R) is plotted in Fig. 4 in blue.

Fig. 4 The simulated experiment shows that the estimated pore radius distribution (green) matches with the synthetic ground truth distribution (blue). The cross-section radius distribution (orange) is from simulation of cutting spherical pores aligned linearly at equal distance.

The simulated cross section measurement takes following steps:

1. Get all pore radii from the same sampling process above with the same random seed

2. Align all pores centers on one axis with no space in between

3. For each 0.01 along the axis, use a perpendicular plane to intersect the ball and record the cross-section radius

The arrangement of spherical pores and cutting method does not affect the result of simulation compared with spherical pores arranged randomly in 3D space. This is because, firstly, the location of the spherical pore does not affect the probability of it being cut, and, secondly, equally spaced cutting is a good approximation of uniform random cutting. The cross-section radii range is limited to (0.05,10) to reflect the fact that cross-section radius measurement in reality must have a lower limit for detection. The a_i are set to be equally log-spaced within this range with 20 classes. The cross section PDF f_r¹(r) is plotted in orange in Fig. 4.

Finally we apply the method to the cross section radius PDF measured from simulation, The estimated pore radius distribution f_R²(R) is plotted in Fig. 4 in green. It shows that the the numerical method yields a good estimation of the ground truth f_R(R) for the test data.

3.3. Comparison with Saltikov-GCO method

Our method to solve the Wicksell's problem is more accurate than the standard Saltikov-GCO method. To demonstrate this, a numerical experiment with following steps is performed:

1. Randomly sample 1000 values from a lognormal distribution as the pore radius. Values below 0.1 and above 10 are discarded.

2. Use simulation as described previously to generate a set of cross-section radii. Values below 0.05 are discarded.

3. Apply our method to estimate the pore radius distribution. The number of bins is 20 and equally log-spaced.

4. Apply Saltikov-GCO method to the cross-section radii to estimate the pore radius distribution. The number of bins is 20 which are equally log-spaced. The implementation from GrainSizeTools³⁰ is adapted to have equally log-spaced binning.

The result is shown in Fig. 5. The pore radius distribution from our method agrees closely to the ground truth. The pore radius distribution from Saltikov-GCO method shifts to the right compared with the ground truth. This behavior of Saltikov-GCO method is also observed in previous literature.²⁷ Our method performs better than the Saltikov-GCO method, probably because our method solves the best fitting pore radius distribution by minimizing the overall residual instead of fitting the pore radius iteratively leading to an unbiased, uniform distribution of error over the entire distribution.

Fig. 5 The result of this numerical experiment shows that our method (orange) can recover the ground truth pore radius distribution (blue) without significant bias of Saltikov-GCO method (green).

3.4. Corrections to the numerical method to reflect realistic SEM cross-section radius measurement

In reality, the actual measurement process can deviate somewhat from our modeled process. One of the strongest assumptions is that “open pores at the section plane are detected by the model”. In reality, when the size of an exposed pore is very small compared with the field of view, such pores are rarely detected. In this case, we need to look at eqn (4) which implies that all the cross-section radius starting from 0 to infinity are detected and thus used to calculate the residual. We can relax this assumption by defining a new parameter r_min as the smallest detectable pore radius. Effectively, eqn (4) is modified as

Consequently,

where

The effectiveness of this correction was tested based on a simulated experiment with following steps

1. A set of pore radii sampled from a normal distribution N(µ = 4, σ = 1) within the range 0.1 to 10 is generated (arbitrary unit)

2. The measuring process is simulated as before and all the cross-section radii larger than 1 are collected in order to simulate the detection limit of cross-section radius

3. The original method is applied to find the pore radius distribution

4. The corrected method with r_min = 1 is also applied to find the pore radius distribution

The result is shown in Fig. 6. The orange step curve is f_r(r) from simulation and when truncated at r = 1 highlighted with gray dashed line. The pore radius distribution f_R(R) from the original method and corrected method are plotted with green and red respectively. The f_R(R) from the original numerical method has a lower probability density than the ground truth near to r = 1 which means there can be non-negligible underestimation for probability density near the detection limit. On the other hand, the f_R(R) by corrected method with r_min = 1 agrees closely with the ground truth f_R(R).

Fig. 6 Comparison of the estimated pore radius distribution from the original method (green) and corrected method (red) taking detection limit into consideration. The result shows that, for original method, there is a diminishing probability density close to the detection limit, while, for the corrected method, the estimated probability density is well aligned with the synthetic ground truth (blue). The cross-section radius distribution (orange) with cutoff at 1(arbitrary unit) simulates measurement with a detection limit for small cross-section radius.

We have shown that resolution limits of SEM detection can have a non-negligible impact on the result from original method and the proposed correction can solve this issue with a good estimation of r_min. In practice, when we apply the corrected method, the r_min is set to be the smallest detected cross-section radius by default. Further correction on this aspect could be achieved with a calibrated detection rate curve.

4. Application

4.1. Method applied to powder samples with sub-micron pores

We applied our method to 4 samples of porous powder from spray drying. These powder samples came from several factory trials with different process parameters. They contain pores with radii ranging from 0.3–20 µm based on our estimation.

For each sample we prepared three SEM stubs for measurement as described previously. The ROI detections were acquired at 250X and the pore detections at 3500X. The models for ROI detection and pore detection were not retrained for these samples. For each sample, the detected cross-section radii of the three SEM stubs were concatenated. We applied the proposed numerical method with 50 classes to estimate the f_R(R).

Both the f_r(r) from SEM measurement and f_R(R) from the proposed numerical method are shown in Fig. 7. Our observations on the result are that:

Fig. 7 SEM image, Synchrotron X-ray reconstruction, and measured distribution data across 4 different samples. Pore radius distribution are measured by our proposed method (orange) and synchrotron X-ray tomography (green) shown in the third column. Our method reveals more features of the distribution than X-ray tomography. There are three distribution modes 0.5 µm, 0.9 µm and 2 µm present in the four samples, which could indicate different pore forming mechanisms.

1. The modes of f_R(R) from our method are slightly shifted to the right compared with f_r(r). This effect is expected as the radius of a pore is always larger than the radius of the cross-section derived from the pore.

2. The method can resolve modes that are not apparent in f_r(r) measured by SEM. For example, in sample 1, we can see a clear mode at r = 0.9 µm in f_R(R), but this mode is not visible for f_r(r).

3. There are three common modes at r ≈ 0.5 µm, r ≈ 0.9 µm and r ≈ 2 µm for all four samples. If we take into account that the material matrix and processing is similar for all samples, these modes might be signatures for three different pore formation mechanisms.

4.2. Comparison with synchrotron X-ray tomography

The best resolving reference method currently available is synchrotron X-ray 3D tomography. The reconstructed tomography slices are shown in Fig. 7. The details about the tomography and post-processing can be found in SI.

We plotted the f_R(R) from synchrotron X-ray 3D tomography against our result in Fig. 7. We see that the f_R(R) measured by 3D tomography only has one mode due to limitation of the resolving power of 3D tomography. The mode lies between the modes measured by our method indicates that it is likely an overlap of two widened neighboring modes due to uncertainty of radius measurement of 3D tomography.

It is worth mentioning that while higher tomographic resolutions are technically possible, they require high-intensity beams that increase the risk of sample damage, such as melting. More importantly, achieving such high resolution restricts the field of view, which limits the number of particles that can be analyzed in a single scan. The synchrotron data used in our comparison was limited to a small field of view (250 µm × 250 µm × 250 µm), comparable to the average size of a particle shown in Fig. 7. Assuming spherical particles with a radius of 50 µm and a packing density of 0.74, the number of powder particles measured is estimated to be 22. Consequently, it is likely that this small sample size introduces bias into the resulting pore radius distribution, even though the total number of pores measured ranges from 6515 to 22 603.

5. Discussion

The new method presented in this work provides a practical and accurate way to characterize the pore radius distribution in powders, particularly in the 100 nm to 10 µm range. Our approach combines high-throughput SEM imaging with a robust numerical algorithm that overcomes the limitations of previous stereological methods.

5.1. Decoupled spatial resolution and statistical robustness

A primary advantage of our method is its high resolution and good statistical robustness when compared to synchrotron X-ray tomography. Synchrotron tomography, while a powerful non-destructive reference, faces a critical trade-off between resolution and field of view. In contrast, our SEM-based method decouples this trade-off. The resolution is set by the SEM imaging technique (magnification and detector), while the statistical significance is independently determined by the number of images and samples analyzed, allowing for both high resolution and a large, representative sample size. Therefore, our method can have high spatial resolution without sacrificing the statistical robustness.

5.2. Unique position in porous powder characterization

Our method complements existing porosimetry techniques with its unique contribution. Its distinct advantage is the ability to characterize the total pore radius distribution—including both open and closed pores—at a sub-micron resolution. It addresses the limitations of the other two primary methods X-ray tomography and MIP which are limited resolution and measurement induced structural damage.

However, our method cannot provide information about closed versus open pore radius distributions since determining whether a particular pore is open or closed requires taking the entire internal particle structure into account. Synchrotron X-ray tomography, even though limited by resolution, gives a 3D structure of the powder particle and thus allows easy classification of pores using technique such as 3D flooding. MIP also measures the open pore size distribution if closed pores are not collapsed during the measurement.

5.3. Limitations and potential improvements

This method relies on several key assumptions and operational steps, each of which presents potential limitations.

First, sample preparation could be biased. Particle segregation during storage or handling can introduce sampling bias. We sought to mitigate this by thoroughly mixing the powder and preparing three separate SEM stubs per sample, which showed high consistency. However, we can't prove that the particles sticking to the SEM stub are an unbiased representation of the whole powder sample. One potential improvement could be to take a small subsample of the powder after thoroughly mixing the whole sample. Then use multiple SEM stubs to stick all the particles in the subsample. In this way, we could ensure that the powder particles sticking to the SEM stubs is an unbiased subsample of the whole powder sample.

Second, the ROI detection model could potentially introduce a systematic bias. To assess this, we categorized false negatives into two types based on their root causes. The first and most prevalent type stems from sample preparation and SEM scanning limitations, where ROIs exhibit low contrast, are situated near image boundaries, or possess a nearly vertical orientation (see SI). Since these factors are independent of the local pore radius distribution, they do not introduce a systematic bias; their primary impact is a reduction in sample size, which can be mitigated by increasing the number of replications. The second, minority type involves false negatives that correlate with the pore structure itself–for instance, particles with specific pore distributions that are fragile and prone to fragmentation during cutting, resulting in fragmented pieces that are less likely to be detected as ROIs. While this type has the potential to induce systematic bias by excluding specific pore structure, its quantitative effect was not measured in this study.

Third, the cross-section detection model could introduce a systematic detection bias. For instance, the model might be more proficient at identifying high-contrast pores versus low-contrast ones, or more regular circles versus slightly elliptical cross-sections, independent of their actual size. This could skew the initial 2D cross-section radius distribution before the algorithm is applied. In this work, the model's performance was only evaluated qualitatively by visually inspecting the annotated images. While this was sufficient for this proof-of-concept, for more rigorous applications, a quantitative study of the model's bias would be essential. Additionally, with standard sample targets, the cross-section detection rate should be correlated with the radius. Both combined should give a more quantitative estimation of error in the cross-section radius distribution.

Finally, the stereological conversion itself assumes all pores are spherical. This is a reasonable assumption for the spray-dried powders studied, as surface tension strongly favors spherical shapes, especially for smaller pores. However, for larger pores, ellipsoidal shapes are observed and they are caused by internal compression during fast solidification. Potential improvements could be applied to the stereological model itself to extend beyond the current assumption of spherical pores to account for ellipsoidal shapes.

5.4. Broader applications for the numerical algorithm

A significant outcome of this work is the development and validation of a robust numerical solution to the Wicksell corpuscle problem. The simulation in Fig. 7 confirms its high accuracy compared to the commonly used Saltikov-GCO method. Unlike iterative approaches like Saltikov-GCO, our non-parametric algorithm successfully avoids characteristic artifacts such as rightward shifts and underestimation of smaller populations. By solving the entire distribution at once via residual minimization, our algorithm provides a more accurate reconstruction of the ground truth 3D distribution.

While we have demonstrated its use for powder porosity, the algorithm itself is fundamentally a stereological tool and is not limited to this application. As noted in the introduction, this accurate, non-parametric algorithm can be directly applied in any field that needs to infer 3D size distributions from 2D cross-sections. For example, in Geology, it could be used to determine the 3D grain size distribution of minerals or the size of vesicles (bubbles) in volcanic rock by analyzing 2D petrographic thin sections. Similarly, in Metallurgy, it offers a more precise way to quantify the 3D grain size from 2D polished and etched metallographic samples; this is critical, as the 3D grain distribution directly influences the mechanical properties, such as strength and ductility, of the final material. In both fields, our method presents a clear improvement in accuracy and robustness over traditional stereological analyses.

6. Conclusion

In this work, we have presented a novel, two-part method for measuring the pore radius distribution in powders, addressing the challenge of characterizing porous structures in the 100 nm to 10 µm range. Our approach combines an automated routine using Scanning Electron Microscopy (SEM) and deep-learning models for high-throughput measurement of 2D cross-section radii, with a new numerical algorithm to accurately convert this 2D data into a 3D pore radius distribution.

The validation of our method demonstrates its effectiveness and accuracy. The numerical algorithm for solving the Wicksell corpuscle problem proves to be more robust than the commonly used Saltikov-GCO method, avoiding the characteristic rightward shift in the distribution and aligning more closely with ground truth data in simulations. When applied to real powder samples, our method yields results that are consistent with data from synchrotron X-ray tomography at a larger scale. Crucially, our SEM-based technique provides enhanced resolution, revealing distinct modes within the pore distribution that are not resolved by synchrotron tomography.

The significance of this method lies in its ability to provide detailed and accurate pore structure characterization without requiring access to limited and costly synchrotron facilities. Furthermore, the inverse transformation algorithm developed here is not limited to powder analysis; it has broader applications in fields such as geology and metallurgy, where the challenge of inferring 3D grain or particle size distributions from 2D cross-sections is frequently encountered.

7. List of symbols and notation

Table 1 lists all the units, abbreviations and mathematical symbols that are frequently used in this work.

Table 1 List of symbols and abbreviations

Units
µm	micrometer
nm	nanometer
Abbreviations
SEM	Scanning electron microscopy
ETD	Everhart–Thornley detector
BSE	Backscattered electron
ROI	Region of interest
PDF	Probability density function
Mathematical symbols
R	Radius of a sphere/pore
r	Radius of the cross-section of a sphere
fR(R)	PDF of R over all spheres in the system
fr(r)	PDF of r over all cutting events

---

Author contributions

Conceptualization: Y. L., A. B., M. A., I. G., S. T., data: Y. L., J. B., formal analysis: Y. C., A. B., investigation: Y. C., J. B., A. B., methodology Y. L., A. B, J. B., supervision: A. B., M. A, writing – original draft: Y. L., A. B, writing – review & editing: Y. C., A. B., J. B, H. J. L.

Conflicts of interest

There are no conflicts of interest to declare.

Data availability

The data supporting this article have been included as part of the supplementary information (SI). The code implementing the numerical method to convert 2D radius distribution to 3D radius distribution can be found in https://github.com/YueCLin/2D-3D-PSD-too. Supplementary information is available. See DOI: https://doi.org/10.1039/d5sm01199j.

Acknowledgements

The authors would like to thank Nestlé Research for providing access to Scanning Electron Microscopy (SEM) facilities and the computational resources used for model training. We also thank the Nestlé Product Technology Centre (NPTC) for providing the powder samples. Finally, we acknowledge the financial support and assistance provided for the Synchrotron X-ray measurements.

References

1. P. A. M. Le-Bail, D. A. D. C. Kautzmann and P. A. M. Le-Bail, J. Food Eng., 2019, 248, 49–65.
2. A. B. Yu, Ind. Eng. Chem. Res., 2003, 42, 1249–1258.
3. K. P. Schuck, S. L. Jeantet and A. B. Schuck, Dairy Sci. Technol., 2012, 92, 347–366.
4. E. W. Washburn, Phys. Rev., 1921, 17, 273.
5. L. D. G. G. Giesche, Part. Part. Syst. Charact., 2006, 23, 9–19.
6. A. B. Abell, K. L. Willis and D. A. Lange, J. Colloid Interface Sci., 1999, 211, 39–44.
7. P. A. Webb and C. Orr, Analytical Methods in Fine Particle Technology, Micromeritics Instrument Corporation, 1997.
8. J. Kaufmann, R. Loser and A. Leemann, J. Colloid Interface Sci., 2009, 336, 730–737.
9. E. P. Barrett, L. G. Joyner and P. P. Halenda, J. Am. Chem. Soc., 1951, 73, 373–380.
10. M. Thommes, K. Kaneko, A. V. Neimark, J. P. Olivier, F. Rodriguez-Reinoso, J. Rouquerol and K. S. Sing, Pure Appl. Chem., 2015, 87, 1051–1069.
11. J. C. Groen, L. A. Peffer and J. Pérez-Ramrez, Microporous Mesoporous Mater., 2003, 60, 1–17.
12. E. Maire and P. J. Withers, Int. Mater. Rev., 2014, 59(1) DOI:10.1179/1743280413Y.0000000023.
13. D. Wildenschild and A. P. Sheppard, Adv. Water Resour., 2013, 51, 217–246.
14. L. B. E. C. C. H. C. C. T. S. Fezzaa, K. Wang and W.-K. Lee, JOM, 2008, 60, 24–29.
15. H. Safari, B. J. Balcom and A. Afrough, Comput. Geosci., 2021, 156, 104895.
16. Wicksell, The corpuscle problem. a mathematical study of a biometric problem, unknown technical report, 1925.
17. J. C. Russ, Practical Stereology, Plenum Press, 1986.
18. J. Ohser and F. Mücklich, Statistical Analysis of Microstructures in Materials Science, Wiley, 2000.
19. G. Jocher and J. Qiu, Ultralytics YOLO11, 2024, https://github.com/ultralytics/ultralytics.
20. D. Depriester and R. Kubler, J. Struct. Geol., 2021, 151, 104418.
21. S. A. Saltikov, The determination of the size distribution of particles in an opaque material from a measurement of the size distribution of their sections, Erevan polytechnical institute technical report, 1967.
22. M. Simmons, P. Langston and A. Burbidge, Powder Technol., 1999, 102, 75–83.
23. P. L. Goldsmith, Br. J. Appl. Phys., 1967, 18, 813.
24. L.-M. C. Orive, J. Microsc., 1976, 107, 235–253.
25. T. Ueda, Powder Technol., 2022, 404, 117462.
26. S. Weibel, Stereological Methods: Practical methods for biological morphometry, 1979.
27. M. A. Lopez-Sanchez and S. Llana-Fúnez, J. Struct. Geol., 2016, 93, 149–161.
28. D. Depriesterb and R. Kubler, Image Anal. Stereol., 2019, 38, 213–226.
29. C. L. Lawson and R. J. Hanson, Solving Least Squares Problems, Classics in Applied Mathematics, 1995 DOI:10.1137/1.9781611971217.
30. M. A. Lopez-Sanchez, J. Open Source Softw., 2018, 3, 863.

---

Footnote

† A convex function in this context is one for which all of the second order partial derivatives are positive. In the case of eqn (5) these second order partial derivatives are

---