Simultaneous Image Reconstruction and Element Decomposition for Iodine Contrast Agent Visualization in Multienergy Element-Resolved Cone Beam CT

Iodine contrast agent is widely used in liver cancer radiotherapy at CT simulation stage to enhance detectability of tumor. However, its application in cone beam CT (CBCT) for image guidance before treatment delivery is still limited because of poor image quality and excessive dose of contrast agent during multiple treatment fractions. We previously developed a multienergy element-resolved (MEER) CBCT framework that included x-ray projection data acquisition on a conventional CBCT platform in a kVp-switching model and a dictionary-based image reconstruction algorithm that simultaneously reconstructed x-ray attenuation images at each kilovoltage peak (kVp), an electron density image, and elemental composition images. In this study, we investigated feasibility using MEER-CBCT for low-concentration iodine contrast agent visualization. We performed simulation and experimental studies using a phantom with inserts containing water and different concentrations of iodine solution and the MEER-CBCT scan with 600 projections in a full gantry rotation, in which the kVp level sequentially changed among 80, 100, and 120 kVps. We included iodine material in the dictionary of the reconstruction algorithm. We analyzed iodine detectability as quantified by contrast-to-noise ratio (CNR) and compared results with those of CBCT images reconstructed by the standard filter back projection (FBP) method with 600 projections. MEER-CBCT achieved similar contrast enhancement as FBP method but significantly higher CNR. At 2.5% iodine solution concentration, FBP method achieved 170 HU enhancement and CNR of 2.0, considered the standard CNR for successful tumor visualization. MEER-CBCT achieved the same CNR but at ~6.3 times lower iodine concentration of 0.4%.


INTRODUCTION
Image guidance is a critical component of cancer radiation therapy. By imaging the patient anatomy at the treatment position, image guidance allows accurate positioning of the tumorous target against the therapeutic radiation beam, thereby increasing dose to the target and sparing dose to nearby healthy organs. Clinical advantages of image guidance in radiotherapy have been previously demonstrated by a number of studies (1).
Cone beam CT (CBCT) installed on a medical linear accelerator (LINAC) is currently the most widely used tool for imaging guidance in radiation therapy (2). Yet, one challenging context of tumor visualization is liver cancer. Due to low softtissue contrast in x-ray CT and similar x-ray attenuation properties between liver tumor and normal liver, it is difficult to directly visualize liver tumor in CBCT, leading to tumor targeting uncertainty. This problem is further exacerbated by other issues such as imaging artifacts due to respiration-induced organ motion during image acquisition. Intravenous contrast enhancement (IVC) is routinely used in treatment planning CT of liver cancer radiotherapy because of improved tumor contrast. However, it cannot be reliably employed on standard CBCT due to poor image quality. A previous study found that IVC-CBCT was effective in only 1/4 of patients with tumor size larger than 120 cm 3 under breath-hold scans. Small tumors or free-breathing CBCT do not show contrast enhancement (3). The dose of repeated injection of the iodine contrast agent over multiple treatment fractions further leads to the concern of toxicity. Hence, there is a strong desire to improve CBCT technology to support image guidance in liver cancer radiotherapy. This is particularly needed in adaptive radiotherapy to enable accurate tumor delineation at each fraction (4,5).
One potential approach to increase sensitivity of CBCT imaging to imaging contrast agent for tumor visualization is to take advantage of energy dimension of the x-ray imaging. As the energy dependence of x-ray attenuation property of iodine is different from other materials in a human body, such as tissue and bone, it may be possible to differentiate iodine contrast from other materials by utilizing information provided by the energy dimension. In fact, extensive studies have been conducted on the CT platform to realize dual-energy or multienergy CT function and to differentiate materials, with iodine imaging being one of the major applications (6)(7)(8)(9). On the CBCT platform side, Zbijewski et al. (10) studied the accuracy of material classification in dual-energy CBCT under different reconstruction algorithms using a table-top system. Lee et al. (11) realized the single-scan dual-energy CBCT function using a multislit filter installed between the x-ray source and the scanned object. The filtered and unfiltered x-ray beams generated projection data at two different energy levels. A novel reconstruction algorithm was developed to employ the joint sparsity between the low-and high-energy CT images. In another study, Li et al. (12) developed a reconstruction algorithm that utilized spatial and spectral correlation among images to improve image quality. Lately, Cassetta et al. (13) achieved fast-kilovoltage peak (kVp) switching function for dual-energy CBCT using the on-board imager of a commercial LINAC and successfully produced virtual monoenergetic and relative electron density images.
In our previous study, Shen et al. (14) successfully controlled the CBCT platform on a commonly used LINAC to realize multienergy CBCT data acquisition via the kVp switching scheme, i.e., taking x-ray projections with the kVp level cycling through different levels among projections. To address the undersampling problem caused by the kVp switching, they also developed a multi-energy element-resolved (MEER) CBCT framework to simultaneously reconstruct CBCT images at different kVp levels, as well as electron density map and the maps of a few major chemical elements. A physics model that correlates the CBCT images at different kVp levels, the electron density image, and the element images was built into the reconstruction algorithm, which served as a strong constraint on the solution to ensure their quality. In this paper, we report our recent developments to investigate the feasibility of using the MEER-CBCT framework for iodine contrast agent identification. The contributions of this study are twofold. First, on the algorithm side, we extended the previous MEER-CBCT reconstruction algorithm Shen et al. (14) to include iodinerelated materials in the dictionary of the algorithm, which allowed us to handle the image reconstruction problem in the presence of the iodine contrast agent. We further converted the nonconvex optimization model in our previous study (14) into a convex form, making it easier and more efficient to solve the problem. Second, we also performed comprehensive simulation and experimental studies to quantitatively determine the lowest possible contrast agent concentration level identifiable in MEER-CBCT images to demonstrate the potential feasibility of this method for iodine contrast imaging.

