**1. Introduction**

Dual-energy computed tomography (DECT) has found applications in clinical and industrial settings. In current DECT, one generally acquires data of low- and high-kVp X-ray spectra over a full-angular range (FAR) of 2*π*, or over at least a short-scan angular range [1–4]. Interest remains in the development of DECT imaging over limited-angular ranges (LARs) that are considerably less than the FAR of 2*π* (or than the short-scan angular range,) because such LAR scans may bear implications for radiation dose reduction, scan time minimization, and collision avoidance between the scanner and the imaged object. Inspired by the directional-total-variation (DTV) work on image reconstruction from LAR data in conventional single-energy CT (SECT) [5], we have investigated image reconstruction previously from LAR data in DECT [6,7] by focusing on the correction only for LAR

**Citation:** Chen, B.; Zhang, Z.; Xia, D.; Sidky, E.Y.; Gilat-Schmidt, T.; Pan, X. Accurate Image Reconstruction in Dual-Energy CT with Limited-Angular-Range Data Using a Two-Step Method. *Bioengineering* **2022**, *9*, 775. https://doi.org/ 10.3390/bioengineering9120775

Academic Editors: Paolo Zaffino and Maria Francesca Spadea

Received: 25 October 2022 Accepted: 29 November 2022 Published: 6 December 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

artifacts and using DTV constraints in the reconstruction of kVp images followed by an image-domain decomposition. Other methods have also been developed for DECT with LAR data, but the angular ranges are generally not smaller than 90◦ [8,9].

In this work, we propose a two-step method to reconstruct quantitatively accurate images in DECT from LAR data by correcting for both BH and LAR artifacts, thus improving the quantitative accuracy of images reconstructed and physical quantities estimated. In the method, from LAR data of low- and high-kVp, a data-domain decomposition (DDD) algorithm [10] is used first for obtaining LAR basis data in which the BH artifacts are compensated for; and a DTV algorithm [5] is then developed and tailored to reconstruct basis images from LAR basis data obtained. The reconstructed basis images can be combined to form virtual monochromatic images (VMIs), i.e., the X-ray linear attenuation coefficients, for visual inspection, and can be used also for estimating physical quantities such as iodine-contrast concentrations and effective atomic numbers within the imaged object [11–14]. We hypothesize that images and physical quantities with both BH and LAR artifacts corrected for in LAR DECT are quantitatively comparable with those obtained in FAR DECT. Therefore, in this work the results obtained for LAR DECT are compared with those obtained from FAR data in DECT in which BH artifacts are corrected for by using the DDD algorithm.

Numerical studies are conducted with a chest phantom [15] and a suitcase phantom [6] containing distinct anatomies and structures of potential relevance in medical and security applications [11,16–19]. Low- and high-kVp data are collected with single-arc (SA) or two-orthogonal-arc (TOA) scans of LAR [6], ranging from 14◦ to 180◦. Using the DDD and DTV algorithms, we estimate basis data and then reconstruct basis images, followed by the formation of VMIs at energies of interest from the basis images reconstructed. In addition to visual inspection and quantitative analysis of VMIs obtained, we also estimate iodinecontrast concentrations in chest images and effective atomic numbers in suitcase images from data of different LARs. Furthermore, we investigate image reconstructions from data acquired with SA and TOA scans of possible implications for potential non-diagnostic imaging applications involving, e.g., C-arm DECT, in which workflow or safety concerns may limit the scan angular range. The two-step method and the study design in the work can also be applied to investigations concerning image reconstruction in DECT and multi-spectral CT using techniques with sandwiched detectors [2], sequential scans [20], or advanced photon-counting detectors [21,22]. DECT with fast-kVp-switching X-ray tubes can also collect approximately overlapping rays [4].

### **2. Materials and Methods**

## *2.1. Scans of Limited-Angular Ranges*

