<figure><img src="image_1.png" /></figure>

# Best projections selection algorithm based on constrained QDEIM for sparse-views X-ray Computed Tomography

Victor Bussy¹, Caroline Vienne¹, Julie Escoda¹, Valérie Kaftandjian²

¹Université Paris-Saclay, CEA, LIST, F-91190, Palaiseau, France,
e-mail: victor.bussy@cea.fr

²Laboratory of Vibrations and Acoustics, INSA-Lyon, Villeurbanne, France,
e-mail: valerie.kaftandjian@insa-lyon.fr

## Abstract

X-ray Computed Tomography (CT) is a non-destructive testing tool increasingly used by manufacturers, with growing interest in in-line testing applications. However, it is still struggling to establish itself as a standard technique due to long acquisition times. In order to meet industrial imperatives, recent research aims at reducing this acquisition time by limiting the number of views while keeping a good CT image quality. The rise of iterative reconstruction methods has made it possible to partially solve the Sparse-View CT problem thanks to the injection of regularisation terms. Furthermore, these methods also allow more freedom in the acquisition trajectories. In this work, we propose a method for reducing the number of projections without impacting the reconstruction quality by acquiring only the most relevant views. Indeed, not all views provide the same amount of information. We, therefore, present a technique for selecting views when the geometry of the inspected object is roughly known *a priori*. Our method is based on the *Q- Discrete Empirical Interpolation Method* (QDEIM) and considers the attenuation of the rays as an additional constraint.

## X-Ray, Tomography, NDT, Sparse-view CT, QDEIM, Best view selection, Online Control

### 1 Introduction and Previous Work

One of the main challenges in X-ray CT is to reduce the number of views to acquire without decreasing the reconstruction quality. Currently, for a 3D CT scan, projections are uniformly acquired along a circular trajectory around the object. However, this strategy is not optimal since most of the acquired information is redundant and improved quality could be obtained by choosing views that provide the most meaningful information; especially in industrial applications, where *a priori* knowledge of the object's geometry is available through a CAD model and could be used to acquire only the best projections. Although there is still no standard method to select optimal projections, many approaches have been suggested, and similarities in the selection criteria can be observed. The first main criterion relies on the idea that projections with salient edges should be preferred. Rays tangential to edges are more likely to capture the delimitation between two zones with higher contrast. This concept can be directly used as in [1][2] but can also be declined in many forms. For example, some authors use wavelet decomposition to find projections that well-capture edges [3][4]. In [5], spectral richness is used as a criterion to increase or reduce the angular step along a circular trajectory. The assumption is that zones with richer spectral importance, typically salient edges, require more projections to be well reconstructed. In [6], the next best view along a circular trajectory is defined as the projection between the two views with the maximum entropy difference. This criterion will favour acquisition around angles where the X-ray projection changes abruptly, therefore around straight edges.

The importance of acquiring projections that are aligned to edges is also found in the Tuy-Smith criterion [7], which states that for a plane feature to be visible, at least one projection has to be within this plane. In their work, Gabriel Herl and al. [8], optimise the completeness of Radon data to ensure a sound reconstruction, but they also consider the penetration length to reduce artefacts. Indeed, the second main criterion for selecting the best views is attenuation reduction. Views with a high Signal-to-Noise ratio (SNR) must be favoured to prevent metal, beam hardening or photon starvation artefacts and reduce the noise in the reconstructed image [9][10].

A convenient way to quantify the effects of attenuation is through Model Observers as a detectability index. Most of the time, they rely on the Modulation Transfer Function, Noise Power Spectrum and a frequency template of the studied features. The detectability index provides a figure of merit that describes how well a feature can be reconstructed and visible [11][12]. Although these methods all offer good results, none allows one to quickly choose the best projections for a complete reconstruction from few *a priori* information.

Based on the previously mentioned concepts, we proposed in a previous work a method for selecting views based on the Q-Discrete Empirical Interpolation Method (QDEIM) [13]. This method, based only on the geometry of the part showed good performance in terms of reconstruction quality with a calculation time below one second. In this paper, we propose a variant of the QDEIM, considering a penalty term related to the physics of the problem and thus to the signal to noise level of the projections. Section 2 describes the QDEIM and its new proposed variant. Section 3 presents simulated results on two objects, and section 4 gives an outlook and perspectives for future investigations.
## 2 Q-Discrete Empirical Interpolation Methods