Image Reconstruction Model
The image reconstruction model generally followed that in our previous study (14). In this study, we converted the model into a convex form to make it numerically more tractable. We also included materials with iodine in the dictionary to handle the reconstruction problem in the presence of the contrast agent.
The general form of MEER-CBCT model can be formulated as follows: min F,r≥0,l≥0 The first term R(F) is established for ME-CBCT reconstruction problem, where F = [f 1 , f 2 , … , f N ]∈ R M × N denotes black CBCT images at N=3 different energy channels. f i ∈ R M represents an image of the x-ray attenuation coefficient at the ith energy channel. M is the number of pixels in the image.
The second term G(F, r, l) characterizes the relationship among F, the image of relative electron density (rED) to water r ∈ R M×M , a diagonal matrix with each diagonal element corresponding to rED of each voxel, and images of elemental compositions (EC) l ∈ R M×D , where D indicates the number of elements of interest. In this study, we focused on three elements (D=3), i.e., hydrogen (H), oxygen (O), and iodine (I). The first two are major elements constructing majority of human nonbone tissues, and I is the effective element of the iodine contrast agent. We did not consider bony tissues, because this study focused on iodine contrast identification and it is expected that iodine contrast agent mainly exists in soft tissues. b is the parameter balancing the contributions from the two terms. Note that r ≥ 0 and l ≥ 0 were naturally requested. In addition, let 1 D and 1 M denote column vectors of length D and M with all entries being in unity, l1 D = 1 M was required, since the compositions of all elements in each voxel should add up to the unity.

Image Reconstruction
We considered a tight-frame-based image reconstruction model (15,16) in this study. As such, the detailed expression of R(F) can be given as where P indicates the x-ray projection matrix characterizing the CBCT data acquisition process, while B represents the acquired CBCT projection data. Note that in a kVp switching data acquisition, images at different kVp levels were projected to different angles. This projection angle information was implicitly contained in the projection matrix P. In this model, a 1 is the regularization parameter, and W is the tight frame (TF) operator. ||s|| F and ||s|| 1 denote the matrix Frobenius norm and l 1 norm, respectively. The first term in Eq. (2) enforced fidelity between the reconstructed CBCT images and acquired projections. The second term encouraged sparsity in the TF-transformed images to help in removing noise and artifacts while preserving edges in F.