In this work, we consider single-arc (SA) or two-orthogonal-arc (TOA) scans in a fan-beam DECT, as shown schematically in Figure 1a,b. The SA scan includes a pair of completely overlapping arcs of LAR *ατ*, whereas the TOA scan includes two pairs of completely overlapping arcs of LARs *α*<sup>1</sup> and *α*2. For each pair of the completely overlapping arcs in the SA and TOA scans, low- and high-kVp data are collected over one of the paired arcs. In this work, we assume that the *x*- or *y*-axis intersects with the middle point of each pair of the completely overlapping arcs, and that the tangential directions at the middle points of the two pairs of completely overlapping arcs in the TOA scan are orthogonal with each other in Figure 1b. We use *ατ* to denote the LAR of an SA scan and investigate image reconstruction from data collected over SAs of LARs *ατ* = 14◦, 20◦, 30◦, 45◦, 60◦, 90◦, 120◦, 150◦, and 180◦. For an SA of LAR *ατ*, we also consider a TOA scan with two arcs of equal LARs satisfying *<sup>α</sup>*<sup>1</sup> = *<sup>α</sup>*<sup>2</sup> = 0.5*ατ*. (The work can readily be generalized to a TOA scan with two arcs of different LARs [23]).

Dual-energy data are generated from a chest phantom and a suitcase phantom in Figure 2 with two different fan-beam geometries used in the numerical study: for the chest phantom, the source-to-rotation distance (SRD) and source-to-detector distance (SDD) are 100 cm and 150 cm, with a linear detector of 70 cm comprising 896 bins, whereas for the suitcase phantom, the SRD and SDD are 100 cm and 150 cm, with a linear detector of 32 cm including 512 bins. The imaged objects are assumed to be completely within the field-ofview of the scan configurations, resulting in no truncation. In the studies involving both phantoms, the angular interval is fixed at 0.25◦ between two adjacent views. Dual-energy data are also collected over two full rotations, or the FAR of 360◦, and images reconstructed from FAR data may be used as references in the work.

**Figure 1.** Schematics of SA (**a**) and TOA (**b**) scans of LARs in fan-beam DECT. The SA scan includes a pair of completely overlapping arcs of LAR *ατ*, and the *x*-axis intersects with the middle point of the two arcs, whereas the TOA scan includes two pairs of completely overlapping arcs of LARs *α*<sup>1</sup> and *α*2, and the *x*- and *y*-axis intersect with the middle points of the two pairs of arcs. For each pair of the completely overlapping arcs in the SA and TOA scans, low- and high-kVp data are collected over one of the paired arcs. In this work, we consider *<sup>α</sup>*<sup>1</sup> = *<sup>α</sup>*<sup>2</sup> = 0.5*ατ*.

**Figure 2.** Row 1: (**a**) water and (**b**) iodine basis images and (**c**) VMI at 100 keV of the chest phantom; and row 2: (**a**) photoelectric effect (PE) and (**b**) Compton scattering (KN) basis images and (**c**) VMI at 40 keV of the suitcase phantom. Display windows for the chest phantom are [0, 1.2] for the two basis images and [0, 0.22] cm−<sup>1</sup> for the VMI, while those for the suitcase phantom are [0, 0.22] and [0.1, 0.65] cm−1, respectively.

### *2.2. Imaging Model*

In DECT, data are collected at ray *j* with two distinct spectra, referred to as low- and high-kVp spectra, and an imaging model can be expressed as [8]

$$\begin{aligned} \mathcal{g}\_j^L &= -\ln \sum\_m^M q\_{jm}^L \exp\left(-\sum\_i^I a\_{ji} f\_{mi}\right), \\ \mathcal{g}\_j^H &= -\ln \sum\_m^M q\_{jm}^H \exp\left(-\sum\_i^I a\_{ji} f\_{mi}\right). \end{aligned} \tag{1}$$

where *g<sup>L</sup> <sup>j</sup>* and *<sup>g</sup><sup>H</sup> <sup>j</sup>* denote model data of the low- and high-kVp scans; *<sup>q</sup><sup>L</sup> jm* and *<sup>q</sup><sup>H</sup> jm* the lowand high-kVp spectra after normalization (including possibly filtered tube spectra and detector response) at energy bin *m*; *aji* the contribution of image pixel *i* to data of ray *j*; and *fmi* the image value at pixel *i* within energy bin *m* of the monochromatic image, i.e., the linear attenuation coefficient.

