Kalisphera: an analytical tool to reproduce the partial volume effect of spheres imaged in 3D

In experimental mechanics, where 3D imaging is having a profound effect, spheres are commonly adopted for their simplicity and for the ease of their modeling. In this contribution we develop an analytical tool, ‘kalisphera’, to produce 3D raster images of spheres including their partial volume effect. This allows us to evaluate the metrological performance of existing image-based measurement techniques (knowing a priori the ground truth). An advanced application of ‘kalisphera’ is developed here to identify and accurately characterize spheres in real 3D x-ray tomography images with the objective of improving trinarization and contact detection. The effect of the common experimental imperfections is assessed and the overall performance of the tool tested on real images.


Introduction
Full-field 3D imaging is causing a revolution in material science by allowing microstructures to be imaged directly. In the case of non-destructive techniques, the evolution of the microstructure can also be followed. Image-processing tools are typically used or specifically developed to make quantitative measurements from such images. The analytical calculation of the metrological properties of image-processing operations is difficult, and a given image-based measurement tool typically strings together a number of different operations. The metrological behavior of such tools is usually characterized by applying it to a number of known, representative, artificial scenarios.
In the world of granular physics or geomechanics, spheres are pervasively used as a simplified microstructure. Insofar as experiments are concerned, glass and steel beads are often adopted [1][2][3]. In numerical simulations, sphere-like particles are deep within the huge emergence of the 'Discrete Element Method': particle-based simulations based on molecular dynamics, where the fundamental unit is a sphere (see [4] and thousands of citing articles). Spheres are also favored representations of particles in analytical developments, where their simplicity allows for closed form solutions for a number of phenomena and grain-assembly properties (e.g. capillary bridges [5], contact indentation [6], solid volume fraction [7], packing stiffness [8]). It is important not to forget that outside the mechanics of particulate media, the microstructure defined by a series of spheres is typically a dual of that defined by a foam.
A full-field 3D density image (e.g. one coming from x-ray tomography) of a sphere suspended in air will in an ideal case measure the sphere's density in the voxels inside the sphere and the air's density away from the sphere. For the voxels that fall on the interface between sphere and air, some intermediate value-proportional to the ratio of the partial occupancy of each phase within the voxel-will be measured. This effect is typically known as the 'partial volume effect' and is illustrated in the top right of figure 1. The accurate reproduction of this boundary-effect is essential when generating synthetic spheres, especially when one is interested in the characteristics of the boundary of the particle (inter-particle contacts, interface with another phase, etc.) The calculation of the 'grayscale' value of partially-occupied voxels on the boundary of a sphere is not trivial, and can intuitively be handled by subdividing a voxel and calculating whether each part of this subdivided voxel is inside or outside the sphere. This 'brute-force' approach may be acceptable in some preliminary studies, but where good metrological quality is required (such as a precise center of mass, or a precise total mass), the number of voxel subdivisions required amounts to an enormous computational cost.
In part 1 of this paper we present an analytical (and therefore fast and accurate) solution to this problem, directly giving the grayscale value of each voxel-this allows 'perfect' spheres to be generated, providing a ground truth.
In part 2 of this paper, we very briefly discuss some examples of current uses of this tool for metrological measurements.
In part 3 we present in detail a technically advanced application, where the synthetic spheres are guided onto a real, noisy 3D image by 3D template matching in order to characterize the (imperfect) spheres therein. This also allows them to be 'removed' (i.e. subtracted) from the image, which in turn allows a better identification of the remaining phases (e.g. water, cement) surrounding them.
Both the codes used to generate the spheres and to identify and describe spheres in 3D images are made available at https://forge.3sr-grenoble.fr/repos/kalisphera/. At the time of writing a Python version is available and a C++ one is under construction; their use is illustrated through working examples.