Material Decomposition
In Eq. (1), G(F, r, l) was established to relate the x-ray attenuation coefficients with rED and EC. This was achieved based on an empirical model (17): Here, k PE i , k R i and k C i are parameters characterizing contributions of photoelectric effect, Rayleigh scattering, and Compton scattering, respectively, to the x-ray attenuation at the ith energy level. These parameters are dependent on the specific CBCT scanner. For a given scanner, they can be obtained via a calibration process (18). In addition, l was encoded inz and ẑ as z = (lz 3:62 ) 1 3:62 and ẑ = (lz 1:86 ) gives the atomic numbers for D different elements.
Through a simple derivation, the forward model in Eq. (3) can be rewritten into a matrix equation form for multiple energy levels, i.e., We further assumed that EC of each voxel can be sparsely represented over a dictionary consisting of EC of different tissues. Hence, we expressed EC as l ≈ VL, where V ∈ R M × E gives the dictionary coefficients with each row being a sparse vector specifying the contribution of each dictionary material. This leads to a model F = r(VLK + K C M ). With this equality connecting EC, rED with x-ray attenuation image, it would be straightforward to define an objective function as the difference between the two sides of this equation, jjF − r(VLK + K C M )jj 2 F , and minimize it in the reconstruction problem. The requirement that the representation over L is sparse can be achieved by minimizing an objective function term of ||V|| 0 , where ||s|| 0 denotes the l 0 norm.
However, one obvious limitation of this approach is the nonconvex form of the objective function due to the product of r and v, as well as the l 0 norm. To circumvent the problem, we combined r and v and defined a new variable X = rṡV. Since it was required that l1 D = 1 M , it followed that V1 E = VL1 D = 1 M , where 1 E is the column vector of length E with all elements being 1. We also relaxed the sparse penalty from the l 0 form to the convex l 1 form. With these changes, the objective function for the material decomposition can be expressed as where x j ∈R M × 1 is the jth row of x and a 2 is the sparsity regularization parameter.

Combined Model
With the image reconstruction and material decomposition terms defined in Eq. (2) and Eq. (5), the complete convex model of MEER-CBCT reconstructoin can be written as: After we solve this problem to obtain x, we can easily compute r and v as

Numerical Algorithm
The model in Eq. (6) is convex and can be solved by the Alternating Direction Method of Multipliers (ADMM) (19). Similar to the iterative scheme proposed in (14), we incorporated ADMM to split the optimization problem into several subproblems and solved F and x in subproblems in an alternating fashion. Specifically, after introducing auxiliary variable, the augmented Lagrangian function was where , h 1 is the Lagrange multiplier, and m 1 is fixed positive parameter. The iterative scheme became where k indexes the iteration steps. The objective function of the F-subproblem was in a form of a summation of three least squares terms. Hence, it had a closedform solution In practice, we solved this using conjugate gradient algorithm instead of directly computing the matrix inverse, which was computationally challenging.
The U-subproblem was a soft shrinkage problem and can be solved for each pixel independently. The closed form solution is As for the subproblem of x, the objective function is a least square term with the l 1 regularization and nonnegative constraint. There is no closed-form solution, and we have to solve it in an inner iteration scheme. Similar to Eq. (6), we used ADMM to solve this X-subproblem by introducing another auxiliary variable Y. The corresponding augmented Lagrangian function now becames where h 2 is the Lagrange multiplier and m 2 is fixed positive parameter. Here, the iteration scheme was to alternatively update three variables: 2 ) where the superscript p indexes the inner loop. Denote , the closed-form solution of X and Y were The convergence of the proposed algorithm is guaranteed (19). The iterative process of the algorithm is summarized in Algorithm 1.

Dictionary
In our previous study (14), the dictionary was constructed with ECs of 71 human tissues of a reference human listed in a previous publication (20). In this study, to handle the image reconstruction problem in the presence of iodine element, we expanded the dictionary by adding one more tissue with 100% iodine element. Note that using a linear combination of this iodine tissue and other tissues in the dictionary, we can express materials with any iodine concentrations. Additionally, we removed Ca component from the dictionary. In clinical practice, iodine contrast agent only exists in soft tissue material, where the amount of Ca component is expected to be negligibly small.

