# Monte-Carlo-Based Estimation of the X-ray Energy Spectrum for CT Artifact Reduction
## 2 Materials and Methods
### Maximum Likelihood Expectation Maximization Method (MLEM)  
#### 2.2.1 Theory and Mathematical Calculations  
To validate the Monte Carlo simulation results, the MLEM method was used. MLEM is a mathematical optimization method that can be applied to the reconstruction of the X-ray energy spectrum inside the flat panel detector from transmission measurements of the X-ray radiation through an object with known shape and material composition [32; 33; 34; 35]. The object usually used for this purpose is a so-called \"step wedge\" phantom, which is composed of a uniform material with different steps, each with a different thickness. Such a sample allows us to write down a discrete system of linear integral equations on the X-ray absorption through the object. When solved, it provides the unknown energy bins that compose the original X-ray spectrum of the source.  
The polychromatic X-ray image formation can be described as follows:  
\[I(x,\ L_{r},\ \epsilon)=I_{0}\int_{\epsilon_{min}}^{\epsilon_{max}}D(\epsilon)S( \epsilon\ )exp\ \bigg{(}-\ \int_{L_{r}}\mu(x,\ \epsilon)dx\bigg{)}\ d\epsilon \tag{1}\]  
with \(I(x,\ L_{r},\ \epsilon)\) the measured X-ray photon intensity at a detector pixel \(r\), \(L_{r}\) the path length from the source to that detector pixel, \(I_{0}\) the total measured X-ray intensity without a sample, \(\epsilon\) the energy of the X-rays, \(D(\epsilon)\) the detector response function, \(S(\epsilon)\) the normalized source spectrum, \(\mu(x,\ \epsilon)\) the attenuation coefficient at location \(x\) and energy \(\epsilon\). By discretizing the spectrum and defining \(W(\epsilon)=I_{0}D(\epsilon)S(\epsilon)\), one obtains:  
\[p(x,\ L_{r})=\sum_{e=1}^{E}w(\epsilon)exp\ \bigg{(}-\ \int_{L_{r}}\mu(x,\ e)dx \bigg{)} \tag{2}\]  
where \(e\) is number of energy bin. Since the composition and dimensions of the step wedge are known, we can create a simulated polychromatic projection of the step wedge, \(\widetilde{p}\), using this model for any estimated spectrum \(\widetilde{w}\). If the step wedge has \(N\) total steps and we denote the thickness of the \(n^{th}\) step as \(l_{n}\), then we have:  
\[p_{n}=\ \sum_{e=1}^{E}w(e)exp\ (-\ l_{n}\ \mu(e)) \tag{3}\]  
with \(p_{n}\) the measurement through the \(n^{th}\) step, assuming parallel beam projections, and \(\mu(e)\) the attenuation coefficient of the material of the step wedge at energy \(e\). When defining \(A_{en}\) as:  
\[A_{en}=\ exp\ (-\ l_{n}\ \mu(e)) \tag{4}\]  
the following linear system of equations needs to be solved for the unknown spectrum \(w\):  
\[p_{n}=\sum_{e=1}^{E}w(e)A_{en} \tag{5}\]  
The matrix \(A\) can be precalculated based on the available data on material and phantom dimensions. The MLEM method is then applied to solve the linear system. This method has a multiplicative update step which guarantees positivity of all coefficients [32]. The \(k^{th}\) iteration step of the MLEM method is then provided by:  
\[w^{k}(e)=w^{k-1}(e)\frac{\sum_{n}A_{en}\frac{p_{n}}{\sum_{\
u}A_{\
u n}w^{k-1} (e^{\prime})}}{\sum_{n}A_{en}} \tag{6}\]  
Note that this method cannot fit to characteristic lines, since the update step cannot introduce discontinuities, and, as such, these have to already be present in the initial guess (\(w^{0}\)).