Kalisphera, an analytical tool to generate correct 3D images of spheres
A three-dimensional raster image of a sphere is made up of voxels which can be put into one of three categories: voxels wholly inside the sphere (typically given a value of one), voxels wholly outside the sphere (typically given a value of zero) and voxels on the boundary of the two (with a value between one and zero), as shown in the voxel, top right in figure 1. The value of voxels on the boundary is a direct ratio of the (unit) volume of the voxel and the volume of the sphere within the same voxel, and consequently this phenomenon is named 'partial volume effect' and is found in all x-ray tomographic images of heterogeneous materials.
A simple approach to determine the value of these voxels is to subdivide each voxel into, for example, a thousand subvoxels ( × × 10 10 10) and determine for each subvoxel, whether it is inside or outside the sphere; summing the number of subvoxels within the sphere and dividing by the number of subvoxels gives the voxel's value. The accuracy of this measurement is obviously directly related to the number of subdivisions and, for reasonable accuracy, the number of computations for each border-voxel starts to become very large.
The proposed approach circumvents this problem by analytically solving the intersection of a sphere with a voxel (i.e. a cube). For the scope of this study, we are interested in spheres that can circumscribe a voxel (i.e. ⩾ ( ) R L 3 /2, where L is the size of the voxel and R the radius of the sphere). A classification of the different possible intersections, based on the number of cube corners within the sphere, is presented in figure 2. For each case there corresponds a different formulation of the integral defining the volume of intersection, as shown in table 1, where the analytic solution functions S i 's are defined in table 2. The choice of the orientation of the axes in figure 2 and in the corresponding table, is entirely arbitrary and any partial volume effect voxel can be reoriented to have its axes coincide with the ones in the figure.
To build then a correct 3D image of a sphere it is sufficient to classify each voxel in one of the classes and determine its value. A script built for this very purpose by the authors is publicly accessible at https://forge.3sr-grenoble.fr/repos/kalisphera/. The software takes as an input the sphere size and position within a box as well as the Gaussian blur (covered later).
With this approach the speed to produce a sphere increases by several orders of magnitude. For example, a single sphere  Definitions and solutions of the intersections between a sphere and an inscribable cube as per classification in figure 2.

Case
Problem Solution , , , , , , , : , , , , , , , , , , , : , , , , , , , , : , , , , , , , , : , , , , , , , , , , , , : with a radius of 12 pixels, like the ones used in section 2, can be produced with the available Python version of kalisphera in about one second on a single 2.2 Ghz CPU. This is to be compared to the four hours needed with the brute force approach, by subdividing each partial volume voxel into 50 × 50 × 50 subvoxels, which still yields an error of 1% on the volume. An image containing multiple spheres can be generated adding images of individual spheres with identical box size, taking the background (void) grayvalue to be equal to zero. Both the background and the sphere grayvalues can then be set to the desired values.

Uses for metrology
The kalisphera analytical tool described above generates a single 3D image of a sphere as it would be captured by a Note: The axes and points letters adopted above are respectful of the definitions in such figure. The fundamental solutions, marked by the symbol ' a ', are listed in table 2. Case Solution  perfect imaging system. As stated in the introduction, this can offer a ground truth for many different types of image-based measurements of microstructure. For example, at the time of writing, this tool is being used to make multiple realizations of 3D spheres in order to investigate two metrological properties of the measurement of contacts between particles (with specific and direct applications to measurements of fabric in granular materials such as sand): (1) To investigate the sensitivity of contact detection techniques when a physically-correct threshold (i.e. one that selects the physically correct volume) is chosen. (2) To test the precision of a 3D level-sets surface convergence to the boundary of a grain in the presence of contacting grains, imaging noise, and blur.
These specific metrological applications will not be discussed in more depth for sake of brevity, and will be the object of subsequent publications. In general, the idea of releasing this tool publicly is to encourage widespread use of correct spheres for testing image-based measurements (of microstructure etc).
The next section will develop an advanced example of the potential applications of this code, employing it to make measurements, based on the hypothesis that a scanned image contains spheres.

Identification and subtraction of spheres in 3D images
This section introduces a kalisphera-based technique for sphere-matching in real 3D images of spheres. The original objective of this technique was to accurately characterize spheres in 3D images of three-phase media (in this case grains-cement-air but also grain-water-air is an immediate possibility) to subtract the sphere from the image including its partial volume effect to be left with an easy-to-quantify map of cement (or water) and air-particularly close to the contact, where cementation (or the occurrence of a liquid bridge) is especially important for the mechanics of granular media. This section will first introduce a general framework for sphere matching, then discuss its application to real images.

