This configuration corresponds to the right side of Figure 1. For M measurements y _i and with an even discretization of the spectrum into N bins, the forward problem can be cast into a set of linear equations y i = ∑ j = 1 N e − l i μ i , j s j , i ∈ 1 , M , j ∈ 1 , N (2)

where l _i is the thickness of the ith measurement and μ _i,j is the linear attenuation coefficient value for the ith measurement at the discrete energy level j.

Eq. 2 is conveniently put in the matrix form as y = A × s (3)

with y = ( y i ) ∈ R M the vector of transmission measurements, s = ( s j ) ∈ R N the vector of the discretized spectrum and A = ( e − l i μ i , j ) ∈ R M × N the forward system matrix.

This configuration is illustrated at the left side of Figure 1. Each measurement contains K layers of distinct materials. Without any loss of generality, we consider that each measurement contains the same number of layers with the same material order. Only the thickness of each material layer is changed, and can be set to zero. With this assumption, the attenuation coefficient no longer depends on the measurement number i. The transmitted intensity y _i for the ith measurement thus writes y i = ∑ j = 1 N ∏ k = 1 K e − l i , k μ k , j s j , i ∈ 1 , M , j ∈ 1 , N , k ∈ 1 , K (4)

where l _i,k is the thickness of the kth layer of the ith measurement and μ _k,j (E) is the linear attenuation coefficient of the material k at energy level j. The forward system matrix in Eq. 3 is modified as A = ( e − ∑ k l i , k μ k , j ) ∈ R M × N .

When multiple piled up materials are allowed for each measurement, the genetic algorithm optimizes on two-dimensional matrices of size M × K instead of one-dimensional vectors in the single material case. Each cell of the matrices contains the current thickness of the column’s corresponding material for the line’s corresponding measurement. In comparison with [10], the main steps of the genetic optimization are modified as follows.

• Initialization: N _pop matrices are initialized with a randomly chosen value in each cell, such that the sum on each row (i.e., the total thickness of the corresponding measurement slab) remains in a chosen interval [l _min, l _max]

When only one material is allowed per measurement, two-dimensional M × K matrices are also used but with only one non-zero value per row, at the column corresponding to the employed material. The main steps of the genetic optimization are modified as follows.

• Initialization: N _pop matrices are initialized with a randomly chosen value in only 1 cell per row, and the other cells are initialized to zero

Prior to the reconstruction algorithm itself, a preliminary step of dimension reduction is realized on a typical spectrum presenting roughly the same characteristics as the unknown spectrum. As shown in Figure 2, by sampling uniformly from the cumulative integral of this typical spectrum’s derivative, an approximately optimal choice of the N energy bins is obtained, allowing an efficient spectrum representation, later employed during the optimization process of the reconstruction.

The spectrum is reconstructed on the optimal energy sampling presented above, using experimental or simulated measurements y. A candidate spectrum s _cand is initialized and then modified during the optimization process. The measurements computed through the forward model y _cand = A ×s _cand are expected to be as close to y as possible.

The optimization problem for the spectrum reconstruction thus writes: argmin s cand ‖ y − A × s cand ‖ 2 (5)

A straightforward method would consist in performing a spectrum discretization over the optimal energy sampling values E ₁, …, E _N where N is usually large to obtain a good resolution of the spectrum. However, a large N increases the degrees of freedom and makes the optimization process harder. It is therefore necessary to reduce the number of parameters used to describe the spectrum. The method chosen here is to sample the spectrum with a small number P of interpolation points. This requires the spectrum to be continuous and deprived of high variation peaks, which is a general feature of high-energy Bremsstrahlung spectra.

As shown in Figure 3, the spectrum is thus fully described by P interpolation points, which coordinates can evolve both in energy and magnitude (e.g., 2P degrees of freedom). At each optimization step, the points are interpolated by a piecewise cubic polynomial function [15], then projected over the N energy intervals. y _cand is then computed and new values of coordinates for the interpolation points can be determined.

A trust-region algorithm for constrained optimization [16] is used to compute the interpolation points corresponding to the approximated spectrum. Minimum and maximum energy levels are enforced, based on a known low energy cut-off and the maximum energy of electrons. The spectrum magnitude is constrained using a simple step function. In addition to this restriction, the spectrum is normalized by its integral during each step of the optimization process.