Initially designed for Reduced Order Models, QDEIM makes it possible to quickly find the optimal positions of sensors to reconstruct a signal from its representation in a low-dimensional tailored space [14]. In [13], we demonstrated how QDEIM could be applied in a different way to select the best views for reconstruction. To apply the QDEIM for view selection, first a set of possible projections along a predefined trajectory is computed (see Figure 1). By abuse of language we will call this set of projections a sinogram. The projections do not need to be of high resolution nor precise, a simple ray tracer is sufficient to simulate them. The figure illustrates an initial classical circular trajectory (the blue points represent the source position along the trajectory) but the same method can be applied to any kind of initial trajectory including limited angle or 3D trajectory (spherical).

<figure><img src="image_2.png" /><figcaption>Figure 1: Illustration of the circular acquisition trajectory around an engine block in the CIVA software. Blue points represent the different source positions.</figcaption></figure>

The method performs a pivoted QR decomposition on the sinogram matrix. At each iteration of the QR factorization, we select the projection, which is the most orthogonal to the previous ones. This algebraic criterion is easy and fast to compute. We have shown in [13] that applying QDEIM to the sinogram matrix was tantamount to greedily selecting the most different projections. This property is attractive because it joins concepts evoked on the information provided by a projection [15]. In this algorithm, at each iteration, QDEIM maximises the determinant of the Gram matrix of the projections, which is a criterion for optimising the information brought by a new acquisition. We have also shown the ability of the algorithm to choose views with tangent rays to the edges. In a completely unsupervised way, the criterion of favouring projections with salient angles also seems to be optimised by the QDEIM. On the majority of the objects studied, we notice that the views most orthogonal to each other are the different faces of the object, which ultimately leads to a selection of views that favour the capture of edges and protuberances from different angles.

To improve the results of the QDEIM, we propose to inject the physical effects due to attenuation process of the X-ray beam. We therefore propose to add a constraint on the pivot of the QR decomposition. This constraint is task-specific, and in this paper we will focus on attenuation reduction. This idea has already been exploited by Clark et al. [16], to introduce a cost constraint for sensor selection. We will use the same approach to introduce a constraint on the attenuation. In order to avoid areas that attenuate the signal too much and reduce the reconstruction quality, we are going to use a constraint aimed at favouring projections where the attenuation is low. We therefore propose the algorithm 1 for the views selection.

<caption>**Algorithm 1: Constrained Pivoted QR algorithm**</caption>

**Data:** $\gamma \in \mathbb{R}$, cost vector $\eta \in \mathbb{R}^m$, $\mathbf{X} = [\mathbf{x}_1 \dots \mathbf{x}_m] \in \mathbb{R}^{n \times m}$
**Result:** $\mathbf{P}$ List of indexes corresponding to the selected columns ordered by decreasing importance
$\mathbf{P} = []$
**for** $\ell = 1$ **to** $m$ **do**
> $i = \underset{j, j \notin \mathbf{P}}{\operatorname{argmax}} (\|\mathbf{x}_j\|^2 - \gamma\eta_j)$
> $\mathbf{X} \leftarrow \mathbf{X}(\mathbf{I} - \mathbf{x}_i^T \mathbf{x}_i / \|\mathbf{x}_i\|^2)$
> $\mathbf{P} \leftarrow [\mathbf{P}, \mathbf{i}]$
**end**

In this algorithm, $\mathbf{X}$ represents the sinogram, $\eta$ is the cost vector of each projection and $\gamma$ a scalar aiming to modulate the compromise between the traditional pivot criterion and the cost introduced in this variant. The algorithm is the same as the traditional pivoted QR, only a $\gamma\eta_j$ term is added in the pivot. The cost vector is chosen to disfavour projections that would create artefacts during the reconstruction. In our implementation, the cost of a projection is proportional to the median of its pixel values, considering only pixels whose incident rays have penetrated through matter. Thus, views that have low attenuation across the entire projection are favoured.
# 3 Simulation Results

## 3.1 Results on an engine cylinder head for reconstruction on a circular trajectory

QDEIM and its constrained variant have been tested on a car engine cylinder head. Given that this part has a strong attenuation and a significant penetration length, this object makes it possible to highlight the interest in the proposed penalty term. In this example, we will sub-sample a circular trajectory of 360 uniformly distributed points and use our method without and then with a constraint. For selecting the views, the simulated projections were made with a simple ray tracer and were 116×91 pixels. For the reconstruction, the projections were 1024x1024 pixels for a volume of 498×256×756 voxels. The magnification factor is 1.2. The projections, as well as the reconstruction, are obtained using the CIVA software [17].

