Single-projection reconstruction technique for positioning monodisperse spheres in 3D with a divergent x-ray beam

The measurement of the position of single-sized spheres in 3D from a single, divergent, radiographic projection is addressed in the present study with the development of a novel method. Generally speaking, the location of the shadow cast by a single sphere on a detector defines a source-detector ray; the position of the particle along this ray is identified by the strong prior knowledge of its radius and the size of the shadow. For a dense assembly of equal-sized particles whose projections overlap, a novel Fourier transform based technique is introduced to give a first 3D determination of the particle centres. The uncertainty of this measurement is calculated from synthetic data with a known noise distribution. A further refinement of this measurement is performed based on the minimisation of the projection residual. The combined approach is validated both on synthetic data, and on real radiographs of a glass bead packing. The effect of noise on the measurement uncertainty is evaluated. The technique is made available to the community in the open source python package radioSphere.


Introduction
A collection of spheres is the simplest form that a granular material can take, yet it exhibits most of the rich behaviour that makes granular mechanics such a fascinating and active field of research [e.g. [1][2][3]. This simple description lends itself well to implementation in computer simulations (see 'discrete element methods' or DEMs stemming from [4]).
For granular experimentalists, glass spheres are also a common system to study. Given the complexity of grain kinematics, imaging methods capable of identifying individual grains are extremely pertinent [5,6] and there now exist grainbased image analysis methods which are able to characterise a granular system from such measurements [7,8]. Individual particle information within a 3D granular system is typically obtained using computed tomography. However, a significant drawback of such a technique is the requirement for the mechanically complex and time-consuming acquisition of radiographs in many directions, which limits time resolution significantly, meaning that dynamical 3D granular flows are essentially out of reach of tomography. One existing method for probing the kinematics of a granular flow is to interrupt the flow very quickly (see for an example of silo flow [9]), record a tomography, and possibly trace backwards in time the position of the beads. Many alternative imaging techniques are beginning to fill this gap: x-ray radiography-based rheography [10] (which does not resolve individual particles directly) and x-ray-or positron-based intruder detection [11,12] (which require few, dense and/or radioactive particles as markers in the flow) are notable examples. Alternatively, grains can be immersed in a viscous fluid with a matching refractive index [13] or characterised by MRI [14] to track them individually or measure the velocity field. Several recent works have used an initial tomography scan of a granular packing and updated particle positions only with a few radiographs of each subsequent imaged state [15], itself a discrete version of [16,17]. There also exist sophisticated techniques to recognise specific sets of shapes in three dimensions from single divergent radiographs [18]. Finally, it is important to mention a similar method for parallel projections of granular media, where boundary conditions are imposed to regularise the displacements in the x-ray direction [19].
Here, a novel technique is described that exploits knowledge of particles shape and size (i.e. a strict requirement for spherical particles of a known single size) to reconstruct the 3D positions of each particle from a single radiograph acquired with a divergent beam. The uncertainties in the measured positions are evaluated both with synthetic data with controlled noise, and with real experimental data (with respect to a tomography image). Since one key application of this method is the analysis of granular flows, solving the problem directly from a single radiograph, rather than an initial tomography, is a challenge to be faced. Indeed, it may not be possible to acquire a tomography scan in a number of flow geometries (and even if it is, it is likely that studying a well-developed flow is more interesting than the immediate vicinity of the static part).
This paper is organised as follows: In section 2 the geometry of the problem is defined. In section 3, an FFT-based sphere detection algorithm called tomopack is developed for a parallel beam, and in section 4, tomopack is used to make a 3D guess of particle positions by scanning through different zoom levels for a divergent cone beam. The tomopack algorithm is then applied to synthetic data generated from DEM sphere packings in section 4.1. The initial 3D guess of the particle locations is then improved with an optimisation algorithm developed in section 5. Finally, the combined tomopack and optimisation techniques are validated on an experimental case in section 6.

Imaging geometry
The imaging system that is modelled herein is a divergent laboratory x-ray (transmission or reflection) source emitting a 'cone beam' which is detected by a relatively large 2D photosensitive detector. Although specifications vary, such sources can achieve emission cone half-angles as high as 70 • , meaning that in principle a 'source-detector distance' (SDD) of the order of the detector size is possible, although not often used in x-ray tomography. Micro-focus x-ray sources can have extremely good focusing of the electron beam onto the target, meaning that the size of the resulting x-ray emission spot can in some cases be below a micrometre in size. A significant advantage of the combination of a focused, divergent x-ray source and a flat-panel detector is that thanks to geometric magnification such imaging systems are able to trade off spatial resolution and field of view very easily by translating a sample closer to the source or detector respectively.
As opposed to the parallel geometry typically offered at a synchrotron, the radiographic projections of a 3D object acquired on a divergent system has a differing level of  Projection of spherical particles with a conical x-ray beam. Two spherical particles of radius r are projected individually onto a detector panel. The blue particle, located closer to the x-ray source, appears larger on the detector panel due to the conical beam. magnification in the direction of the beam-this is necessarily taken into account in tomographic reconstruction techniques such as the famous FDK algorithm [20]. The technique proposed in this work aims to position spheres in three dimensions: simply put, with knowledge of the imaging geometry and the physical size of a sphere in a radiograph, the peak of the absorption spot on the detector defines a source-detector line (as shown in figure 1) on which the centre of the sphere must lie; positioning the sphere along this line can then be achieved by detecting the size of the absorption spot on the detector, which is related to the level of magnification of the sphere, and thus its position along the line.
In order to describe the problem mathematically, the coordinate system is defined in figure 1. The centre of the ith sphere is denoted with capital letters as X i in three dimensional space. The first coordinate X is along the centre line of the x-ray beam path, and the remaining directions Y and Z are perpendicular to it and each other.
Projected coordinates on the detector panel are denoted in lowercase, x = {y, z} in mm. The sphere radius is denoted r and is assumed to be identical for all spheres, although it would be desirable to generalise the approach in the future to a set of discrete particle sizes. The projection of a sphere in a divergent beam is schematically represented in figure 2. With different magnification factors (distances between the source and the sphere), the width of the absorption peak recorded is altered, but the magnitude of the peak is not affected.
It is important to note that the projections treated here are geometrical in nature, meaning that the scalar 'measured' on the detector is the distance travelled through the material (i.e. L in mm). This is clearly not the raw output from a radiograph in the lab, requiring at the very least the natural log of the intensity-normalised radiograph ln(I/I 0 ) and a suitable calibration. In the case where polychromaticity of the beam is significant, an absorption calibration with a large object composed of the same material as the particles that will be studied is suggested, allowing an appropriate calibration function of a known distance against the measured I/I 0 . This calibration procedure is reported in the appendix for the experimental validation in section 6. Although keeping projections as calibrated values of L is convenient for a quantitative comparison to projected or synthetic data, the signal-to-noise ratio is defined in greylevels in terms of I/I 0 , as the difference in values between the background and the maximum value through a single particle (its diameter) divided by the standard deviation of the background.

Finding sphere centres in a radiograph
Initially a parallel beam is discussed, instead of a divergent beam. The connection between this derivation and the use of the technique for a divergent beam is made in section 4. For a parallel x-ray beam, the projection of one disk centred at Y = Z = 0 (and therefore y = z = 0), in units of mm, will be recorded on the detector as ψ(x), where Hence the projection of the entire pack of N particles is easily written as Such a function is shown in figure 3. Let us note that it can be rewritten as where ⋆ denotes a convolution, and I denotes the indicator of the projected particle centres as depicted in figure 3. A first question to address is to estimate I(x) from the known p(x). In order to solve this 'deconvolution' problem, it is straightforward to go to Fourier space, and Example of the convolution process to construct an image following equation (3). A set of 30 randomly positioned 1 mm diameter particles is projected in a parallel projection onto a detector panel. Left: ψ, the projection of a fictional particle located at Y = Z = 0 (colour is L in mm). Centre: I, the indicator which represents the spatial position of each particle, with non-dimensional units such that a value of unity represents a single particle with centroid at that pixel. Right: p = ψ ⋆ I, the constructed image (colour is L in mm).
for all k. Hence, it appears trivial to write the solution as Although this expression is mathematically true, the inverse ofψ(k) is ill-behaved and any algorithm simply based on this expression appears to be highly unstable. The ill-behaviour of ψ(k) −1 can be traced back to its (quasi-)divergences at some wavenumbers as shown in figure 4 (i.e. values becoming vanishingly small), hence the Fourier transform I(k) contains 'gaps', in the sense that some wavenumber amplitudes should be treated as unknown, in order to avoid noise amplification in the measured data.

The tomopack algorithm
To be able to compute the inverse Fourier transform f(x), all f(k) should be known including its gaps. To find these missing values, prior information on f(x) can be used. Mathematically, it is positive, mostly 0, but contains a collection of δ functions whose amplitude is quantified (an integer number). If the unknown Fourier transform amplitudes are set to an arbitrary value, then those properties have no chance to be obeyed. It is straightforward to project any given f(x) onto one such indicator having the desired properties. For this purpose, the projector Π is defined, which operates on any f(x) to producef(x) such thatf (x) = round(pos(f(x))) Algorithm 1: Algorithm for a staggered iterative scheme (FFT means discrete Fourier transformation, and iFFT its inverse transformation).
Input: Projection data, p(x), shape function ψ(x), trusted values Ktrust, convergence factor tol, sub-relaxation parameter λ Output: where 'round' indicates rounding to the nearest integer, and 'pos' represents just positive values, such that values of f(x) < 0 are projected to 0.
Thus the following algorithm is proposed. First, wavenumbers are classified into two categories k ∈ K trust when |ψ(k)| > ψ min , and k ∈ K discard if not. The former may be trusted but not the latter ones. Then, the algorithm summarised in algorithm 1 proceeds sequentially by enforcing thef(k) values only for trustworthy wavenumbers where the information is known, and then projecting f(x) onto an indicator which is physically admissible. These two steps are repeated until a fixed point solution is obtained.
It is to be noted that the projection Π[ f ] may induce a large change on f, eventually leading to trapping the solution at an unsatisfactory position, but where the correction of the next iteration is cancelled during the projection. For this reason, a softer condition is chosen in the form of a sub-relaxation parameter 0 < λ < 1 such that the last line in the loop only accepts a fraction, λ, of the correction proposed by the projection. The full correction is obtained for λ = 1, but a smaller value (0.5) is preferable. In the example shown in figures 3 and 4, the position of 30 randomly located particles is estimated using this algorithm. After tomopack is applied, the recovered indicator is used to compute an estimated radiograph, and the signed residual shown in the bottom right of figure 4. All particles were correctly identified within 1 pixel of their initial location.

Conical effects
The use of a cone beam necessarily implies that the further from the centre-line of the x-ray beam a particle is imaged, the more it will appear deformed in its projection on the detector panel. The following parameters define the projection geometry for a given particle of radius r at location X p = [X p , Y p , Z p ] as is relevant to this deformity: ϕ, the half angle described by the projection of a particle on the detector panel, Figure 5. Performance of tomopack in finding a single particle as a function of aspect ratio. Grey area indicates ± one standard deviation.
θ, the half angle describing the width of the cone beam, and AR, the projected aspect ratio of the particle.
The behaviour of tomopack as a function of the aspect ratio is shown in figure 5, with decreasing positional accuracy with increasing aspect ratio. This issue could be ameliorated by either choosing an imaging geometry with small aspect ratios, or by remapping the flat detector panel onto an imaginary spherical detector panel centred on the source so that all particles appear as spheres. Failing this, another geometry suited to the sample, e.g. cylindrical for a collection of particles in a column could partially alleviate the problem. In the work presented below, the first approach has been used, which will additionally apply to most imaging possible with a standard micro-focus x-ray source with θ ≲ 25 • and ϕ ≲ 5 • .

Initialising a 3D guess from tomopack
If the tomopack algorithm were to be used to detect singlesized spheres in a divergent beam, the key question is its tolerance to deviations in magnification. Simple numerical experiments are performed with a single particle, where the structuring element ψ is kept constant and the particle is gradually moved in the X direction; tomopack is run and the sphere is considered as detected if the resulting indicator is larger than 0.25 in the known position of the particle centre. The shaded area in figure 6 shows, for different beam angles, the change in size of the projected particle on the detector for which it is still detected by ψ. This reveals that Figure 6. Limits on particle detection with tomopack by varying zoom level with a constant ψ with a single centred 1 mm particle. Limits of sphere detection (in units of radius) with varying half-beam angle for particle radius .
the tomopack algorithm is sensitive to changes in size of the projected particle (regardless of beam angle), with a ±0.3 % deviation in radius tolerated.
This relatively narrow, but non-zero tolerance for discrepancies in magnification between ψ and p in effect means that in a divergent projection of a granular assembly, different ψ can be used to scan the range of expected magnifications (i.e. particle positions in the X direction) in the experiment. Each particle therefore is expected to appear for a number of different ψ magnifications (corresponding to ±0.3% variation in particle size), with the best match being in the middle of the range. A rough 3D guess of particle positions can thus be obtained by scanning the divergent projection of a mono-sized granular assembly with gradually varying ψ. Since the algorithm is sensitive to changes in projected size, the ±0.3% change of projected diameter corresponds to different displacements along the beam axis as a function of the beam angle.
The sensitivity of tomopack can be exploited to differentiate particles at different levels of magnification by varying the magnification of ψ. In this case a series of indicators I is computed for different magnifications of ψ which can then be analysed for the presence of particles. Figure 7 shows a ψ scan either side of the correct value for a synthetic case of a single centred sphere, and presents the max projection of the resulting I image series in X, Y and Z directions. The Y and Z projections in particular show how the value of the indicator function increases and localises into a point around the correct value, facilitating its identification. At the first order the highest values in the I-series can be selected as detected particles (possibly imposing a non-overlapping constraint), however in the current implementation, the particle locations are identified by matching the characteristic converging cones in I by convolution.
Therefore a divergent radiograph of a sphere packing be can analysed by varying the X position of the centred synthetic sphere projection used to generate ψ and identifying the best X position for each particle. For each particle, the best X position, taken together with the detector coordinates y, z-which can easily be converted to Y and Z-yields an initial guess of the 3D position of the particle.
The accuracy of this guess will depend on the size of the X steps in the ψ scanning and the limited accuracy of particle position on the detector.

Validation with synthetic tests
Validation of the ability of the tomopack technique to find 3D particle locations was performed against artificial radiographs produced from discrete element method simulations using MercuryDPM [21]. Packings of 1 mm diameter particles were produced at a solid fraction of 0.6 (near the random close packing limit for monodisperse spheres) to simulate a dense packing of grains, with negligible overlaps between particles (much less than one pixel), which is sheared over time to generate many realisations of grain locations. To generate an artificial radiograph, the set of grains that have centres within a test domain are selected and projected using the same mechanism as described above.
The artificial radiographs were produced assuming that ϕ = 1 • , the pixel size on the detector was 0.1 mm, the detector resolution was 512 × 512 pixels and a zoom factor (SDD/SOD) of 5, which implies a half cone beam angle of θ ≈ 20 • at the edges of the sample. Particles were only sampled if they were within a distance of ±2 mm from the X axis, and at varying distances in the X direction. In this way, different typical optical depths could be investigated. The optical depth in figure 8 is defined as the distance in mm in which particle centroids should lie in the X direction to be sampled. Additionally, artificial noise was added to the artificial radiographs which was normally distributed with a mean of zero and a known standard deviation. This standard deviation is reported in figure 8 as the mean noise level.
A particle is defined to be 'lost' if the measured centroid of the particle is not within a distance of half a radius of the true location, thus treating X and (Y, Z) on the same footing. With respect to finding particles, it can be seen that the tomopack algorithm functions at > 99% efficiency for optical depths up to 5 mm with low noise levels. As the noise level and/or optical depth increases, the performance of the tomopack algorithm decreases.

Real-space optimisation with 'sensitivity fields'
In order to improve a 3D guess of particle positions (starting from above, or from a previous tomography scan where particles have been labelled), an optimisation to minimise squared residuals is carried out.
Algorithm 1 outlines the implemented procedure, which iteratively attempts to explain the current residuals as a Algorithm 2: Algorithm for local and individual optimisation of guessed 3D positions to minimise a residual computed on the projection.

Input:
Measured projection data Pm, N particle position initial guesses X 0 , particle radius r in mm, perturbation to apply ϵ in mm, convergence factor tol Output: Updated particle positions X initialisation; Compute complete synthetic projection of current particle positions: q(X n old , r, WholeDetector); Compute current squared residual: for n = 1 to n = N do Compute detector ROI for particle n Compute reference local projection: Perturbation in direction d: Local residual for perturbation: On ROI, solve for weights vector w: Update current guess: Correct X by adding displacements as function of overlap end |δX| ← |X old − X|; X old ← X end combination of 3 synthetic residual fields, which are the perturbations of the current guess of each particle in each direction. Figure 9 illustrates one step in this approach, showing the signed residual on the whole radiograph between the input radiograph and projection of the current guess in the top left and in the remaining plots the local sensitivity fields in X, Y, Z directions with a perturbation ϵ of 15 mm in X and 1 mm in Y and Z. As per the algorithm, a combinations of these three fields will be sought to best match the residual for this iteration.
Since many calls to a particle projection function are made in this procedure, it is solved locally on a Region Of  Interest (ROI) of the detector pixels, by identifying the pixels concerning a given particle with an added margin, and only projecting the detector pixels concerned, and thus only solving the problem on those pixels. To calculate an appropriate perturbation in the X direction, ϵ X , the following relationship can be used to relate this perturbation to the equivalent perturbations in the detector plane, ϵ X = ϵ Y SOD r . In this work, where not indicated otherwise, an initial perturbation of ϵ Y = ϵ Z = 1 pixel on the detector panel is assumed. Figure 10 shows the evolution of the step size |δX| and the squared residual per pixel in a synthetic case of a noisy radiograph containing a single centred sphere of radius 1 mm, with an incorrect initial guess that is offset by 0.5 mm in X, Y, and Z. This shows that the proposed algorithm is able to converge for relatively poor initial guesses in this ideal condition.
Studying a more complicated case-the synthetic reference case with 173 spheres-an initial guess is needed. In this case the approximate 3D position resulting from the X-direction scan with the tomopack algorithm is used as an initial guess. Figure 11 shows the updated residuals after 100 iterations of the algorithm with two different colour bar ranges. The mean positioning error is 0.076 mm, with SD 0.069 mm obtained on the input radiograph with an SNR of 60.

Experimental validation
In order to prove the robustness and applicability of this algorithm, it is also validated on a real experiment.
Radiographic acquisition is performed at the detector highest speed setting (i.e. at 60 Hz which imposes 4 × 4 binning and thus an effective pixel size of 0.508 mm on the detector), in order to validate this technique for imaging of dynamic processes. As usual in x-ray imaging, the 'dark field' of the detector is measured and subtracted from all subsequent measurements.
The RX-Solutions Easytom x-ray scanner in SIMAP (Grenoble) was used for this validation. Interestingly, two Hamamatsu x-ray sources are available on this machine, as listed in table 2. Both are suitable for this validation, and although a larger maximum half-angle θ is available on the transmission source, the reflection source is selected for its much higher flux, since fine focus is not needed for this validation, given the relatively large particle size.
Since both tomopack and the sensitivity optimisation have been discussed (and programmed) with P in mm, the acquired experimental data needs to be converted into this description. To this end, a larger calibration sphere with diameter 7 mm, is also scanned allowing the attenuation mm −1 curve to be fitted beyond 1 particle diameter -there is a strong assumption that the material of the calibration sphere and the material of the spheres studied are the same. To simplify this calibration, the 130 kV beam is strongly hardened with a 0.5 mm Cu filter. The photon flux available on this x-ray source means that the detector can be run in low-sensitivity mode, which helps to reduce shot noise. The source-detector distance was 242.6 mm and the source-object distance around 23 mm. The maximum cone beam half-angle is approximately 21.5 • and the limits of the beam are evident in the radiographs, as dark zones. ϕ for the middle of the sample is 2.5 • .
The following datasets are acquired, each time allowing a 5 min stabilisation of the x-ray source.
• A single 'flat field' I 0 (the image of the beam with no object) averaging 64 images. • A single image of the 7 mm soda-lime glass calibration sphere. • 360 radiographs as the empty sample holder is rotated (continuously) around 360 • , again averaging 64 images.  A few of the acquired radiographs, as well as the experimental setup are illustrated in figure 12. The radiograph of the 7 mm calibration sphere is normalised by I 0 , whereas the radiograph of the column of 2 mm spheres in the holder is directly normalised by the empty sample holder at the same rotation angle. The natural log of the normalised images is computed, and for the calibration sphere the path length L inside the sphere vs. these image values is also computed (see appendix). This fitted function (L in mm vs log(I/I 0 )) is then applied to the natural log of the images acquired of the column of spheres, finally resulting in a projection P in mm.
Since radiographs around 360 • have been acquired, the data is also tomographically reconstructed, which offers a convenient validation of the quality of the 3D positions obtained on the first radiograph with tomopack, as well as after the optimisation. The tomography data is analysed with spam [8]: reconstructed grey values are thresholded and particles separated 4 from ballandrollerstore.com. using a markers-based watershed, thereafter centres of mass of the particles are computed in pixels and converted to mm with the known projected pixel size of 0.05 mm.
Thereafter each radiograph (at different angles and averaging amounts) is normalised by the radiograph of the sample holder acquired by averaging 64 measurements, and the fit applied to the log of the image to make this an experimentallymeasured p of L in mm. The normalisation of low-average radiographs with the high-average holder reveals some movement artefacts (see figure 14), which will induce bias in the higher-noise experimental images. It is expected that the measured positioning errors for the higher-noise image are an overestimation compared to a more suitable normalisation. The background of p (zones within the source cone but without particles) is fitted with a bi-linear function, which is then subtracted from p, with the explicit goal of improving X direction iterations for the sensitivity field. A tomopack guess scan is obtained by varying ψ, and passed to a sensitivity field optimisation, which is run until changes of position fall below 5 µm.
3D particle centres are compared by subtracting the mean position from all three datasets (labelled, tomopack, optimised) and relabelling the centres in the labelled image to their closest corresponding particle from the tomopack scan. The rigid-body motion of the labelled centres that minimises Results are shown for all averaging levels (and thus SNR levels) in figure 13, revealing that the tomopack 3D guess already offers a good estimation of 3D positions with respect to the centres obtained from tomography. As expected, errors increase with noise level, and the direction normal to the detector, X, is the most error-prone. For the lowest noise level, absolute errors averaged over all particles for the tomopack guess are X, Y, Z = 0.061, 0.016, 0.025 mm and after optimisation 0.068, 0.008, 0.011 mm. As a reminder, these experimental errors are to be compared to the particle radius (1 mm), the detector pixel size (0.508 mm) and the projected pixel size in the middle of the sample (0.05 mm), which is the voxel size of the reconstructed tomography volume. Using the projected pixel size as a reference, the error in the X direction from tomopack is slightly above this dimension, and the optimisation step worsens the guess slightly, most likely due to the slightly inhomogeneous background. The Y and Z errors from tomopack are respectively a third and half a project pixel-the difference likely due to the sample being longer in the Z dimension and thus the particle aspect ratio increasing away from the centre. After the optimisation step (not sensitive to aspect ratio), both error are about a fifth of the projected pixel size, which is a very satisfactory result.
Furthermore, it is worth noting that the purchased particles are ball bearings with a precision grade of 100, meaning a 'nominal ball diameter tolerance' of ±0.0127 mm, which also puts the positioning errors into context. It is expected that there will be a positioning error in the X direction due to incorrectly assumed radii, inversely proportional to the beam angle. With a perfectly calibrated p, particles with incorrect radii would appear as circular patches in the residual.  Interestingly, as the SNR decreases from 62 to 32, the measured errors are essentially constant, indicating that some artefacts not described by the model are limiting accuracy (scattering artefacts are visible in the top row of figure 14). Below an SNR of 30, noise begins to limit accuracy-the tomopack guess is more sensitive to the noise (in the case of SNR = 24 the optimised position error is still close to the lowest noise one). At an SNR of 12 -which means imaging at 60 Hz-the errors for tomopack are X, Y, Z = 0.202, 0.025, 0.027 mm and the optimised positions are 0.127, 0.024, 0.037 mm.
For reference, the residual images obtained with the tomopack guess, and after the sensitivity field optimisation are presented in figure 14 for the highest and lowest noise cases (1 and 60 Hz imaging respectively). In the case of the lowest noise (top), residuals are very low, however on the right, it can be seen that they are far from zero: for large values of L (i.e. with significant overlaps) there is an underestimation of L in the experimental image, and around the sample there is a scattering corona. The underestimation of L is likely due to beam hardening. In the higher noise case, it is clear that there is an artefact induced by the different averaging of the sample holder, which causes problems visible on the ± 1 radius scale of the left and middle images, however it seems that the optimisation is able to converge successfully despite this significant bias.
Although the analysis of a time-series could be performed by making displacements guesses and only using the optimisation, in order to evaluate the robustness of the combination of tomopack and the optimisation step, the entire set of 360 radiographs is analysed individually. This means that a large number of different particle overlap configurations are tested. Until SNR 24, excellent tracking is obtained, after which the loss of a few particles in a few views damaged the tracking obtained. As an illustration of the results that can be obtained, the measured trajectories for the lowest-noise case are shown in figure 15. The figure shows good tracking for all particles, and mostly circular trajectories, and confirms that the uncertainty in the X direction is higher than the others.

Conclusions
The combination of the tomopack algorithm and the optimisation method mean that mono-disperse spheres of a known size can be placed in 3D space from a single divergent radiograph, as can be acquired on any laboratory x-ray scanner. This adds tremendous time resolution to measurements of particle position, as well as removing the requirement to rotate the sample, at the cost of some positioning uncertainty, especially in the direction of the beam. Between the synthetic validation cases and the simple experiment presented it turns out that the measurement of a radiograph containing the distance travelled through the scanned object is quite delicate: the application of a measured attenuation profile on the same material is apparently not sufficient, and small variations of the source can cause problems (especially for the optimisation step), better knowledge and thus modelling of flat fields is doubtlessly needed, as per [22].

Perspectives
The perspectives coming from this work are split into two categories: possible uses of this tool as-is, and further developments that could improve performance and release some assumptions.
In the case of a parallel beam, the tomopack technique could be a very effective (again with ψ scanning) tool to count particles and their sizes. This may also be achieved with divergent acquisitions at different zoom levels. The combined tomopack and optimisation, could and will be used to make previously-impossible measurements of 3D particle kinematics in a number of fields such as the flow of grains down an inclined plane (kinematic interaction with an obstacle), hydrodynamic suspensions of particles (kinematics of viscous resuspension), or sheared granular media (measuring granular temperature and vortex structures).
The perspectives for development of this technique are manifold. First of all, it must be mentioned that the requirement of single-sized spheres is a strong limitation that can be lifted in a number of different ways: in the case of a two-(or more) source and detector imaging system, the limitation of monodispersity can likely be removed with ease. Otherwise, if the initialisation of 3D positions from a tomopack scan is abandoned and replaced with particles labelled in an initial tomography, the sensitivity field optimisation process should be able to successfully converge for displacement unknowns for spheres of different sizes, and may be able to yield good results also optimising three rotation unknowns for sufficiently unique shapes, in the style of [15]. The initialisation from a tomography obviously requires all the particles to be in the field of view, and cannot tolerate particles being lost, so will require small transformations between images.
More prosaically, some short-term improvements to tomopack that could be implemented with relative ease are the measurement of an experimental ψ to be used, rather than a synthetic one (possibly at different positions), which would include blurring as well as scattering effects, which might render the particle identification more robust. Furthermore, for very large θ angles, there is distortion of particles towards the edges of the detector, which can be faced with some unwarping of the image (or warping of ψ). The effect of radiograph de-noising by filtering (perhaps directly in Fourier space) might increase the signal-to-noise ratio of input data in a way that helps both tomopack and the subsequent optimisation.
The sensitivity optimisation has been found to be lead astray by missing particles, so some way of enforcing edges in the optimisation may help. Furthermore, small offsets in the calibrated p in units of L cause offsets in the measured X position, meaning that x-ray imaging artefacts such as scattering, source movement, and beam hardening all have deleterious effects on the optimised result. Better characterisation of the source will certainly be a way to make improvements in this direction.

Data availability statement
The data that support the findings of this study are openly available at the following URL/DOI: https://gricadgitlab.univ-grenoble-alpes.fr/ttk/radioSphere.

Appendix A. Fitting attenuation to path length L
For convenience, a solution is provided for the path length L of the sphere radius r located at Y = Z = 0 for a ray at angle θ. The angle β, indicated in figure 16, can be calculated as sin β = SOD sin θ r .
With this definition, the path length L can be calculated as L = 2r cos(β).
This computation has been performed for a calibration glass sphere, and when the beam is filtered the Beer-Lambert law well fits the relationship between this path length measure and the measured attenuation on the detector panel, as shown in figure 17.