Evaluations
We performed both simulation studies and experimental validations using a phantom with inserts containing water and iodine solution (175 mgI/ml) of different concentrations, see Figure 1. In the experimental studies, the physical phantom was made of acrylic. Its diameter was 20 cm. Nine inserts were placed inside the phantom with one at the phantom center and eight evenly placed on the periphery with their centers on a circle of a radius of 7.5 cm. Each insert was a 50-ml plastic lab tube with a diameter of 3 cm. The solution for the insert at the phantom center was pure water. The other eight inserts contained iodine solution with concentrations of 0.1%, 0.5%, 1.0%, 1.5%, 2.0%, 3.0%, 5.0%, and 10% in weight (%w/w). The hardware platform used in this study was the on-board imaging system of a Varian TrueBeam LINAC (Varian Medical System, Palo Alto, CA, USA). The source-to-isocenter distance was 100 cm, and the source-to-detector distance was 150 cm. We programmed the CBCT system to implement the kVp switching scanning protocol using its developer mode controlled by a customized "xml" file. We scanned the phantom to acquire x-ray projections. In a CBCT scan with 600 projections evenly distributed in a full gantry rotation, the energy levels were set to sequentially cycle through 80, 100, and 120 kVp. To ensure future clinical translation of the develop method, we employed xray beams in the 80-120-kVp range that are commonly available in clinical CT scanners. We also attempted to spread the kVp values to allow the largest possible spectrum separation, which is beneficial for material decomposition. These considerations led to the selected 80-, 100-, and 120-kVp beams. The tube currents for these three kVp levels were 1.4, 0.8, and 0.5 mAs, respectively. These mAs levels were empirically chosen. A relatively lower mAs was used for the channel with a higher kVp level to make the noise level approximately similar among different energy channels.
In the simulation studies, we constructed a digital phantom with the same dimension as the physical phantom. However, we replaced the background acrylic material with pure water for simplicity. The digital phantom was then voxelized with a voxel size of 1 mm 3 . We defined the material composition and density of all voxels based on known phantom information. x-ray projections of the digital phantom was then calculated following the same kVp switching scheme as in the actual experiment using our in-house developed Monte Carlo simulation tool (21). All other settings in the simulation matched those of the experiment.
In both simulation and experimental studies, we first reconstructed the MEER-CBCT images using Algorithm 1. To benchmark our method, we compared the reconstruction results against those of the clinical standard Filtered Back Projection (FBP) reconstruction method and those of iterative reconstruction method using single energy data. For FBP reconstruction, projections at 80, 100, and 120 kVps were simulated/acquired individually with 600 projections covering a whole gantry rotation to allow sufficient number of projections to reduce streak artifacts. This means that the comparison actually favored the FBP reconstruction case, as in the MEER-CBCT case, each energy only had one-third number of x-ray projections. For iterative reconstruction with single energy, the setup was the same as MEER-CBCT. The algorithm solved the problem in Eq. (1) with b = 0 to determine images F. This setting essentially treated the three energy channels independently and ignored the interchannel relationship expressed in Eq. (4).
To quantitatively evaluate the results, we considered two metrics. The first one was contrast enhancement defined as x roi − x bg , where x roi is the mean value of the image intensity of a region of interest (ROI) and x bg is the mean value of the image intensity of the background. The ROIs were selected as square regions within the eight circular inserts at the periphery region, and the background was the square region in the insert located at the center of the phantom. In the experimental study, the ROI selection was carefully performed to avoid bubbles in each insert. We used this metric to measure contrast enhancement caused by the iodine contrast agent. The second metric was contrast-to-noise ratio (CNR), a key quantity characterizing the detectability of an object (22). CNR was defined as  where s roi and s bg are the standard deviations of image intensity at the ROI and the background, respectively. Additionally, since we know the ground truth iodine concentrations, we also evaluated errors of iodine concentrations derived by MEER-CBCT. All the numerical computations analyzing results were conducted on a desktop with CPU (Intel i7-6700, 3.4 GHz) and MATLAB 9.2 (R2017a).