In this example, we will set $\gamma=5$. This value was empirically chosen to have the best compromise between the algebraic criterion and the attenuation. The norm of the projections has been normalised between 0 and 100, and the constraint vector between 0 and 1, so that the factor $\gamma$ roughly represents the percentage ratio between the two criteria. Several values have been tested for the reconstruction, and $\gamma=5$ shows promising results and an interesting choice of projections. Figure 2 shows the first four views selected by the unconstrained and constrained method.

<figure><img src="image_3.png" /><figcaption>Figure 2: Sorted first four projections selected by unconstrained QDEIM (first row, from left to right: 269°, 167°, 343°, 89°) and constrained QDEIM ($\gamma=5$, second row, from left to right: 75°, 153°, 238°, 326°)</figcaption></figure>

The first view is slightly shifted from the unconstrained version, and an area with less attenuation is visible on the left of the image. We also notice that the constrained version's first view faces four holes vertically aligned on the right of the image. This detail emphasises that the constrained version also respects the criterion of salient edges. The second view is also similar to the unconstrained version but shifted, so the size of the area with a small penetration length is larger on the left while maintaining a thorough view of the holes on the right. The third projection, however, is very different. This oblique view has areas with short penetration lengths on both sides of the cylinder head and is salient to exhaust holes. This view is particularly interesting as it brings new information with a nicer SNR than the third view of the unconstrained QDEIM. As for the fourth view, although the one selected by QDEIM is interesting from a geometric point of view, it will not allow a reasonable reconstruction of the interior of the part because of the strong attenuation of the incident X-ray beam. The constrained version chooses an oblique projection with a lower average penetration length which will surely bring more information during the reconstruction.
<figure><img src="image_4.png" /><figcaption>Figure 3: Comparison of the reconstruction for different sampling strategies</figcaption></figure>

For the reconstruction, Figure 3 shows cross-sections of the reconstructed engine cylinder head with 30 projections. From left to right, we display the results obtained with the equiangular sampling, QDEIM, and the constrained QDEIM with $\gamma=5$. The optimised samplings improve the overall quality of the reconstruction with contours of the material better defined. The voxel values of the background are closer to zero and the fine structures of the part have higher contrast, especially for the constrained version, compared to the reconstruction result based on uniform sampling.

Figure 4 shows the profile lines of the reconstructions (see red dashes in Figure 3). The boundaries between material and void areas are sharper in Constrained-QDEIM. Uniform sampling on a circular trajectory has less contrast, making it more challenging to interpret and analyse the reconstruction for inspection.

<figure><img src="image_5.png" /><figcaption>Figure 4: Line profile through the engine cylinder head</figcaption></figure>

In this example, we have chosen to select attenuation as the criterion for the penalty term. This criterion is relevant because many of the physical effects harmful to the reconstruction quality are precisely due to it. However, the proposed method is very general, and any task-specific term can be used. If we suppose that the shape or location of possible defects are known *a priori*, one can simulate a sinogram with these, calculate properties on the projections, such as the contrast of the defects, and use this as a criterion. However, in our approach, we do not wish to use the information on the defects. We will therefore use the same constraint criterion in the rest of the paper.

## 3.2 Results on a model aircraft engine for defect detection on a spherical trajectory

In this section, we will look at another example, the reconstruction of a model aircraft engine. This time, we will focus on the reconstruction of two defects inserted into the engine's stiffeners. The defects are modelled with CIVA, the first is a spherical defect, the second a crack. This object has less attenuation than the engine block previously studied, which is why the factor $\gamma$ has been reduced to 2. The engine has a 'L' shape with two very defined protuberances on two perpendicular axes. The salient views to the edges are therefore on different planes (see CAD model in Figure 5a). Thus, we have chosen to sub-sample a spherical trajectory around the part in this example (Figure 5b). The precomputed sinogram for selecting views was made of 900
low resolution projections of 128×128 pixels obtained from the theoretical model of the part (without defect), while the realistic projections used for the reconstruction were simulated from the part with the two inserted flaws, with a polychromatic X-ray beam and a flat panel detector of resolution 512×512 pixels. The corresponding volume, reconstructed with a SART algorithm, has a resolution of 512³ voxels. The uniform sampling is made from a *fibonacci lattice*, which allows to approximately sample uniformly on a sphere [18].

<figure><img src="image_6.png" /><figcaption>Figure 5: CAD Model of the aircraft engine and illustration of the spherical trajectory in the CIVA software</figcaption></figure>