Ideal sphere-matching
The problem of the identification and geometrical description of spheres in an image, given their analytical description, can be classified as a template matching problem [9], although other techniques such as the Hough transform [10] are also applicable. A commonly adopted approach for this kind of problem is to maximize a measure of the similarity of a reference image (typically a subset) to an image-template which can be varied to match. The cross correlation between the two images provides a simple yet robust measure of similarity. This measure is often normalized by the mean intensity of the figure to account for the possibility of non-homogeneous brightness, to obtain a measure named normalized cross correlation coefficient (NCC). For two 3D images it is typically defined as: where ( ) I x y z , , 1 and ( ) I x y z , , 2 are the grayscale values of the voxel of coordinates ( ) x y z , , for the image subset and the template respectively, while u, v, z are the integer displacements of the template sphere in the x,y and z directions respectively. Written in this form, the NCC value is 1 for two identical images.
Different techniques are typically adopted to interpolate this coefficient to obtain sub-pixel (-voxel) accuracy, but in our case it is possible to simply generate a new template with subpixel displacement. It also straightforward to vary the radius of the sphere. For the problem at hand the potential can therefore be rewritten as: where S is the raster image of the sphere, X Y Z , , c c c are the coordinates of its center, and R is its radius. Figure 3(a) shows the evolution of − 1 NCC (so that a value of zero means a perfect match), using as a reference image a single kalisphera sphere with a radius of 12 pixels, and varying the radius and position of the matching sphere (the − − x y z direction displacements generate identical curves). When using the same position and radius as the reference image, it is encouraging that the calculated − 1 NCC is exactly zero (perfect match). Both the variation in displacement and in radius generate smooth, monotonous potentials down to the desired value. The effect on the NCC of the positional error is symmetrical (as could be expected), unlike the radius where the intersecting volumes are different for an undersized or oversized matching sphere. Figure 3(a) shows clearly that even starting from template with a poor guess of sphere center and radius, it is possible (by iteratively optimizing the − 1 NCC) to converge to the exact match. While the function is not convex, it does not present other stable minima than the correct position and size, as it can be observed from the plot of the full landscape in figure 3(b). Furthermore, the smoothness and increasing slope towards the exact match make this a function particularly well suited for gradient-sensitive minimization algorithms. For this work a Broyden-Fletcher-Goldfarb-Shanno algorithm [11] was used (more specifically the one implemented in scipy [12]). While this optimization algorithm is not ideal for non-convex functions, the minimized function is convex relatively close to the global minimum. Alternative optimization algorithms are pre-implemented in the code provided, at the discretion of the user. Typical convergence curves are shown in figure 4. Here, starting from increasingly inaccurate guesses, the optimization procedure allows the estimation of the center of the sphere to converge with the actual value. We remark that the ordinate axis in the figure is logarithmic and that the lower limit of − 10 2 (pixels) coincides with the initial optimization step adopted. This was chosen as a trade off between the desired accuracy and the speed of convergence. For accuracies significantly below one 100th of a pixel, several more steps are required, possibly because of the noise induced by the numerical truncation of the analytical solutions.
It is evident from this image that, while the proposed technique allows the detection and characterization of spheres starting from virtually no prior knowledge of the image (in the uppermost case, the reference sphere and the template one do not touch), the process is considerably sped up by even rough guesses of these features. Unlike the accurate measures the proposed technique aims to provide, rough estimates are usually easy to obtain. For example it is possible to segment the image and determine the center of mass and volume of the objects therein. In some works (e.g. [13]) template matching is performed in the transform domain, although in this context its computational cost is comparably low.