Convergence and Parameter Sensitivity Analysis
Before presenting the image reconstruction results, we first show the convergence property of the numerical algorithm and results on sensitivity analysis of parameters in the algorithm. We empirically demonstrated the convergence of the proposed ADMM algorithm in Figure 2. Specifically, we examined the objective value in Eq. (6). We also studied the relative change in the restored images F and the variable x that was introduced to convexify the optimization problem between two successive iteration steps during the iteration process. Here, we only considered the simulation study. All the three quantities showed monotonically decaying trends. The objective function saturated at the end, indicating that the iteration reached convergence. We stopped the iteration at the step number 40, where the objective function value did not decrease significantly any further, and the relative changes of x and F were less than 10 -2 and 10 -3 , respectively.
There are multiple parameters in Algorithm 1, whose values are expected to affect the final results. Hence, it is important to study the sensitivity of the results to these parameters and select their optimal values. There are five parameters in Algorithm 1, b, a 1 , a 2 , m 1 , and m 2 . The two parameters m 1 and m 2 are expected to only affect algorithm convergence rate, but not the results. Therefore, we studied the sensitivity of parameters b, a 1 , and a 2 . It would be computationally challenging to scan the entire range of these parameters. Hence, we first sampled a few combinations of these parameters in the possible value range and gradually adjusted the parameters in a trial-and-error way to determine the optimal parameter set that yielded the highest CNR of different inserts in the reconstructed image with 80 kVp. After that, we fixed two parameters at their optimal values and studied the dependence of CNR as a function of the third parameter. Similar to the convergence study, the simulation case was used here. The results are shown in Figure 3. Based on this study, the optimal parameter values were b = 50 and a 1 = 1.0. The CNR was found to be not sensitive to a 2 , and we set it to a 2 = 0.001 in this study. These parameter values were used in subsequent image reconstruction studies.

Reconstruction Results
Figures 4 and 5 present the reconstructed MEER-CBCT images in the simulation and the experimental studies, respectively. The first and the second rows in each figure are images of the phantom images reconstructed by the FBP algorithm and the MEER-CBCT algorithm, and columns are for different kVps. Comparing the results using the FBP reconstruction method, the images produced by the MEER-CBCT method achieved improved image quality, as indicated by visually reduced noise level, which can be ascribed to the inclusion of regularization terms in spatial and energy dimension in the reconstruction process.
For quantitative comparison, we first present the results of contrast enhancement in Figure 6. As expected, the contrast enhancement increased approximately linearly with respect to the iodine agent concentration. The contrast enhancement became higher for lower kVp levels due to increased photoelectric interactions at low-energy range. The enhancement levels in the images reconstructed by the FBP algorithm and by the MEER-CBCT algorithm generally agreed with each other, indicating that the use of regularization terms in our method did not suppress image contrast. Figure 7 presents the results of CNRs. As expected, CNRs increased with iodine agent concentration and reduced kVp levels. In both simulation and experimental studies, for a given concentration of the iodine solution, the corresponding CNR obtained by our method was significantly higher than that of the FBP reconstruction method. For instance, the CNR of MEER-CBCT was approximately 6.5 times of that of the FBP method at 2.5% iodine contrast concentration in the experimental study (13.0 and 2.0, respectively). This can be ascribed to the use of image domain regularization, as well as the correlation among images at different kVps that was made possible by Eq. (4). We considered CNR = 2.0 as the threshold of tumor detectability. Based on the experimental study [ Figure 7 (right)], the case with 2.5% of iodine contrast solution concentration yielded this level of CNR (averaged over all energy channels) for the conventional FBP reconstruction method. At this concentration level, Figure 6 indicates that the averaged contrast enhancement over all energy channels was about~170 HU. In clinical practice, the reported contrast between liver tumor and normal liver varies in the range of 50 HU to~300 HU depending on specific protocols of CT scan and contrast injection (23)(24)(25). The enhancement at this 2.5% iodine contrast case fell in this range. Compared with the FBP reconstruction method, MEER-CBCT achieved the same CNR level of 2.0 at~0.4% of iodine solution concentration, about 6.3 times reduction of concentration. Figure 8 presents comparisons of contrast enhancement results and CNRs between MEER-CBCT images and those reconstructed by iterative reconstruction algorithm for each energy channel independently in the simulation case. It was observed that the two methods maintained the same level of contrast enhancement. MEER-CBCT was able to improve CNR by approximately a factor of 2 because of the incorporation of interenergy relationship.
One advantage of MEER-CBCT algorithm is the capability of resolving iodine contrast distribution. Figure 9A, B presents the computed iodine concentration image in the simulation and the experimental studies, respectively. We plotted the iodine concentration in each insert over the known ground truth value in Figure 9C. The results derived by MEER-CBCT agreed well with ground truth values in general. The deviation from ground truth started to appear at low iodine concentration cases. Note the plot is on a log-log scale. The mean relative errors for the simulation and the experimental studies were 15% and 20%.