Figure 6 displays the results obtained with the different sampling strategies, compared to the reference corresponding to the CAD model with the two inserted defects. The spherical defect on the left stiffener and the crack on the right stiffener are barely visible with the uniform strategy, while the QDEIM allows these defects to be reconstructed with high contrast.

<figure><img src="image_7.png" /><figcaption>Figure 6: Comparison of a reconstructed cross-section for different sampling strategies and 40 projections</figcaption></figure>

To assess more quantitatively the improvement due to the sampling strategy, we show in Figure 7 the Peak Signal to Noise Ratio (PSNR) in dB of the two inserted defects. Even if the optimised samplings are not specifically designed to improve the reconstruction of defects, we show that by improving the complete reconstruction, defects are thereby better reconstructed. As the engine has a low thickness, the penalty term is not as important as in the first example. However, the few differences will propagate in the QR decomposition and still lead to a clear improvement of the method. In Figure 7, we notice that the constrained QDEIM well reconstructs the spherical defect with only 10 projections. In comparison, other methods need over 20 projections to have the same reconstruction quality.

We report in Table 1 the PSNR values for both defects computed in the volume reconstructed from 40 projections with respectively a uniform sampling, a QDEIM sampling and a QDEIM sampling with attenuation constraint. We also compare the SSIM values computed around these defects, which allows for a comparison in terms of structural similarity to the reference image [19] as well as a contrast value, which is often considered as a reliable detectability criterion for a defect. This contrast C is computed as follows:

$$C = \left| \frac{I_{\text{zone}} - I_{\text{defect}}}{I_{\text{zone}} + I_{\text{defect}}} \right| \qquad (1)$$

Where $I_{\text{zone}}$ represents the mean intensity around the defect, and $I_{\text{defect}}$ the mean defect's intensity. An optimised trajectory can improve the reconstruction by orders of magnitude. The contrast around the defects is almost double with QDEIM compared to a uniform sampling.
<figure><img src="image_8.png" /><figcaption>Figure 7: Comparison of reconstruction quality for different number of projections for every proposed sampling strategy for a spherical trajectory</figcaption></figure>

<table><caption>Table 1: Reconstruction quality of the engine's defects for 40 projections</caption><thead><tr><th rowspan="2"></th><th colspan="2">Uniform</th><th colspan="2">Without constraint</th><th colspan="2">Constrained</th></tr><tr><th>Sphere</th><th>Crack</th><th>Sphere</th><th>Crack</th><th>Sphere</th><th>Crack</th></tr></thead><tbody><tr><th>PSNR [dB]</th><td>36.22</td><td>37.85</td><td>40.80</td><td>42.08</td><td>47.23</td><td>53.27</td></tr><tr><th>SSIM</th><td>0.901</td><td>0.921</td><td>0.954</td><td>0.966</td><td>0.976</td><td>0.983</td></tr><tr><th>Contrast [%]</th><td>38.5</td><td>18.3</td><td>62.0</td><td>34.6</td><td>62.7</td><td>39.4</td></tr></tbody></table>

## 4 Conclusion

In this paper, we propose a fast and efficient sinogram subsampling method to select the most relevant views for tomographic reconstruction. In order to take into account specific constraints and noise when selecting views, we propose adding a penalty term to the QDEIM method. We have chosen to use a constraint on the median attenuation of the projections. This technique has been used for the reconstruction of an entire automotive engine cylinder head and has shown improved reconstruction on an undersampled circular trajectory. Constrained-QDEIM was also used on a model aircraft engine with simulated defects on a spherical trajectory. The proposed technique succeeds in greatly improving the contrast around the defects. The results show the interest in this constraint term and thus motivate us to continue the study of this method, in particular, to find how to determine the parameter $\gamma$ optimally.

## References

[1] Z. Zheng, K. Mueller, Identifying sets of favorable projections for few-view low-dose cone-beam ct scanning, in: In 11th Fully3D meeting, 314–7, 2011, pp. 314–317.

[2] M. Venere, H. Liao, A. Clausse, Genetic algorithm for adaptive tomography of elliptical objects, Signal Processing Letters, IEEE 7 (2000) 176–178. doi:10.1109/97.847361.

[3] M. A., H. W., S. R., Trajectory optimization in computed tomography based on object geometry, 11th Conference on Industrial Computed Tomography, Wels, Austria (ICT 2022) (2022). URL https://www.ndt.net/search/docs.php?id=26627