For this algorithm, the P abscissa of the interpolation points are fixed, which improves optimization performances and decreases execution time. However, other algorithms might be more efficient with the abscissa taken as additional degrees of freedom. In the case of this study, these abscissa values are taken as evenly spaced inside the energy interval, in logarithmic scale. The optimization algorithm is then applied to the P magnitudes of the interpolation points, with a fixed number of iterations.

The work done in [10] led to a great reduction in orders of magnitude of the forward matrix condition number, and highlighted an exponential distribution of thicknesses for the optimal measurement set using one material. Using this one-material genetic algorithm, the optimal arrangement of thicknesses plotted on Figure 4A was obtained, with a condition number of 1.8 × 10⁶.

Figure 4B displays the results obtained using the first type of transmission measurements with several materials, described in Section 2.2.1. When multiple piled up materials are allowed per measurement, a significant improvement of the previous results is observed, with a further reduced condition number of 1.61 × 10⁴.

Figure 4C shows the results obtained using a single material per measurement. With this second version of the algorithm presented in 2.2.2, an even greater improvement is observed with a highly reduced condition number of 2.11 × 10³.

Regardless of the performance of the different algorithms, another main goal of the proposed algorithms is to ensure a good stability of the optimized measurements. To compare the stability of the multi-material algorithms developed in this study, both were executed several times with the same parameters. The set of materials obtained for each optimization was then plotted.

Results are shown in Figure 5. In Figure 5A where multiple materials per measurement are allowed, the mean thickness value and the standard deviation of each piled up material for each measurement are plotted. However, in Figure 5B where only one material is allowed per measurement, the plots have to be separated because a given measurement can correspond to different materials for different optimizations.

The result of the stability comparison is clear. As illustrated on Figure 5A, the algorithm version with piled up materials is quite unstable, with a significant standard deviation on material thicknesses between experiments and a strong variance on the condition number. Conversely, for the algorithm with only one material per measurement, the stability is satisfying. As shown in Figure 5B, each measurement is attributed almost always the same material with the same thickness across experiments. The reduction of the research space size thus plays a key role in the stability improvement.

In this study, the nominal configuration for the spectrum reconstruction consists in: M = 12 measurements, N = 100 energy intervals, K = 4 different materials (iron, copper, tantalum and lead), only one material allowed per measurement in the genetic measurement set optimization (hyperparameters: 500 generations, 1000 in population size), config. Three for the simulation configuration, and p = 10 interpolation points for the reconstruction algorithm.

For this nominal configuration, the reconstruction algorithm has been applied multiple times. Reconstructed spectrums (dotted lines) are plotted on Figure 7, along with the objective theoretical spectrum (solid black line).

Finally, ablation studies were performed to evaluate the influence of every part of the spectrum estimation pipeline on the reconstruction accuracy. Figure 8 highlights the relative improvements obtained with each major step of the reconstruction method.

More precisely, reconstructions have been performed independently in the nominal configuration in the cases where:

• No anti-scatter grid was used for the measurements (Figure 8A)

As shown in Figure 7, the overall reconstruction accuracy is satisfying in the nominal configuration, but two areas can be distinguished. For high energy (E > 40 keV) there is no restriction as the constraint step function is set to unity, and the reconstruction shows great accuracy even for as few as 12 interpolation points. However, in the low energy range, the points are constrained to 0 and are thus not optimized. In practice, this is necessary because of the much lower contribution of low energies in the transmission measurements. Indeed, when dense materials are subjected to an X-ray beam of given spectrum, most of the low-energy photons are absorbed and are not detected at the end. The difficulty to reconstruct the low energy part of spectra is thus an intrinsic issue for the inverse problem at hand.

Figure 8A emphasizes the substantial impact of the presence of an anti-scatter grid in the experimental setup on the reconstruction accuracy. When nothing is done to reduce the scattered rays, the measurements are very noisy, leading inevitably to an inaccurate reconstruction.

Similarly, the influence of using multiple materials is illustrated in Figure 8B. Because only iron is allowed in the first optimization problem, the genetic algorithm cannot converge to a satisfying minimum of the condition number of the system. Thus, even with low noise, the reconstruction is unstable and inaccurate.

Finally, the interest of using the “trust-constr” algorithm for the second optimization problem formulated above is clearly highlighted in Figure 8C. When the expectation-maximization algorithm is used instead, as it has usually been done in other studies [10], the quality of the reconstruction is very poor for low and medium energy ranges, and the spectrum peak is not retrieved at the expected energy level.