DISCUSSION
One important aspect for clinical application of CBCT for radiotherapy image guidance is x-ray radiation dose, as  excessive radiation dose can increase risk of secondary cancer. In the MEER-CBCT approach, we acquired x-ray projection data in a single rotation via the kVp switching approach. The radiation dose is hence expected to be approximately the average of the radiation dose of CBCT scans of each individual kVp scan, which is considered acceptable. One contribution of this study is convexifying the reconstruction model presented in our previous work (14). The advantages of this approach included a unique solution to the optimization problem, and more importantly, a numerical algorithm that can solve the problem more efficiently. In the simulation study, our algorithm was 3.5 times faster than the one used to solve the nonconvex model in (14) with the same performance. Meanwhile, because of the convex nature, our model does not rely on the choice of initial guess, which is a favorable feature for practical applications.
Despite the exciting results in terms of achieving the same level of CNR as the FBP algorithm with substantially reduced iodine dose, the current study has a few limitations. First, the study was performed using a relatively small size phantom, where the impact of x-ray scatter was not significant. It is well known that x-ray scatter is a major concern affecting quantitative accuracy of CBCT imaging due to the large x-ray illumination field and detector size (26). In the liver site, the large body size would increase scatter component and reduce primary component of the x-ray beam at the detector, making the impact of scatter more profound than that of the phantom case. Hence, future challenges include proper removal of scatter before CBCT reconstruction. Over the years, advanced hardware-based or computation-based scatter correction methods have been successfully developed (26)(27)(28). These methods could be used to preprocess the x-ray FIGURE 6 | Comparison of contrast enhancements between our method and FBP reconstruction method in simulation study (left) and experimental study (right). Result of 10% iodine concentration not displayed to focus on relevant data range. projection data before using our algorithm for MEER-CBCT reconstruction.
Second, many factors affecting contrast agent visualization in CBCT images have not been considered in this initial study, which hence pose challenges to translate our method to clinical practice. The first factor is kinetic behavior of contrast enhancement of liver tumor. In fact, the liver tumor contrast enhancement shows a complex kinetic behavior occurring during a time interval with a length comparable with the CBCT data acquisition time (23). The CBCT acquisition has to be performed at a proper time after contrast injection, and the data acquisition time would average the contrast enhancement  level. The second factor is organ motion. Liver is subject to respiratory motion. This motion during the CBCT scan, if unaddressed, would further blur the resulting images and diminish the contrast enhancement in the reconstructed images. We are in the process of developing a liver phantom with a realistic contrast enhancement mechanism and motion to study the impact of these factors. Novel image reconstruction techniques, such as reconstruction algorithms with temporal dimension included (29) or using deep learning (30), may be potentially employed to overcome these challenges.

CONCLUSION
To improve the detectability of iodine contrast agent in CBCT for image guidance of liver cancer radiotherapy, we developed a MEER-CBCT framework that acquired x-ray projections in a kVp switching scan on a conventional CBCT platform of a LINAC. MEER-CBCT image reconstruction method simultaneously reconstructed x-ray attenuation images at all kVp levels, the image of rED and images of EC. The composition of each voxel was subject to a constraint of a sparse representation of materials in a dictionary containing typical human tissues and iodine. We converted the nonlinear formalism of MEER-CBCT reconstruction problem to a linear form to ease the burden solving this problem. In both simulation and experimental studies, MEER-CBCT achieved similar contrast enhancement as the clinical standard FBP reconstruction method but significantly higher CNR. At 2.5% iodine solution concentration, FBP method achieved~170 HU enhancement and CNR of~2, considered the acceptable CNR for successful liver tumor visualization. MEER-CBCT yielded the same CNR but at~6.3 times lower iodine concentration of 0.4%.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors upon request.

AUTHOR CONTRIBUTIONS
XJ, CS, and MY contributed to the study concepts and study design. CW and HJ contributed to data acquisition and image reconstructions. CW, HJ, and CS contributed to the data analyses and interpretation. All authors contributed to the manuscript preparation and editing. All authors contributed to the article and approved the submitted version.