[4] O. Barkan, J. Weill, S. Dekel, A. Averbuch, A mathematical model for adaptive computed tomography sensing, IEEE Transactions on Computational Imaging 3 (2017) 551–565.

[5] M. Haque, M. O. Ahmad, M. Swamy, M. K. Hasan, S. Lee, Adaptive projection selection for computed tomography, IEEE transactions on image processing : a publication of the IEEE Signal Processing Society 22 (08 2013). doi:10.1109/TIP.2013.2280185.

[6] G. Placidi, M. Alecci, A. Sotgiu, Theory of adaptive acquisition method for image reconstruction from projections and application to epr imaging, Journal of Magnetic Resonance, Series B 108 (1) (1995) 50–57. doi:https://doi.org/10.1006/jmrb.1995.1101. URL https://www.sciencedirect.com/science/article/pii/S1064186685711016

[7] H. K. Tuy, An inversion formula for cone-beam reconstruction, SIAM Journal on Applied Mathematics 43 (3) (1983)
546–552.

URL http://www.jstor.org/stable/2101324

[8] G. Herl, J. Hiller, A. Maier, Scanning trajectory optimisation using a quantitative tuybased local quality estimation for robot-based x-ray computed tomography, Nondestructive Testing and Evaluation 35 (3) (2020) 287–303.
arXiv:https://doi.org/10.1080/10589759.2020.1774579, doi:10.1080/10589759.2020.1774579.
URL https://doi.org/10.1080/10589759.2020.1774579

[9] N. Grozman, R. H. Schmitt, A. Buratti, Investigating the influence of workpiece placement on the uncertainty of measurements in industrial computed tomography, in: Proceedings of 9th Conference on Industrial Computed Tomography (iCT), 2019.

[10] A. Amirkhanov, C. Heinzl, M. Reiter, E. Gröller, Visual optimality and stability analysis of 3dct scan positions, IEEE Transactions on Visualization and Computer Graphics 16 (6) (2010) 1477–1486. doi:10.1109/TVCG.2010.214.

[11] J. W. Stayman, S. Capostagno, G. J. Gang, J. H. Siewerdsen, Task-driven source-detector trajectories in cone-beam computed tomography: I. theory and methods., Journal of medical imaging (Bellingham, Wash.) 6 (2) (2019) 025002–025002. doi:10.1117/1.jmi.6.2.025002.
URL https://lens.org/016-724-596-220-826

[12] A. Fischer, T. Lasser, M. Schrapp, J. Stephan, P. Noël, Object specific trajectory optimization for industrial x-ray computed tomography, Scientific Reports 6 (2016) 19135. doi:10.1038/srep19135.

[13] V. Bussy, C. Vienne, V. Kaftandjian, Fast algorithms based on empirical interpolation methods for selecting best projections in sparse-view x-ray computed tomography using a priori information, NDT & E International (2022) 102768doi:https://doi.org/10.1016/j.ndteint.2022.102768.
URL https://www.sciencedirect.com/science/article/pii/S0963869522001670

[14] Z. Drmač, S. Gugercin, A new selection operator for the discrete empirical interpolation method—improved a priori error bound and extensions, SIAM Journal on Scientific Computing 38 (2) (2016) A631–A648. doi:10.1137/15m1019271.
URL https://doi.org/10.1137%2F15m1019271

[15] I. G. Kazantsev, Information content of projections, Inverse Problems 7 (6) (1991) 887–898.
doi:10.1088/0266-5611/7/6/010.
URL https://doi.org/10.1088/0266-5611/7/6/010

[16] E. Clark, T. Askham, S. L. Brunton, J. N. Kutz, Greedy sensor placement with cost constraints (2018).
doi:10.48550/ARXIV.1805.03717.
URL https://arxiv.org/abs/1805.03717

[17] M. Costin, D. Tisseur, C. Vienne, R. Guillamet, H. Banjak, N. Bhatia, R. Fernandez, Civa ct, an advanced simulation platform for ndt, in: Proceedings of 6th Conference on Industrial Computed Tomography (iCT), 2016.

[18] B. Keinert, M. Innmann, M. Sänger, M. Stamminger, Spherical fibonacci mapping, ACM Transactions on Graphics 34 (2015) 1–7. doi:10.1145/2816795.2818131.

[19] Z. Wang, A. Bovik, H. Sheikh, E. Simoncelli, Image quality assessment: from error visibility to structural similarity, IEEE Transactions on Image Processing 13 (4) (2004) 600–612. doi:10.1109/TIP.2003.819861.