Matching spheres in real images
To apply the sphere-matching technique outlined above to real images of particles from x-ray tomography, the tool must be able to converge in the presence of the complications of real images: the three following crucial factors are explored: (1) Presence of other, touching, spheres (2) Image quality (random noise, and blur from the imaging system's point spread function) (3) Imperfections in the geometry of the real spheres (asphericity)

Presence of other spheres.
As mentioned above, the aim of this application is the characterization of spheres in a real image (e.g. for contact detection or to remove spheres from an image to facilitate segmentation in multi-phase materials). The obvious interest is to apply this to an image of a microstructure composed of many particles. The presence of other particles in contact with the particle of interest (the one on which it is trying to converge) will disturb the relatively perfect matching profile shown in figure 3(a). This is shown in figure 5, where the − 1 NCC profile for one template sphere ( = R 12 px) being moved from left to right is reported for different configurations of contacting spheres.
The effect of contacts above and below the scanning direction is minimal (solid black and dashed azure lines), since the only disturbance is in the pixels shared at the contact. However, this disturbance means that − 1 NCC never reaches zero. The effect of other particles on the path of the scanning direction (red dotted line-blue dash-dotted line) is significant: the occurrence of two additional local minima in front of and behind the central particle corresponding to the centers of the surrounding particles. The convergence on the desired center occurs if the initial guess is less than a radius away while the accuracy of the optimization is essentially unaffected.
The match profile for the variation of the template's radius in figure 6 (solid black line) shows that a false, local minimum is provoked when the radius of the template is large enough  to contain part of the surrounding particles 1 . This means that achieving convergence on a particle of interest when trying to optimize the variable of the radius is feasible only where the initial value of the radius is close to the real one 2 .
In order to overcome this difficulty, the comparison of the template and the subset of the reference image are changed in order to take the boundary of the grains into account. There are a number of means of achieving this through edge detection operators such as Gaussian, Laplacian, and Sobel functions, and in this case we adopt a variance filter. This means that in the computation of the NCC instead of comparing the two original images, compares two variance filtered images 3 .
This yields a single minimum of 1-NCC, as shown in the dashed red line in figure 6. The use of the variance also has the advantage of accelerating the convergence process while maintaining the accuracy of the measure unaltered. Nonetheless, the adoption of this operator increases the sensitivity of the technique to noise, as detailed in the next section.

Image quality.
Several categories of defects and artefacts can occur in x-ray images (e.g. beam hardening, ring artefacts [14]), and a vast and expanding literature aims at their mitigation. In this section we are interested in the effect on the technique of two of the most persistent ones: blurring and random noise, although all the tests on real data include ring artefacts and beam hardening which have been compensated for as much as possible in pre-processing.
Blurring. Blurring can be interpreted as the result of the convolution of the ideal image with the point spread function of the optical system. While methods to deconvolve the image and reduce the blur have been proposed, they often require accounting for numerous features and settings of the imaging system. Instead, in this technique we choose to accept the image as-is and to reproduce the inherent blur on the matching template. The point spread function is approximated by a Gaussian blur operator applied on the generated spheres. The standard deviation of the Gaussian blur is another parameter (along with the sphere position and radius) that can be optimized, although in practice it has been found (as expected) to be constant for a 3D x-ray tomography image.
Random noise. To assess the effect of noise on the accuracy of the technique, random noise with Gaussian distribution of the intensity was added to each voxel of the image. To describe the level of noise in images a commonly adopted measure is the signal to noise ratio, defined as:  Radial intensity profiles for a sphere surrounded by 6 spheres. The potential has a spurious minimum when the radius incorporates the surrounding spheres. The application of an edge detecting filter (variance) eliminates the spurious minimum while maintaining the accuracy of the approach unaltered. 1 The NCC in figure 6 is computed on a box of fixed size and large enough to contain all the spheres, in order to isolate the effect of the sphere radius and remove that of the box size from the figure. This justifies the difference with the values in figures 3 and 5. Using a box with increasing size yields a qualitatively analogous result, and the values of the minima get closer together. 2 For sphere packings with a narrow range of radii (low dispersity or monodisperse packings) this may not be an insurmountable problem as it is sufficient to restrict the range of the radius.
where μ i is the mean gray value of the image and σ n the standard deviation of the background. In experimental mechanics images with < SNR 10 are rarely accepted. As shown in figure 7(a) for SNR levels typical in x-ray images, the noise does not affect the accuracy of the technique and only marginally reduces the convergence speed when the original image is processed. Unsurprisingly, edge-sensitive operators as the variance filter are more strongly affected by noise, as shown in figure 7(b). At ⩽ SNR 5 this filter can generate local minima far from the actual value, rendering the optimization process more difficult. The use of the variance filter is therefore a trade off between the removal of spurious minima in polydisperse samples (and accelerated convergence) and the increased sensitivity to noise. Its adoption should be functional to the type of sample (i.e. level of polydispersity and accuracy of the initial guess) and to the SNR of the images at hand. An hybrid of the two techniques can also be imagined, taking advantage of the variance filter to approach the solution, followed by further steps on the original image.

Non-sphericity of particles.
While techniques to produce spheres with deviations from the theoretical radius smaller than − 10 7 of R are known [15], most commercial grade spheres have a range of asphericity of ± [0.07-2.5] μm. Assessing the implications of these imperfections on the accuracy of the predictions of the technique is of course pivotal to any serious metrological study. To simulate the imperfections of the spheres, ellipsoids of revolution were numerically generated with a technique analogous to that described as 'brute force' in section 1, subdividing each partial volume voxel into 50 3 subvoxels. The asphericity can be defined as = − As c a where c and a are the major and minor radii, respectively, as sketched in figure 8(a).
The matching profile for the template's position along the major radius in figure 8 shows an area of essentially flat NCC in the proximity of the center of the sphere. This hinders the convergence, often leading to estimation errors. This area of uncertainty is of the same size as the asphericity and both in the original and variance filtered image, as shown in figure 8(b). In both cases, the radial intensity profile can be shown to have a unique minimum around ( + + ) a b c /3.

Validation on real images
In this section the performance of the kalisphera spherematching tool developed above is evaluated on two real images  of a sphere packing. The two 3D images are x-ray tomographies acquired on the laboratory scanner in Laboratoire 3SR (Grenoble) of a packing composed of a few thousand sapphire spheres (of three different radii-400, 300 and 250 μm) purchased from Sandoz Fils S.A., with a tolerance on the diameter of ±1.25 μm, and imaged at a pixel size of 8.7 μm/pixel. This microstructure is imaged in two different configurations differing by just a rigid-body displacement of the microstructure. This means that the displacement of each grain between the images is identical. Initial guesses are obtained by thresholding, separating (with a 3D watershed) and tracking the grains between images with ID-Track [16]. From these initial guesses spheres are matched using, in this case, six iterations of the BFGS algorithm (generating kalisphera spheres with the measured average pore value in the image and the sphere's specific gray value 4 ) and using these images with representative grayvalues to calculate the variance. Each iteration step takes roughly twice the time needed to generate a single sphere (detailed in section 1) multiplied by the number of axes on which the optimization is run (i.e. x, y, z, R), plus one. As stated above, the Gaussian blur applied to the kalisphera spheres is fixed, in this case to a standard deviation of 0.81 pixels, based on some initial iterations. In the assumed, ideal case of no sphere rearrangement between the two imaged steps, the measured distribution of displacement is a Dirac delta function. Figure 9 shows, in solid black, the displacement of the centers of the spheres obtained with the kalisphera spherematching tool, which shows a tight distribution around a displacement of 59.55 pixels (note that the entire extent of the x-axis in figure 9 is 1.5 pixels). Practically all the measured displacements fall within 0.3 pixels of the mean: the Gaussian fit of their distribution (dashed red) has a standard deviation of 0.1 pixels meaning that at σ ±3 , 99.7% of the measurements are within 0.3 pixels. The asphericity of the spheres used is ± 1.25 μm on the diameter, meaning ±0.15 pixels at this resolution. Assuming that the asphericity of the real spheres can be modelled as elliptical, the considerations in section 3.2.3 would imply that expected prediction errors due to asphericity are in the same same order as the observed ones and therefore that the majority of the error in the positions measured (which are incurred twice for the measurement of displacements) can largely be justified by the physical sphericity of the spheres tested 5 rather than by the quality of the kalisphera-based template matching technique. For completeness, figure 9 also reports, in dot-dashed blue, the distribution of displacement obtained using ID-Track (where the centers of the grains are obtained by the centers of mass of the separated objects), which gives comparable results for the quality of spheres used. Further tests show that when purposely 'throwing off' the initial guess values by a few pixels, the position of the spheres's centers obtained from the template matching is unaffected.
A multi-core implementation of this template matching procedure with working examples is also provided by the authors at the code repository given above.

Applications of kalisphera sphere-matching
Satisfied with the behavior of the kalisphera sphere-matching tool tested on real images in section 3.3, this technique is applied in this section to two cases of real geomechanics experiments with x-ray tomography. The first case is the study of calcite cementation between glass spheres, where a large number of spheres has been purposefully imaged at a relatively low resolution, with the objective of having a more mechanically representative specimen, in order to be able to make a micro-mechanical validation of recent work [17,18]. The low number of pixels across a diameter of the spheres in this case, means that the partial volume voxels make up an important part of the total number of voxels for each sphere (see figure 10, top row), meaning that thresholding the grains, or the cement, will induce significant errors on the boundary.
As the top row of figure 10 shows, in the presence of significant noise, and blur, the kalisphera sphere-matching converges well (middle image). When the contribution of each sphere's graylevel above that of the air is subtracted from the original input data (rightmost image), what is left behind is a cement distribution which is much easier to quantify, since spheres have been subtracted with their partial volume pixels leaving behind just cement and air (which can be quantified for example with a simple histogram-based threshold).
The bottom row of figure 10 shows another application with a partially-wet granular medium, from the work of [19]. As is visible from the images, there are more pixels for a sphere diameter in this case. The challenge in this case is essentially the same as above, with the added complication that the contrast between the air and water in these images is relatively low. A recurring problem in the quantitative analysis of multi-phase images such as these is the categorization of each pixel as belonging to grain, water, or air (i.e. the trinarization of the image). This is often problematic at the boundary between the densest (grain) and least dense phase (air), since 5 Spheres with a better sphericity are available from the aforementioned manufacturer. the partial voxels have the same grayvalue as the intermediate phase (e.g. water, cement). As shown in figure 10 the proposed method allows these voxels to be distinguished and therefore correctly categorized.

Conclusions
In this contribution we present a novel technique to produce 3D images of spheres analytically describing their partial volume effect as the solution of the possible intersections between a cube (representing a voxel) and a sphere. This technique is named kalisphera, and is expected to be of significant use for the creation of realistic, controlled images of simplified microstructures to test 3D image measurements tools on.
A sophisticated kalisphera-based sphere matching tool is also presented, based on the optimisation of a correlation potential allowing rapid convergence on ideal spheres with resolutions verified down to − 10 2 px. The analytical solution has proven necessary to achieve the desired accuracy and speed. The effect of the most common experimental imperfections of real-world imaging on the technique is quantified and the technique is enhanced for better performance on polydisperse materials. The overall reliability of the technique is tested against real x-ray tomography experimental data. Finally, the use of this tool is illustrated on two different domains of geomechanics, where the removal of the spheres constituting the granular skeleton simplifies the identification of the remaining phases (cement and water in the examples).
Work is ongoing to improve the technique: the certainty of the sphere impenetrability restricts the possible positions of the spheres. Its implementation in the minimization of the cross correlation potential would prevent erroneous overlaps and enhance the accuracy of the estimations. Extension to other geometrical shapes is under evaluation, currently with the main aim of reducing the effect of sphere imperfections on the accuracy of the method. From a coding viewpoint, a Python version is currently available in the repository provided above and an optimized C-code is under development. this paper. Alessandro Tengattini wishes to thank the University of Sydney International Scholarship scheme. Laboratoire 3SR is part of the LabEx Tec 21 (Investissements d'Avenirgrant agreement ∘ nANR-11-LABX-0030). The research leading to these results has received funding from the European Research Council under the European Union's Seventh Framework Program FP7-ERC-IDEAS Advanced Grant Agreement no 290963 (SOMEF).