In the absence of the basis-decomposition error, *fmi* can be written as the combination of two basis images, i.e.,

$$f\_{\rm mi} = \mu\_{0m} b\_{0i} + \mu\_{1m} b\_{1i} \,\tag{2}$$

where *bki* denotes basis image *k* at pixel *i* and *μkm* the linear attenuation coefficient at energy bin *<sup>m</sup>* for basis material *<sup>k</sup>* (*<sup>k</sup>* = 0 or 1). Image **<sup>f</sup>***m*, with *fmi* as its elements, obtained with Equation (2) is referred to also as the virtual monochromatic image (VMI).

In the work, assuming *μkm*, *q<sup>L</sup> jm*, and *<sup>q</sup><sup>H</sup> jm* are known, the two-step method is proposed for accurately reconstructing basis images *bki*, or, equivalently, VMI *fmi*, from data collected over an SA or TOA of LARs in fan-beam DECT.

#### *2.3. Numerical Phantoms Studied*

We consider in the work two phantoms, i.e., the chest phantom [15] and suitcase phantom [6] shown in Figure 2, motivated by their possible implications in medical and security imaging, two distinct DECT imaging applications, and their distinctly different anatomic structures for evaluating algorithm performance. The chest phantom contains four regions of interest (ROIs) 1–4 with iodine-contrast agents at concentrations of 5 mg/mL, 10 mg/mL, 15 mg/mL, and 20 mg/mL, respectively, and other ROIs with mixed materials, such as muscle, lung tissue, and bone; whereas the suitcase phantom includes three ROIs 0–2 of single-element calibration materials, i.e., carbon, aluminum, and calcium, and four more ROIs 3–6 of mixed materials, corresponding to water, ANFO (Ammonium Nitrate and Fuel Oil [11]), teflon, and PVC, respectively.

In DECT, basis images may be selected according to the task considered. For the chest phantom, to estimate iodine concentrations, we select material-based basis images of water and iodine concentration of 20 mg/mL, with the corresponding *μkm*'s obtained from the NIST database [24]. For the suitcase phantom, in order to estimate effective atomic numbers, we select interaction-based basis images of the photoelectric effect (PE) and Compton scattering (KN) with *μkm*'s that are 1/*E*3, where *E* denotes X-ray energy, and obtained with the Klein–Nishina formula [1], respectively. The basis images and VMIs of the chest and suitcase phantoms are formed on image arrays of 200 × 256 and 150 × 256 square pixels of size 0.7 mm, as displayed in rows 1 and 2, respectively, in Figure 2.

#### *2.4. Image Reconstruction Approach*

In an attempt to compensate for the BH effect inherent in *g<sup>L</sup> <sup>j</sup>* and *<sup>g</sup><sup>H</sup> <sup>j</sup>* , we rewrite Equation (1) as

$$\begin{aligned} g\_j^L &= -\ln \sum\_m^M q\_{jm}^L \exp\left(-\mu\_{0m} l\_{0j} - \mu\_{1m} l\_{1j}\right), \\ g\_j^H &= -\ln \sum\_m^M q\_{jm}^H \exp\left(-\mu\_{0m} l\_{0j} - \mu\_{1m} l\_{1j}\right). \end{aligned} \tag{3}$$

where *lkj* = <sup>∑</sup>*<sup>I</sup> <sup>i</sup> ajibki*, *<sup>k</sup>* <sup>=</sup> 0 or 1, denotes the sinogram of basis image *<sup>k</sup>*, also referred to as basis data, which is independent of energy *m*. Therefore, applying the DDD algorithm [10] to Equation (3), we can obtain basis sinograms *lkj* from knowledge of *g<sup>L</sup> <sup>j</sup>* and *<sup>g</sup><sup>H</sup> <sup>j</sup>* for each ray *j*. It has been shown empirically [25] that the DDD algorithm can recover accurately basis sinograms from *g<sup>L</sup> <sup>j</sup>* and *<sup>g</sup><sup>H</sup> <sup>j</sup>* . Using existing algorithms such as the FBP algorithm, one can reconstruct readily accurate basis images from full knowledge of basis sinograms *lkj* in a FAR or short scan. In the work, because knowledge of *lkj* can be available only over a SA or TOA of LARs, the FBP algorithm yield basis images of significant artifacts. We thus develop and tailor the DTV algorithm to reconstruct basis images with minimized LAR artifacts from knowledge of *lkj*'s available only over a SA or TOA of LARs.

Using vectors **<sup>b</sup>***<sup>k</sup>* and **<sup>L</sup>***<sup>k</sup>* (*<sup>k</sup>* = 0 or 1) of sizes *<sup>I</sup>* and *<sup>J</sup>*, respectively, to denote basis images and their sinograms with elements *bki* and *lkj* in concatenated forms, we formulate the reconstruction problem of basis images from their sinograms as a convex optimization problem

$$\begin{aligned} \mathbf{b}\_k^\star &= \underset{\mathbf{b}\_k}{\text{argmin }} \frac{1}{2} \parallel \mathbf{L}\_k - \mathcal{A} \mathbf{b}\_k \parallel\_2^2\\ \text{s.t. } ||\mathcal{D}\_x \mathbf{b}\_k||\_1 &\le t\_{kx} \cdot ||\mathcal{D}\_y \mathbf{b}\_k||\_1 \le t\_{ky} \text{ and } b\_{ki} \ge 0, \end{aligned} \tag{4}$$

where matrix A of size *J* × *I* denotes the discrete fan-beam X-ray transform with element *aji*; · <sup>2</sup> the -2-norm of a vector; and ||D*x***b***k*||<sup>1</sup> and ||D*y***b***k*||<sup>1</sup> are the image's directional total variations (DTVs) [5] of the basis image **b***<sup>k</sup>* along the *x*- and *y*-axis, respectively.

The DTV algorithm used to reconstruct basis images from knowledge of the basis sinograms in DECT through solving Equation (4) shares the same general structure as that of the algorithm for image reconstruction from LAR data in conventional SECT [5]. The pseudo-code is thus summarized in Appendix A for clarity.

#### *2.5. Visual Inspection and Quantitative Analysis of Images*

As VMIs are of visualization interest in DECT, we first obtain VMIs at energy levels of interest from basis images reconstructed by using Equation (2) and then visually inspect LAR artifacts in the VMIs. Additionally, two quantitative metrics, normalized root-mean-square error (nRMSE) and Pearson correlation coefficient (PCC) [5,26,27] are calculated. Metric nRMSE evaluates quantitative difference, while metric PCC assesses visual correlation, between a VMI obtained from LAR data and a reference image obtained from FAR data. In particular, higher PCCs suggest better visual correlation between the VMI and its reference image. The VMI and its reference are identical when PCC → 1 and nRMSE → 0.

In the chest phantom study, we seek to estimate iodine-contrast concentration within ROIs 1–4 shown in the basis images in row 1 of Figure 2. Using the estimated basis image of 20-mg/mL iodine-contrast agent, we estimate the concentration of iodine-contrast agent within ROIs 1–4 with a linear fitting [6]. Constants in the linear relationship are determined by using pixel values within iodine-contrast ROIs 1–4 in the reference image of the chest phantom obtained from the FAR data by use of the two-step method, and fitting into the corresponding known concentrations. In the work, the calibrated slope and intercept of the linear fitting were computed as 19.3 mg/mL and −0.0074 mg/mL. In general, the linear fitting, as compared to the default setting of 20 and 0 as slope and intercept, yields more accurate estimation of the iodine concentration, because the mean pixel values within ROI 0 in the 20-mg/mL iodine basis image could be non-zero. This occurs as a result of the incomplete basis set in the material decomposition model by using 2 materials. On the other hand, in the study involving the suitcase phantom, we seek to estimate the effective atomic number of materials [11]. As the basis images are estimated as PE and KN components, their ratios are used in an affine transform with the effective atomic number in the log-log domain [6]. The effective atomic numbers are then computed for ROIs 3–6 of the suitcase phantom, as shown in row 2 of Figure 2. Constants in the affine transformation are determined by using the pixel values within single-element ROIs 0–2 in the reference image of the suitcase phantom obtained from the FAR data by use of the two-step method, and fitting into the corresponding known atomic numbers.

#### **3. Results**

#### *3.1. Numerical Study Design and Data Generation*

In our numerical studies with noiseless and noisy LAR data, the TASMIC model [28] was used for generating filtered tube spectra of given low- and high-kVps. Taking into account the detector's energy-integrating response, we then obtain *q<sup>L</sup> jm* and *<sup>q</sup><sup>H</sup> jm* by multiplying the filtered tube spectra with corresponding X-ray energies *E*. For both phantoms, the low- and high-kVp spectra are set at 80 and 140 kVp, with a 5-mm Al filter used in both.

For each of the chest or suitcase phantom in an SA or TOA scan described in Section 2.1 above, basis sinograms *lkj* are first generated from basis images shown in Figure 2, and noiseless low- and high-kVp data *g<sup>L</sup> <sup>j</sup>* and *<sup>g</sup><sup>H</sup> <sup>j</sup>* can be generated subsequently by use of Equation (3) with *lkj*, and knowledge of *μkm*, *q<sup>L</sup> jm*, and *<sup>q</sup><sup>H</sup> jm* determined. The aims of the noiseless data study are (1) to verify that the two-step method, including the DDD and DTV algorithms, can recover numerically accurate basis images and VMIs first from FAR-scan data and (2) to use the two-step method verified to explore empirically its performance upper bound , i.e., the performance in the best case scenario without any inconsistencies, such as noise and decomposition error, in the data, as a function of LARs for yielding accurate reconstruction of VMIs and physical quantity estimation in DECT with LAR scans.

**Table 1.** NEQs per detector bin in air scans of either the low- or high-kVp scans for the chest and suitcase phantoms with LARs ranging from 14◦ to 180◦, as well as with the FAR of 360◦.


Using noiseless data as the means of the Poisson noise model, we obtain noisy data containing Poisson noise. For both chest and suitcase phantoms, Table 1 shows the noiseequivalent quanta's (NEQs) of each detector bin for the SA or TOA scans studied, which are determined such that the means in SA or TOA scans studied have a fixed total number of quanta of <sup>∼</sup>6.9 <sup>×</sup> 109 in an air scan, amounting to 75% of that in a FAR scan with 360 projection views, 512 detector bins, and <sup>∼</sup><sup>5</sup> <sup>×</sup> 104 NEQs per detector bin [23]. The purpose of the noisy data study is to yield some preliminary insights into the reconstruction robustness of the two-step method. Clearly, its reconstruction accuracy depends not only on the LAR extent but also on the characteristics and level of data noise. No additional data or image processing is applied in the study with noisy data, although such processing may improve the quality of VMI visualization and physical quantity estimation.

Constraint parameters *tkx* and *tky* have an impact on image reconstruction by defining the feasible solution set of Equation (4). In the study below with consistent noiseless data, the DTV values of the phantom basis images in Figure 2 are selected as the values of parameters *tkx* and *tky*, in order to form the tightest feasible solution set that still contain the desired solution (i.e., the truth basis images in this case). In the study with noisy data, the values of parameters *tkx* and *tky* are selected in terms of visual evaluation of reconstructed VMIs with minimum artifacts [5,6]. In general, parameter selection is accomplished through surveying the parameter space within relevant ranges and optimizing a well-defined imagequality metric, e.g., image visualization for artifact reduction or quantitative estimation of iodine-contrast concentration, for studies with inconsistent data, including those with real data where the truth images are not available. In our experience, *tkx* and *tky* selected in the noisy data studies are generally smaller than those in the corresponding noiseless data studies. In the work, the *tkx* values selected are in general smaller than *tky* in the SA scans, so as to suppress horizontal streaks along the *y*-axis, while both *tkx* and *tky* selected in the TOA scans are slightly larger than those in the SA scans, as the improved conditioning of the system matrix leads to fewer artifacts overall. Basis images are also reconstructed from **L***<sup>k</sup>* estimated by using the FBP algorithm with a Hanning kernel and a cutoff frequency at 0.5, which are then combined into the VMIs with Equation (2). The FBP algorithm is used only for demonstrating the LAR artifacts associated with the phantoms and data conditions in the work.
