Non-convex primal-dual algorithm for image reconstruction in spectral CT
Buxin Chen, Zheng Zhang, Dan Xia, Emil Y. Sidky, Xiaochuan Pan
NNon-convex primal-dual algorithm for image reconstruction in spectral CT
Buxin Chen a , Zheng Zhang a , Dan Xia a , Emil Y. Sidky a , Xiaochuan Pan a,b, ∗ a Department of Radiology, The University of Chicago, Chicago, IL 60637, USA. b Department of Radiation and Cellular Oncology, The University of Chicago, Chicago, IL 60637, USA
Abstract
The work seeks to develop an algorithm for image reconstruction by directly inverting the non-linear datamodel in spectral CT. Using the non-linear data model, we formulate the image-reconstruction problemas a non-convex optimization program, and develop a non-convex primal-dual (NCPD) algorithm to solvethe program. We devise multiple convergence conditions and perform verification studies numerically todemonstrate that the NCPD algorithm can solve the non-convex optimization program and under appropri-ate data condition, can invert the non-linear data model. Using the NCPD algorithm, we then reconstructmonochromatic images from simulated and real data of numerical and physical phantoms acquired witha standard, full-scan dual-energy configuration. The result of the reconstruction studies shows that theNCPD algorithm can correct accurately for the non-linear beam-hardening effect. Furthermore, we applythe NCPD algorithm to simulated and real data of the numerical and physical phantoms collected with non-standard, short-scan dual-energy configurations, and obtain monochromatic images comparable to those ofthe standard, full-scan study, thus revealing the potential of the NCPD algorithm for enabling non-standardscanning configurations in spectral CT, where the existing indirect methods are limited.
Keywords:
Spectral CT, photon-counting CT, image reconstruction, non-convex optimization,primal-dual algorithm
1. Introduction
In CT imaging with polychromatic X-rays, thebeam-hardening effect necessarily leads to a non-linear data model (Alvarez and Macovski, 1976).Works exist on investigating indirect methods ordirect algorithms for image reconstruction in dual-energy and spectral CT based on the non-lineardata model. An indirect reconstruction methodfirst converts non-linear data into linear data fromwhich an image is reconstructed subsequently withan algorithm, such as the filtered-back projection(FBP) algorithm, designed for inverting a lineardata model (Alvarez and Macovski, 1976; Roessland Proksa, 2007; Zou and Silver, 2008; Zhanget al., 2014; Sawatzky et al., 2014; Rigie and LaR-iviere, 2015). It generally requires multiple mea-surements made for each of the X-ray paths with afull set of multiple spectra for performing data con-version. The measurement condition required by ∗ Corresponding author
Email address: [email protected] (Xiaochuan Pan) the indirect methods may not be satisfied by non-standard scanning configurations of interest, thuslimiting their application to configurations withoverlapping rays. Conversely, direct algorithms re-construct an image directly by the inversion of thenon-linear data model in a discrete form (Cai et al.,2013; Long and Fessler, 2014; Nakada et al., 2015;Zhao et al., 2015; Barber et al., 2016; Chen et al.,2017; Mory et al., 2018), and they may reconstructimages from data measured only with a subset ofmultiple spectra for some X-ray paths, thus allow-ing for enabling non-standard configurations (Chenet al., 2017, 2018b).The objective of the work is to develop a directalgorithm for image reconstruction by formulatingit as a solution to a non-convex optimization pro-gram involving the non-linear data model in spec-tral CT. While some of the previous works (Cai We simply use terms “data model” and “data model in adiscrete form” interchangeably hereinafter, without causingconfusion.
Preprint submitted to Computerized Medical Imaging and Graphics February 23, 2021 a r X i v : . [ phy s i c s . m e d - ph ] F e b t al., 2013; Long and Fessler, 2014; Nakada et al.,2015; Zhao et al., 2015) have reported the develop-ment of direct algorithms for image reconstructionin spectral CT, it is unclear whether the direct al-gorithms can (mathematically or numerically) ac-curately solve their non-convex optimization pro-grams and invert the non-linear data model. Mean-while, there has been significant advancement insolving convex optimization programs, which maybe exploited for developing algorithms to solve non-convex optimization programs.There has also been an effort in adapting exist-ing algorithms to non-convex optimization, such asthe inertial proximal algorithm for non-convex op-timization (Ochs et al., 2014) and the primal-dual(PD) hybrid gradient method for nonlinear opera-tors (Valkonen, 2014). It remains to be investigatedwhether the algorithms can be applied to invertingthe non-linear data model in spectral CT.We have recently developed iterative algorithmsby adapting a first-order primal-dual (PD) algo-rithm (Chambolle and Pock, 2010; Sidky et al.,2012) to solve convex optimization programs of dif-ferent forms relevant to image reconstruction intomographic imaging systems(Zhang et al., 2016a;Xia et al., 2016; Zhang et al., 2016b). We havealso been investigating to adapt instances of thePD algorithm to non-convex optimization programsfor the application to image reconstruction in spec-tral CT, including our previous work (Barber et al.,2016; Barber and Sidky, 2016) that involved form-ing a convex quadratic local bounding function tothe non-convex data discrepancy term and adaptinginstances of the PD algorithm to provide a descentstep for that quadratic approximation.We continue to seek alternate adaptations of theinstances of the PD algorithm to non-convex op-timization programs that are simpler in an algo-rithmic sense (Chen et al., 2018a). This work firstpresents a derivation of a new instance of the PD al-gorithm, referred to as the convex PD (CPD) algo-rithm and then an intuitive adaptation of the CPDalgorithm, with minimally modified update steps,to a non-convex PD (NCPD) algorithm for a non-convex optimization program. A significant differ-ence between the NCPD algorithm and our previ-ous work (Barber et al., 2016; Barber and Sidky,2016) is that the step-size parameters are fixed andcomputed ahead of time in the former, while theyhave to be calculated at every iteration for the lat-ter.The issue of non-convexity in the work is dealt with by “convexifying” the non-convex optimiza-tion program, using a global linear expansion ofthe non-linear data model about a set of stationarybasis images selected. This strategy has been ap-plied to adapting the ASD-POCS algorithm (Sidkyand Pan, 2008) to a non-convex optimization inspectral CT (Chen et al., 2017, 2018b). We in-vestigate here the strategy for adapting the newCPD instance of the general PD algorithm to theNCPD algorithm for image reconstruction in spec-tral CT. Specifically, the NCPD algorithm is basedupon the CPD algorithm, thus eliminating issuesassociated with the ASD-POCS algorithm: (1) itis only a heuristic algorithm for convex optimiza-tion programs and (2) it can be applied only toconvex optimizations of a specific form with data- (cid:96) -constrained TV-minimization.Following the derivation and demonstration ofthe NCPD algorithm in Sections 2 and 3 for solv-ing the non-convex program and under appropri-ate data condition, for inverting the non-linear datamodel in spectral CT, we then perform studies onreconstruction of the monochromatic images by us-ing the NCPD algorithm from full-scan simulated-and real-data in spectral CT in Sections 4 and 5 toshow accurate correction of the non-linear beam-hardening effect, along with a study in Section 6that demonstrates the potential of the NCPD algo-rithm for enabling non-standard scanning configu-rations in spectral CT.
2. Methods: derivation of the NCPD algo-rithm
In spectral CT, a monochromatic image is of in-terest that represents the spatial distribution of X-ray linear-attenuation coefficient at a particular en-ergy within the subject scanned. Let vector f m ofsize I denote the monochromatic image with ele-ment f mi within voxel i at discrete energy bin m ,where m ∈ , · · · , M , i ∈ , · · · , I , and M and I are total numbers of energy bins and voxels, respec-tively. We use vector b k of size I to denote basis im-age k with element b ki within voxel i , k ∈ , · · · , K ,and K is the total number of basis images consid-ered. An aggregated basis image vector b of size K c = KI can be formed by concatenating b k .We assume that monochromatic image f m can beexpressed as f m ( b ) = (cid:88) k µ mk b k , (1)2here decomposition coefficients µ mk are assumedto be known. Therefore, the reconstruction task of f m is tantamount to that of basis image b .Let g [ M ] sj denote measured data, after air-scannormalization and logarithmic transform, for ray j with spectrum s , where j ∈ , · · · , J s , J s is thetotal number of rays illuminated with spectrum s , s ∈ , · · · , S , and S is the total number of differentspectra involved. Using vector g [ M ] s of size J s todenote measured data for spectrum s , we can formaggregated vector g [ M ] of size J = (cid:80) s J s , com-posed by concatenating g [ M ] s of size J s , and simplyrefer to g [ M ] as measured data. The task of imagereconstruction considered in the work is to deter-mine basis image b from knowledge of measureddata g [ M ] . Using Eq. (1), we can write the non-linear datamodel in a discrete form in spectral CT as (Chenet al., 2017) g [ NL ] sj ( b ) = − ln (cid:88) m q sjm exp (cid:32) − (cid:88) k µ mk (cid:88) i a sji b ki (cid:33) , (2)where g [ NL ] sj ( b ) denotes the non-linear-model datafor ray j with spectrum s ; q sjm normalized X-ray spectrum s at energy bin m for ray j , i.e., (cid:80) m q sjm = 1; a sji the discrete X-ray transformspecified by ray j and voxel i with spectrum s .As in Zhao et al. (2015); Chen et al. (2017), weperform a Taylor-series expansion of Eq. (2) at ¯ b = : g [ NL ] sj ( b ) = (cid:88) k µ sjk (cid:88) i a sji b ki + ∆ g [ NL ] sj ( b ) , (3)where the zeroth-order term is zero at ¯ b = 0, and µ sjk , the spectrum-weighted average attenuationcoefficient independent of energy, and ∆ g [ NL ] sj ( b ),the high-order, non-linear term in the expansion,are given by µ sjk = (cid:88) m q sjm µ mk , (4)∆ g [ NL ] sj ( b ) = g [ NL ] sj ( b ) − (cid:88) k µ sjk (cid:88) i a sji b ki . (5) For simplicity, we refer to q sjm as X-ray spectrum, eventhough it should be understood conceptually as the combi-nation of X-ray-source spectrum and detector-spectrum re-sponse. For spectrum s , we form vectors g [ NL ] s ( b ) and∆ g [ NL ] s ( b ) of size J s , respectively, with elements g [ NL ] sj ( b ) and ∆ g [ NL ] sj ( b ). Subsequently, we can alsoform aggregated vectors g [ NL ] ( b ) and ∆ g [ NL ] ( b ) ofsize J , respectively, by concatenating g [ NL ] s ( b ) and∆ g [ NL ] s ( b ). With the aggregated vectors defined,we rewrite Eq. (2) in a matrix-vector form as g [ NL ] ( b ) = H b + ∆ g [ NL ] ( b ) , (6)where matrix H of size ( J × K c ) is given by H = R A R A · · · R K A R A R A · · · R K A ... ... . . . R S A S R S A S · · · R SK A S , (7) R sk is a diagonal matrix of size ( J s × J s ) with diag-onal element µ sjk , and A s denotes a matrix of size( J s × I ) with element a sji .The objective of the work is to develop the NCPDalgorithm for reconstructing basis image b , andsubsequently monochromatic image f m ( b ), frommeasured data g [ M ] by inversion of the non-lineardata model in Eq. (6). Using the non-linear data model in Eq. (6),we formulate the image reconstruction problem inspectral CT as an optimization program as b (cid:63) = arg min b D g [ NL ] ( b )s . t . || ( |∇ f m (cid:48) ( b ) | ) || ≤ γ & f m (cid:48) ( b ) (cid:23) , (8)where data term D g [ NL ] ( b ) is given by D g [ NL ] ( b ) = || g [ M ] − g [ NL ] ( b ) || / || ( g [ M ] − ∆ g [ NL ] ( b )) − H b || / , (9) || ( |∇ f m (cid:48) ( b ) | ) || the total variation (TV) ofmonochromatic image f m (cid:48) , ∇ a matrix of size 3 I × I for spatial finite difference, (cid:23) a vector-element-wise non-negativity constraint of monochromaticimage, and γ a non-negative parameter.The optimization program in Eq. (8) is similar inform to the ones used in our previous work (Zhanget al., 2016a; Xia et al., 2016; Zhang et al., 2016b),however, the TV and non-negativity constraints areapplied to the monochromatic image, which is a lin-ear combination of the basis images, instead of the3asis images themselves, at a given energy. Suchconstraints are used because the monochromaticimage directly represents physical quantity, i.e., theX-ray linear attenuation, of application interest. Inthe numerical studies of the work, both the TVand non-negativity constraints are enforced on amonochromatic image at an energy of 100 keV.Because the Hessian of data term D g [ NL ] ( b ) canbe shown to have negative eigenvalues, D g [ NL ] ( b )and consequently the optimization program in Eq.(8) are non-convex. Basing on the CPD algorithmdeveloped in Appendix A and Appendix B, wederive below the NCPD algorithm for solving thenon-convex optimization program in Eq. (8). The NCPD algorithm is developed heuristicallybelow because no algorithms exist, to the best ofour knowledge, that solve mathematically exactlythe non-convex program in Eq. (8). In particular,the NCPD algorithm is adapted from the CPD al-gorithm presented in Appendix A and AppendixB for solving a convex program in Eq. (A.2) in-volving the linear data model in Eq. (A.1). It canbe observed that Eqs. (8) and (A.2) are identi-cal in form except for the difference between terms∆ g [ NL ] ( b ) in Eq. (6) and ∆ g c in Eq. (A.1). It isthe similarity in form of the two programs that mo-tivates the adaptation of the NCPD algorithm fromthe mathematically exact and numerically accurateCPD algorithm derived.It can be observed that if ∆ g [ NL ] ( b ) is replacedwith constant term ∆ g c , Eq. (8) becomes convexand can thus be solved by use of the CPD algo-rithm in Appendix A and Appendix B. Motivatedby this observation, we devise below the heuristicNCPD algorithm for numerically solving Eq. (8) byusing ∆ g [ NL ] ( b ( n ) ), computed with b ( n ) at currentiteration in Eq. (5), as an estimate of ∆ g c in line 6of Algorithm 1 in Appendix A, i.e.,6 : p ( n +1) = p ( n ) − σ m (cid:48) [ g [ M ] − ∆ g [ NL ] ( b ( n ) ) −H ¯ b ( n ) ]1 + σ m (cid:48) . (10)We thus obtain the NCPD algorithm withpseudo codes identical to those listed in Al-gorithm 1 in Appendix A except that line 6is now given by Eq. (10) above. The heuristiccomponent in the NCPD algorithm is thus limitedonly to Eq. (10) because the CPD algorithm isbased upon a mathematically rigorous footing. As a part of the NCPD algorithm development,we devise below convergence conditions, which areused for numerically demonstrating that the NCPDalgorithm can solve the non-convex optimization inEq. (8) and that under appropriate data condition,can invert the non-linear data model in Eq. (6).
It should be pointed out that different conver-gence conditions, as discussed below, can be de-signed, depending on whether measured data areconsistent (i.e., g [ M ] = g [ NL ] ( b [truth] )) or incon-sistent (i.e., g [ M ] (cid:54) = g [ NL ] ( b [truth] )) with the datamodel. For measured data that are either consistent orinconsistent with the data model, we first devisegeneral convergence conditions of the NCPD algo-rithm as d (cid:101) D ( n ) g → (cid:101) D ( n )TV → d (cid:101) D ( n ) b → , (11)as n → ∞ , where the dimensionless metrics aregiven by d (cid:101) D ( n ) g = | D g [ NL ] ( b ( n ) ) − D g [ NL ] ( b ( n − ) | / || g [ M ] || (cid:101) D ( n )TV = | || ( |∇ f m (cid:48) ( b ( n ) ) | ) || − γ | /γd (cid:101) D ( n ) b = || b ( n ) − b ( n − || / || b ( n − || . (12)It can be observed that the convergence conditionsin Eq. (11) are similar to those in Eq. (A.19) ofthe CPD algorithm in Appendix A.As shown in Eq. (A.20) in Appendix A, theCPD algorithm also has mathematical convergenceconditions as (cid:103) cPD ( n ) → (cid:101) T ( n ) → (cid:101) S ( n ) → , (13)where dimensionless metrics are given by (cid:103) cPD ( n ) = cPD ( n ) cPD (1) , (cid:101) T ( n ) = T ( n ) T (1) , and (cid:101) S ( n ) = S ( n ) S (1) , (14)with conditional primal-dual (cPD) gapcPD ( n ) (Xia et al., 2016), transversalityT ( n ) (Hiriart-Urruty and Lemar´echal, 2013),and dual residual (or splitting gap) S ( n ) (Goldsteinet al., 2013) defined in Eq. (A.23) in Appendix A.While the metrics can still be computed at eachiteration of the NCPD algorithm, it is unclear4hether, due to the heuristic step in Eq. (10), theconditions in Eq. (13) remain mathematical conver-gence conditions of the NCPD algorithm. Noticingthat the NCPD algorithm becomes the CPD algo-rithm as ∆ g [ NL ] ( b ( n ) ) stabilizes to becoming con-stant, we use them instead as empirical convergenceconditions for the NCPD algorithm. For measured data consistent with the non-lineardata model, i.e., g [ M ] = g [ NL ] ( b [truth] ), we can de-vise an additional convergence condition as (cid:101) D ( n ) g = D g [ NL ] ( b ( n ) ) / || g [ M ] || → , (15)as n → ∞ . Clearly, for a case in which mea-sured data are inconsistent with the data model,i.e., g [ M ] (cid:54) = g [ NL ] ( b [truth] ), this condition cannotbe used because data term D g [ NL ] ( b ( n ) ) is gener-ally non-zero. The convergence conditions discussed above areused for demonstrating that the NCPD algorithmcan numerically accurately solve the non-convex op-timization program in Eq. (8). However, it remainsto be demonstrated that the NCPD algorithm cannumerically accurately invert the non-linear datamodel in Eq. (6). The demonstration is performedunder the conditions that g [ M ] = g [ NL ] ( b ) and that b [truth] is known, as the NCPD algorithm is devel-oped based upon the non-linear data model. There-fore, we devise a unique image-convergence condi-tion directly to check on the accuracy of recovering b [truth] as (cid:101) D ( n ) b = || b ( n ) − b [truth] || / || b [truth] || → , (16)as n → ∞ .In all of the numerical studies in Sections 3-6 be-low, we obtain the results when the respective con-vergence conditions are achieved up to four-bytefloating-point precision and simply state that theconvergence conditions are achieved.
3. Methods: numerical demonstration of theNCPD algorithm’s Convergence
We conduct quantitative studies below numeri-cally to demonstrate that the NCPD algorithm cansolve the non-convex optimization program in Eq. (8) and under appropriate data condition, can in-vert Eq. (6), in terms of achieving the respectiveconvergence conditions.A circular fan-beam scanning configuration isconsidered that has source-to-center-of-detectorand source-to-center-of-rotation distances of 1500mm and 1000 mm, along with a 400-mm linear de-tector consisting of 256 1.56-mm detector bins. Weobtain µ mk from the NIST database (Hubbell andSeltzer, 2004) for water and bone basis materials,and create two X-ray spectra of 80 kVp and 140kVp by using the TASMICS worksheet (Hernandezand Boone, 2014).In the studies, data are generated from a numeri-cal disk phantom with water and bone basis images(i.e., the truth basis images) of size 128 × b [truth] is known, truth-monochromaticimage f [truth] m (cid:48) , and its TV (i.e., truth-TV) value, canthus be calculated from the truth basis image. We numerically demonstrate below that theNCPD algorithm can solve the non-convex opti-mization program in Eq. (8) for data that are eitherconsistent or inconsistent with the non-linear datamodel in Eq. (6). (a) Water truth(b) Water recon. (c) Bone truth (d) Bone recon.Figure 1: (a) Truth water & (c) truth bone basis images, and(b) water & (d) bone basis images reconstructed by use ofthe NCPD algorithm, of the disk phantom. Display window:[0.0, 1.5].
Using the fan-beam scanning configuration andnumerical disk phantom described above in the non-linear data model in Eq. (6), we generate measureddata g [ M ] at 160 views uniformly distributed over2 π , and refer to the data as the full-scan consistentdata, i.e., g [ M ] = g [ NL ] ( b ).5pplying the NCPD algorithm, along with γ thatis the truth-TV value, to the full-scan consistentdata, we compute the convergence conditions inEqs. (11) and (15) from images reconstructed andplot them as functions of iteration n in panels (a)-(c) and (g) of Fig. 2 (solid, dark). It can be observedthat the convergence metrics consistently decrease,revealing that the NCPD algorithm can solve thenon-convex optimization program in Eq. (8) by nu-merically achieving the convergence conditions inEqs. (11) and (15). n dD ( n ) g (a) 10 n D ( n )TV (b) 10 n dD ( n ) b (c)10 n cPD ( n ) (d) 10 n T ( n ) (e) 10 n S ( n ) (f)10 n D ( n ) g (g) 10 n D ( n ) b (h) Figure 2: Convergence metrics (a) d (cid:101) D ( n ) g , (b) (cid:101) D ( n )TV , (c) d (cid:101) D ( n ) b , (g) (cid:101) D ( n ) g , and (h) (cid:101) D ( n ) b , and empirical convergencemetrics (d) (cid:103) cPD ( n ) , (e) (cid:101) T ( n ) , and (f) (cid:101) S ( n ) , of the NCPD al-gorithm (solid, dark), as functions of iteration n , from full-scan consistent data. For comparison, their counterparts forthe CPD algorithm (dashed, red) from Fig. B.12 are alsoover-plotted. We also compute and show in panels (d)-(f) ofFig. 2 (solid, dark) empirical convergence condi-tions in Eq. (13). Considering the NCPD algorithmis adapted from the CPD algorithm, we are inter-ested in comparing their convergence properties. Assuch, we also inlay in Fig. 2 the corresponding con-vergence results (dashed, red) of the CPD algorithmfrom Fig. B.12 in Appendix B. It can be observedthat the convergence property of the NCPD algo-rithm is analogous to that of the CPD algorithm.Therefore, the convergence property of the CPD al-gorithm provides a nice benchmark for that of theNCPD algorithm.
We now show that the NCPD algorithm can alsosolve the non-convex optimization program in Eq.(8) for inconsistent data. We generate full-scan in-consistent data g [ M ] by adding Poisson noise, cor-responding to 10 photons per detector bin in anair scan, to full-scan consistent data g [ NL ] ( b [truth] )above. Therefore, g [ M ] (cid:54) = g [ NL ] ( b ), and Eq. (15)cannot be used as a convergence condition as (cid:101) D ( n ) g does not approach to zero as n → ∞ .Applying the NCPD algorithm, along with γ thatis 98% of the truth-TV value, to the full-scan in-consistent data, we compute the convergence con-ditions in Eq. (11) and plot them as functions ofiteration n in panels (a)-(c) of Fig. 3. We also com-pute and show in panels (d)-(f) of Fig. 3 empiricalconvergence conditions in Eq. (13). n dD ( n ) g (a) 10 n D ( n )TV (b) 10 n dD ( n ) b (c)10 n cPD ( n ) (d) 10 n T ( n ) (e) 10 n S ( n ) (f) Figure 3: Convegence metrics (a) d (cid:101) D ( n ) g , (b) (cid:101) D ( n )TV , and (c) d (cid:101) D ( n ) b , and empirical convergence metrics (d) (cid:103) cPD ( n ) , (e) (cid:101) T ( n ) , and (f) (cid:101) S ( n ) , of the NCPD algorithm, as functions of n , from full-scan inconsistent data of the disk phantom. Even if the NCPD algorithm is demonstrated tosolve the non-convex optimization program in Eq.(8) above, evidence remains to be established as towhether the NCPD algorithm can invert the non-linear data model in Eq. (6). As discussed inSection 2.4.3 above, we seek the evidence only forconsistent data, i.e., g [ M ] = g [ NL ] ( b ) with b [truth] known, because the NCPD algorithm is developedbased upon the non-linear data model.Applying the NCPD algorithm to the full-scan simulated, consistent data described in Sec-tion 3.1.1 above, we reconstruct basis images b ( n ) ,calculate metric (cid:101) D ( n ) b in Eq. (16), and plot it (solid,6ark) in panel (h) of Fig. 2. It can be observed thatmetric (cid:101) D ( n ) b consistently decays, revealing that un-der the full-scan data condition, the NCPD algo-rithm can invert the non-linear data model in Eq.(6). We show in columns (b) and (d) of Fig. 1 recon-structed water and bone basis images of the numer-ical disk phantom, which are numerically identicalto their respective truth counterparts in columns(a) and (c) of Fig. 1. As the verification studyshows in Section 6, the NCPD algorithm can also,under short-scan data conditions, invert the non-linear data model in Eq. (6).The verification study not only yields evidencethat the NCPD algorithm can invert the non-lineardata model in Eq. (6) but also identifies data condi-tions under which the evidence can be established.Moreover, the result of a verification study likelyprovides a “best case scenario” in terms of recon-struction accuracy by use of the NCPD algorithm.
4. Results: monochromatic image recon-struction from simulated data
We now perform reconstruction studies onmonochromatic images from simulated data by us-ing the NCPD and CPD algorithms and demon-strate that the beam-hardening artifacts in theCPD image can be corrected effectively for in theNCPD image. For a benchmark purpose, we alsoreconstruct monochromatic images from the samedata by using an indirect method (Zou and Silver,2008), in which the FBP algorithm, along with aHanning kernel and a cut-off frequency of 0.3 cm − ,is used for reconstruction of the basis images fromthe estimated basis projections. In the simulation studies below, we use a circularfan-beam scanning configuration identical to thatdescribed in paragraph 2 of Section 3, except thatthe linear detector now consists of 1024 0.39-mmdetector bins.We create a numerical disk phantom to mimic thestandard DE-472 physical phantom (GAMMEX,2015) used in dual-energy CT imaging. The diskphantom contains water and bone (i.e., K = 2) ba-sis images of 512 ×
512 square pixels of size 0.49 mm.When each of the two spectra is used to generatedata from the phantom over a full rotation, we ob-tain a set of dual-energy data, which are referred toas full-scan data in the work.
Using Eq. (6) with the configuration and spec-tra parameters described, we generate from the nu-merical disk phantom full-scan simulated-noiselessdata g [ M ] at 640 views uniformly distributed over2 π . In the study, because g [ M ] = g [ NL ] ( b [truth] )and because truth image b [truth] is known, conver-gence conditions in Eqs. (11), (15), and (16) canbe calculated and are thus used. In this simulationstudy, truth-monochromatic image f [truth] m (cid:48) is knownbecause truth image b [truth] is known, and so is itsTV (i.e., truth-TV) value that is used as the valueof parameter γ .Applying the NCPD algorithm to the full-scansimulated-noiseless data, we reconstruct the basisimages and subsequently monochromatic image f m (cid:48) at energy 100 keV with Eq. (1), which are dis-played in column (a) of Fig. 4. We also show incolumns (b) & (c) of Fig. 4 images reconstructedfrom the same data by use of the CPD algorithmand the indirect method. It can be observed thatthe beam-hardening artifacts, including dark bandsconnecting the inserts and the cupping effect in thebackground material, in the CPD image can be cor-rected effectively for by use of the NCPD algorithm.The indirect method appears also to correct effec-tively for the beam-hardening artifacts. However,it yields an image of spatial resolution lower thandoes the NCPD algorithm. This can be observedfrom the edges of the inserts and enlarged pin-holestructures of negative contrast, as highlighted bythe arrows in Fig. 4. (a) NCPD (b) CPD (c) IndirectFigure 4: Monochromatic (top row) and their ROI (bottomrow) images at energy 100 keV reconstructed by use of (a)the NCPD and (b) the CPD algorithms, and (c) the indirectmethod, from simulated-noiseless data. The rectangle ROIis shown in column (b). The arrows highlight the pin-holestructures of negative contrast. Display window: [-150, 150]HU. .3. Monochromatic image reconstruction fromsimulated-noisy data We create full-scan simulated-noisy data byadding Poisson noise to the full-scan simulated-noiseless data in Section 4.2. The added noise levelcorresponds to 10 photons per detector bin in anair scan at each detector bin prior to the logarith-mic step. Because measured data g [ M ] now containnoise that is inconsistent with the non-linear datamodel in Eq. (6), i.e., g [ M ] (cid:54) = g [ NL ] ( b [truth] ), animage is reconstructed when the convergence con-ditions only in Eq. (11) are achieved.In the studies with noisy data here and real databelow, the truth-TV value may be neither appro-priate nor generally available, and we need to de-termine the value of parameter γ for a reconstruc-tion. For example, the value of γ could be de-termined that empirically yields a highest qual-ity/utility metric designed within a γ -value rangeof interest. Please see additional discussion of im-age reconstructions with different γ values in Sec-tion 5.1.2 below.In the study, we use the NCPD algorithm to re-construct, from the full-scan simulated-noisy data,a number of basis images with different values ofparameter γ , and then monochromatic images f m (cid:48) at energy 100 keV. Through visual inspection ofthe monochromatic images, we select subjectively γ = 423 . γ = 423 .
7. The images reconstructed by use ofthe CPD algorithm and indirect method are shownalso in columns (b) & (c) of Fig. 5. Again, thebeam-hardening artifacts in the CPD image can becorrected effectively for by use of the NCPD algo-rithm and the indirect method. However, the latteryields images with a noise level comparable to, butspatial resolution lower, than that of the former, ashighlighted by the arrows in Fig. 5.
5. Results: monochromatic image recon-struction from real data
In the study with real data below, we collect datafrom the GAMMEX DE-472 CT phantom by use ofa clinical CT scanner (Chen et al., 2018b), whichhas a curved-detector array of 896 1-mm bins form-ing a fan angle of 49.2 ◦ and a magnification fac-tor of 1.7. The GAMMEX DE-472 CT phantom (a) NCPD (b) CPD (c) IndirectFigure 5: Monochromatic (top row) and their ROI (bottomrow) images at energy 100 keV reconstructed by use of (a)the NCPD and (b) the CPD algorithms, and (c) the indirectmethod, from full-scan simulated-noisy data. The rectangleROI is depicted in Fig. 4. The arrows highlight the pin-holestructures of negative contrast. Display window: [-150, 150]HU. consists of calcium and iodine inserts (GAMMEX,2015) and is used widely in dual-energy CT studies.For each spectrum of 80 kVp and 135 kVp, data arecollected at 1200 views uniformly distributed over2 π . Again, we refer to the dual-energy data col-lected as full-scan real data. The images are recon-structed on an array of I = 512 ×
512 square pixelsof size 0.68 mm.
In the study with real data, the truth-TV valueis unknown. Therefore, as discussed above, we usethe NCPD algorithm to reconstruct, from the full-scan real data, a number of basis images with dif-ferent values of parameter γ , and then monochro-matic images f m (cid:48) at energy 100 keV. Through vi-sual inspection of the monochromatic images, weselect subjectively γ = 436 . γ =436 .
6, we subsequently obtain monochromatic im-age f m (cid:48) at energy 100 keV, which are displayed incolumn (a) of Fig. 6. For visual comparison, wealso show the images obtained with the CPD algo-rithm and the indirect method in columns (b) & (c)of Fig. 6, respectively. It can be observed that thebeam-hardening artifacts in the CPD image can be8orrected effectively for by use of the NCPD algo-rithm and the indirect method. However, the in-direct method yields an image that appears noisierand less sharp than that of the NCPD algorithm. (a) NCPD (b) CPD (c) IndirectFigure 6: Monochromatic (top row) and their ROI (bottomrow) images at energy 100 keV reconstructed by use of (a)the NCPD and (b) the CPD algorithms, and (c) the indi-rect method, from full-scan real data. The rectangle ROI isdepicted in Fig. 4. Display window: [-150, 150] HU. γ val-ues As γ plays a role in specifying the solution set ofthe non-convex program in Eq. (8), we carry out aseries of reconstructions by using the NCPD algo-rithm from the same real data with different valuesof γ , and display in Fig. 7 monochromatic images f m (cid:48) reconstructed at energy 100 keV (top row) fromthree selected γ values. It can be observed that theimages obtained with different values of γ are visu-ally different, as expected. In an application, as aspecified quality/utility metric is computed basedgenerally upon the image reconstructed, differentvalues of γ would result in different reconstructedimages and thus different values of the metric. Thisobservation also suggests that the parameter couldbe selected when yielding an empirically high valueof a quality/utility metric of interest (Zhang et al.,2016a; Xia et al., 2016; Zhang et al., 2016b). Furthermore, it may be of practical interest in in-specting how the reconstruction evolves as the iter-ation increases, especially in the form of monochro-matic images. We show in Fig. 8 monochromaticimages at energy 100 keV reconstructed with γ =436 . (a) γ = 378 . γ = 436 . γ = 1455 . γ =378.4,(b) 436.6, and (c) 1455.2, from full-scan real data by use ofthe NCPD algorithm. The rectangle ROI is depicted in Fig.4. Display window: [-150, 150] HU. along with the reconstruction obtained when theconvergence conditions in Eq. (11) are achieved. (a) iter. 30 (b) iter. 100 (c) iter. 500(d) iter. 1000 (e) iter. 2000 (f) ConvergentFigure 8: Monochromatic images at energy 100 keV obtainedat iteration (a) 30, (b) 100, (c) 500, (d) 1000, and (e) 2000,respectively, and (f) the convergent reconstruction from full-scan real data by use of the NCPD algorithm. Display win-dow: [-150, 150] HU.
6. Results: monochromatic image recon-struction from short-scan data
The NCPD algorithm may be exploited for inves-tigating accurate image reconstruction (or inversionof the non-linear data model in Eq. (6)) from a por-tion of the full-scan data corresponding to data col-lected with a non-standard scanning configurationof practical interest. In this work, we consider a9pecific non-standard scanning configuration, sim-ply referred to as the short-scan configuration indual-energy CT. In this configuration, two regularshort-scans are performed each of which is with adistinct spectrum. The two different spectra can beobtained, e.g., through a slow-kV switch technique,and the switching time defines the scanning angulargap between the end of the first regular short scanand the start of the second regular short scan, asillustrated in Fig. 9.
Figure 9: Illustration of (a) the full-scan configuration andshort-scan configurations with (b) 0 ◦ (b) and (c) 15 ◦ angulargaps, with the two circular curves depicting scanning angularranges of two distinct kVp spectra. First, from the full-scan simulated-noiseless datain Section 4.2, we extract two sets of short-scansimulated-noiseless data, corresponding to short-scan configurations with an angular gap of 0 ◦ or15 ◦ , mimicking an instantaneous or 15 ◦ -delayedkV-switch operation. In this simulation study,truth-monochromatic image f [truth] m (cid:48) is known be-cause truth image b [truth] is known, and so is itsTV (i.e., truth-TV) value that is used as the valueof parameter γ in the study.For each set of short-scan simulated-noiselessdata, using the NCPD algorithm, we reconstructthe basis images by using the NCPD algorithm,and then obtain the monochromatic image at 100keV. As shown in columns (b) & (c) of Fig. 10,the short-scan images are visually comparable tothat obtained with the full-scan data in column(a) of Fig. 10, showing the potential of the NCPDalgorithm for enabling the short-scan configura-tion. This verification study with noiseless consis-tent data on the new, non-standard short-scan con-figurations also surveys the data sampling condi-tions of the configurations and guides designing theconfigurations of promise before moving on to thecharacterization studies with inconsistent data. Itsets “best case scenarios” and performance upper-bar for the upcoming study with real data.Second, from the full-scan real data in Section 5,we extract two sets of short-scan real data, for the (a) Full (b) Short - 0 ◦ (c) Short - 15 ◦ Figure 10: Monochromatic (top row) and their ROI (bot-tom row) images at energy 100 keV reconstructed from (a)the full-scan simulated-noiseless data and from short-scansimulated-noiseless data acquired with an angular gap of (b)0 ◦ or (c) 15 ◦ . The rectangle ROI is depicted in Fig. 4.Display window: [-150, 150] HU. two short-scan configurations in Fig. 9. In this real-data study, we use the approach discussed above toselecting γ = 436 .
6. We display in columns (b) &(c) of Fig. 11 monochromatic images at 100 keVfor the two short-scan configurations. It can be ob-served that they are visually comparable to that ob-tained with the full-scan data in column (a) of Fig.11, demonstrating that the NCPD algorithm canenable the short-scan configuration in dual-energyCT. (a) Full (b) Short - 0 ◦ (c) Short - 15 ◦ Figure 11: Monochromatic (top row) and their ROI (bottomrow) images at energy 100 keV reconstructed from (a) thefull-scan real data and from short-scan real data acquiredwith an angular gap of (b) 0 ◦ or (c) 15 ◦ . The rectangle ROIis depicted in Fig. 4. Display window: [-150, 150] HU.
7. Discussion
We have developed the NCPD algorithm directlyto invert the non-linear data model in spectral CT.10onsidering the heuristic and iterative nature ofthe NCPD algorithm, we have designed and per-formed verification studies numerically to demon-strate that the NCPD algorithm can solve the non-convex optimization program and under appropri-ate data condition, can invert the non-linear datamodel in spectral CT. While the NCPD algorithm isadapted heuristically from the CPD algorithm thatis developed on a mathematically exact footing, itsheuristic component is limited only to the adapta-tion step. As such, the NCPD algorithm appearsto possess desirable features: it becomes the CPDalgorithm in the absence of the spectral effect, andthus the CPD algorithm provides a useful bench-mark of the NCPD convergence property; and itcan enable prototyping of different scanning config-urations, thus possibly yielding insights into designsof scanning configurations of practical workflow in-terest in spectral CT. Therefore, the NCPD algo-rithm may be exploited for helping design image-reconstruction algorithms of practical value.We have reconstructed images by using theNCPD algorithm from simulated and real data.The study results reveal that the NCPD algo-rithm can correct effectively for the non-linearbeam-hardening artifacts observed in images recon-structed with an algorithm for a linear data model.In the real-data study, we have also demonstratedthat the involved parameters such as γ can impactimage reconstruction, as expected, in the presenceof data components such as noise that are incon-sistent with the non-linear data model upon whichthe NCPD algorithm is based. We have also carriedout additional studies with different energy spectraand with numerical and physical phantoms of vari-ous levels of anatomy complexity of interest. Whileresults of these studies are not shown, observationsconsistent with those presented were obtained.The NCPD algorithm can accommodate data inwhich some of the rays are not measured with thefull set of multiple spectra, and can thus allow forenabling spectral CT with non-standard scanningconfigurations, as illustrated in the study exam-ples of short-scan dual-energy CT. The NCPD algo-rithm can also readily be applied to spectral dataacquired with photon-counting or sandwiched de-tectors. In the case of photon-counting detector,spectrum q sjm in the non-linear model in Eq. (2)would denote the response function within energychannel s for incoming photons at energy m andray j . Ideally, the response function would haveits support within the energy channel. Practically, due to physical factors such as charge sharing, theresponse function would have spill-overs to neigh-boring channels. As a result, direct algorithms suchas the NCPD algorithm that deal directly with thenon-linear data model in Eq. (2) may be suffi-ciently flexible to take into consideration such non-idealities.As shown in Eq. (8), the constraints on image TVand non-negativity are applied to monochromaticimage f m (cid:48) , instead of the basis images themselves,at an energy specified by energy bin m (cid:48) . Whilebeing assumed to be independent of spectral en-ergy (see Eq. (1)), the basis images reconstructedthrough solving Eq. (8) are possibly to be depen-dent on the choice of energy level m (cid:48) . It remains tobe investigated as to how constraints on monochro-matic images at different or combined energies inEq. (8) impact the basis images and consequentlymonochromatic images obtained.We have based the NCPD algorithm upon thespecific CPD algorithm. Additional algorithmsmay be developed for solving non-convex optimiza-tion programs by basing upon other algorithms forsolving the convex optimization programs. Again,whether such algorithms can accurately solve thenon-convex optimization programs and invert thenon-linear data model considered needs to be in-vestigated.The NCPD algorithm can be extended poten-tially to address other non-convex optimization pro-grams in spectral CT. It would be a future researchtopic of significance to evaluate the practical utilityof the NCPD algorithm in rigorous assessment stud-ies. The general approach taken, and the develop-ment of the CPD and NCPD algorithms, may alsoyield insights into solving non-convex optimizationprograms (Ochs et al., 2014; Valkonen, 2014), andinverting non-linear data models (Liu et al., 2019;Pan et al., 2015), in other imaging technologies andapplications.
8. Conclusion
We have developed an algorithm, referred toas the NCPD algorithm, to reconstruct imagesthrough solving a non-convex optimization programfor the non-linear data model in spectral CT. Re-sults of simulated- and real-data studies suggestthat the NCPD algorithm can yield accurate re-construction from full-scan data and that it hasthe potential to enable non-standard scanning con-figuration such as the short-scan configurations in11pectral CT. The work may also yield insights intosolving non-convex optimization programs and in-verting non-linear data models that arise in otherimaging modalities.
Appendix A. Derivation of the new CPDinstance of the general PD al-gorithm
We consider a linear data model g [ L ] ( b ) = H b + ∆ g c , (A.1)and a convex optimization program b (cid:63) = arg min b D g [ L ] ( b )s . t . || ( |∇ f (cid:48) m ( b ) | ) || ≤ γ & f (cid:48) m ( b ) (cid:23) , (A.2)where ∆ g c denotes a constant data vector indepen-dent of b , D g [ L ] ( b ) = || g [ M ] − g [ L ] ( b ) || / || ( g [ M ] − ∆ g c ) − H b || / , (A.3)and notations and symbols are consistent withthose in Section 2. It can be observed that thelinear term in Eq. (A.1) is identical to the linearcomponent of the non-linear model in Eq. (6) andthat the form of the convex optimization programin Eq. (A.2) is identical to that of the non-convexoptimization program in Eq. (8) except for the dif-ference between ∆ g c in the former and ∆ g [ NL ] ( b )in the latter. It is the similarity observed betweenthe two cases that motivates the development ofthe NCPD algorithm by adapting the CPD algo-rithm that is derived below for solving the convexoptimization program in Eq. (A.2).While the general PD algorithm can mathemati-cally exactly solve the convex optimization programin Eq. (A.2), a direct, brute-force implementationof some of the key operations, which are expressedonly in mathematically compact forms (e.g., theproximal operations,) in the general PD algorithmcan often lead to ambiguity issues concerning notonly computational accuracy and efficiency but alsoadditional parameters required, thus compromis-ing their original theoretical rigor and applicabil-ity. Also, the general PD algorithm provides lit-tle concrete knowledge of its convergence and ef-ficiency properties when solving a specific convexoptimization program. Therefore, we develop theCPD algorithm, as a new instance of the general PD algorithm, by deriving analytic, close-form so-lutions to the proximal operations specific to theconvex optimization program in Eq. (A.2) so thatthe solution can readily be implemented accuratelyand efficiently.The convex optimization program in Eq. (A.2)can be fitted into a primal minimization (Cham-bolle and Pock, 2010; Sidky et al., 2012) asmin x { F ( K x ) + G ( x ) } , (A.4)where, with g (cid:48) = g [ M ] − ∆ g c , F ( K x ) = F ( u , v , w ) = F ( u ) + F ( v ) + F ( w ) ,F ( u ) = 12 || u − g (cid:48) || , F ( w ) = δ P ( w ) ,F ( v ) = δ diamond( α m (cid:48) γ ) ( | v | ) , G ( x ) = 0 , (A.5) x = b , y = K m (cid:48) b = uvw , K m (cid:48) = H α m (cid:48) U m (cid:48) β m (cid:48) V m (cid:48) , u = H b , v = α m U m (cid:48) b , w = β m V m (cid:48) b , (A.6) x or b is the primal variable; y is the splitting vari-able that consists of three independent components, u , v , and w , obtained by applying matrix K m (cid:48) to b ; matrices U m (cid:48) = ( µ m ∇ , µ m ∇ , ..., µ mK ∇ )and V m (cid:48) = ( µ m I , µ m I , ..., µ mK I ); α m (cid:48) = ||H|| / ||U m || , β m = ||H|| / ||V m || . Here, α m (cid:48) and β m (cid:48) are multiplicative factors for the TV and non-negativity constraints in Eq. (18), which are com-puted in this work as the ratio of the matrix normsto scale and balance the magnitude of the matricesin stacked matrix K m (cid:48) (Sidky et al., 2012). Further, δ diamond( αγ ) ( x ) and δ P ( x ) are indicator functions inthe form of δ S ( x ) = (cid:40) x ∈ S ∞ x / ∈ S , (A.7)where S is a convex set. Here diamond( αγ ) and P are two sets defined as diamond( αγ ) = { x | || x || ≤ αγ } and P = { x | x i ≥ , ∀ i } . Similarly, δ ( x ) isalso an indicator function for a set of zero vectors,defined as { x | x i = 0 , ∀ i } .Using p , q , and r to denote the respective dualvariables of u , v , and w , we can obtain F ∗ and G ∗ ,the convex conjugate functions of F and G , as F ∗ ( p , q , r ) = F ∗ ( p ) + F ∗ ( q ) + F ∗ ( r ) ,G ∗ ( s ) = δ ( s ) , (A.8)12here F ∗ ( p ) = 12 || p || + g (cid:48)(cid:62) p ,F ∗ ( q ) = α m (cid:48) γ || ( | q | ) || ∞ ,F ∗ ( r ) = δ P ( − r ) . (A.9)The dual optimization of Eq. (A.4) takes the gen-eral form ofmax y {− F ∗ ( y ) − G ∗ ( −K (cid:62) y ) } , (A.10)which can be written explicitly here asmax p , q , r − || p || − g (cid:48)(cid:62) p − α m (cid:48) γ || ( | q | ) || ∞ − δ P ( − r ) − δ (cid:0) −H (cid:62) p − α m (cid:48) U (cid:62) m (cid:48) q − β m (cid:48) V (cid:62) m (cid:48) r (cid:1) . (A.11)The difference between the primal and dual objec-tives in Eqs. (A.2) and (A.11), without the indica-tor functions, is referred to as the conditional PD(cPD) gap (Sidky et al., 2012), which is given bycPD( b , p , q ) = 12 || g (cid:48) − H b || + 12 || p || + g (cid:48)(cid:62) p + α m (cid:48) γ || ( | q | ) || ∞ . (A.12)A proximal mapping of function F ∗ ( p ) is definedasprox σ [ F ∗ ]( p ) = arg min p (cid:48) (cid:26) F ∗ ( p (cid:48) ) + || p − p (cid:48) || σ (cid:27) . (A.13)For deriving the CPD algorithm of our interest, wefirst determine the proximal mappings of F ∗ and G .Noticing Eq. (A.8), the proximal mapping of F ∗ isgiven byprox σ [ F ∗ ]( p , q , r ) = prox σ [ F ∗ ]( p ) + prox σ [ F ∗ ]( q )+ prox σ [ F ∗ ]( r ) , (A.14)where prox σ [ F ∗ ]( p ) = p − σ g (cid:48) σ , (A.15)prox σ [ F ∗ ]( q ) = q − σ q | q | POLB α m (cid:48) γ ( | q | σ ) , (A.16)prox σ [ F ∗ ]( r ) = neg( r ) ≡ (cid:40) r i if r i ≤
00 if r i > . (A.17)The proximal mapping of G is given asprox τ [ G ]( x ) = x . (A.18) The CPD algorithm and its pseudo codes.
Substi-tuting the proximal mappings derived above intothe corresponding update steps in the general PDalgorithm (Sidky et al., 2012), we obtain the CPDalgorithm with the pseudo-code in Algorithm 1for solving the optimization problem in Eq. (A.2),where POLB γ ( x ) indicates a (cid:96) -projection onto the (cid:96) -ball of size γ (Duchi et al., 2008); neg( x ) denotesa non-positivity thresholding on each component of x . Algorithm 1
Pseudo-code of the CPD algorithm α m (cid:48) ← ||H|| / ||U m (cid:48) || ; β m (cid:48) ← ||H|| / ||V m (cid:48) || ; L m (cid:48) ← ||K m (cid:48) || ; σ m (cid:48) ← /L m (cid:48) ; τ m (cid:48) ← /L m (cid:48) ; θ ← Initialize b (0) , p (0) , q (0) and r (0) to zero ¯ b (0) ← b (0) repeat p ( n +1) = p ( n ) − σ m (cid:48) ( g [ M ] − ∆ g c − H ¯ b ( n ) )1 + σ m (cid:48) q (cid:48) ( n ) = q ( n ) + σ m (cid:48) α m (cid:48) (cid:80) k µ mk ∇ ¯ b ( n ) k q ( n +1) = q (cid:48) ( n ) − σ m (cid:48) q (cid:48) ( n ) | q (cid:48) ( n ) | POLB α m (cid:48) γ ( | q (cid:48) ( n ) | σ m (cid:48) ) r ( n +1) = neg( r ( n ) + σ m (cid:48) β m (cid:48) ( (cid:80) k µ mk ¯ b ( n ) k )) b ( n +1) = b ( n ) − τ m (cid:48) K (cid:62) m (cid:48) p ( n +1) q ( n +1) r ( n +1) ¯ b ( n +1) = b ( n +1) + θ ( b ( n +1) − b ( n ) ) until convergence conditions are satisfied General convergence conditions.
We devise generalconvergence conditions of the CPD algorithm as d (cid:101) D ( n ) g → (cid:101) D ( n )TV → d (cid:101) D ( n ) b → , (A.19) (cid:103) cPD ( n ) → (cid:101) T ( n ) → (cid:101) S ( n ) → , (A.20)as n → ∞ , where the dimensionless metrics aredefined as d (cid:101) D ( n ) g = | D g [ L ] ( b ( n ) ) − D g [ L ] ( b ( n − ) | / || g [ M ] || (cid:101) D ( n )TV = | || ( |∇ f (cid:48) m (cid:48) ( b ( n ) ) | ) || − γ | /γd (cid:101) D ( n ) b = || b ( n ) − b ( n − || / || b ( n − || , (A.21) (cid:103) cPD ( n ) = cPD ( n ) cPD (1) , (cid:101) T ( n ) = T ( n ) T (1) , and (cid:101) S ( n ) = S ( n ) S (1) , (A.22)conditional primal-dual (cPD) gap cPD ( n ) (Xiaet al., 2016), transversality T ( n ) (Hiriart-Urruty13nd Lemar´echal, 2013), and dual residualS ( n ) (Goldstein et al., 2013) are given bycPD ( n ) = cPD( b ( n ) , p ( n ) , q ( n ) ) , T ( n ) = ||H (cid:62) p ( n ) + α m (cid:48) U (cid:62) m (cid:48) q ( n ) + β m (cid:48) V (cid:62) m (cid:48) r ( n ) || , S ( n ) = || σ p ( n ) − p ( n − q ( n ) − q ( n − r ( n ) − r ( n − −K m (cid:48) ( b ( n ) − b ( n − ) || . (A.23)The conditions on d (cid:101) D ( n ) g and d (cid:101) D ( n ) b check thechange rates of the data-divergence and recon-structed image, whereas conditions on (cid:103) cPD ( n ) , (cid:101) T ( n ) , (cid:101) S ( n ) , and (cid:101) D ( n )TV form a set of first-order opti-mality conditions (Sidky et al., 2012; Hiriart-Urrutyand Lemar´echal, 2013) for solving Eq. (A.2). Data-convergence condition.
For measured dataconsistent with the linear data model, i.e., g [ M ] = g [ L ] ( b [truth] ), we can devise an additional conver-gence condition as (cid:101) D ( n ) g = D g [ L ] ( b ( n ) ) / || g [ M ] || → , (A.24)as n → ∞ . Clearly, for a case in which mea-sured data are inconsistent with the data model,i.e., g [ M ] (cid:54) = g [ L ] ( b [truth] ), this condition cannotbe used because data term D g [ L ] ( b ( n ) ) is gener-ally non-zero. It remains unknown, however, asto whether it can accurately invert the linear datamodel in Eq. (A.1). Image-convergence condition.
The convergenceconditions discussed above are used for demon-strating that the CPD algorithm can solve theconvex optimization program in Eq. (A.2). How-ever, evidence remains to be established that theCPD algorithm can invert the linear data model inEq. (A.1). We seek to establish the evidence underthe conditions that g [ M ] = g [ L ] ( b ) and that b [truth] is known, as the CPD algorithm is developed basedupon the linear data model. Therefore, we devisean image-convergence condition as (cid:101) D ( n ) b = || b ( n ) − b [truth] || / || b [truth] || → , (A.25)as n → ∞ . Appendix B. Numerical demonstration ofthe CPD algorithm
While the CPD algorithm is known mathemati-cally exactly to solve the convex optimization pro-gram in Eq. (A.2), it remains to be demonstrated to solve the program by achieving the general conver-gence conditions in Eqs. (A.19) and (A.20), andadditional convergence conditions in Eqs. (A.24)and (A.25) if measured data are consistent with thedata model. This would (1) ensure that the CPDalgorithm (containing the proximal operation in aclosed-form with a finite number of algebraic calcu-lation derived) is implemented correctly, (2) yieldconcrete knowledge about the convergence propertyof the CPD algorithm, and (3) provide an impor-tant benchmark of the convergence property for theheuristic NCPD algorithm adapted from the CPDalgorithm.
Appendix B.1. The CPD algorithm solving theconvex optimization program
Using the scanning configuration described inSection 4.1 above and the water and bone basis im-ages of a numerical disk phantom of size 128 × g [ M ] = g [ L ] ( b [truth] ) with the linear datamodel in Eq. (A.1) at 160 views uniformly dis-tributed over 2 π . We have assumed ∆ g c = .Also, truth-monochromatic image f [truth] m (cid:48) is knownbecause truth image b [truth] is known, and so is itsTV (i.e., truth-TV) value.Applying the CPD algorithm, along with γ thatis the truth-TV value, to consistent data g [ M ] , wecompute the convergence conditions in Eqs. (A.19),(A.20), and (A.24) and plot them as functionsof iteration n in panels (a)-(c), (d)-(f), & (g) ofFig. B.12, respectively. It can be observed that theconsistently decaying convergence curves (dashed,red) obtained demonstrate that the CPD algorithmcan solve the non-convex optimization program inEq. (A.2) in terms of achieving the convergence con-ditions.We have also performed a study with inconsistentdata g [ M ] that are obtained by addition of Poissonnoise to consistent data g [ L ] ( b [truth] ) used in theverification study above. As a result, Eq. (A.24)can no-longer be used as a convergence condition inthe study in which measured data (i.e., noisy datahere) are inconsistent with the linear data model,i.e., g [ M ] (cid:54) = g [ L ] ( b [truth] ). The study results, whichare not shown here, again demonstrate that theCPD algorithm can solve the non-convex optimiza-tion program in Eq. (A.2) in terms of achieving thegeneral convergence conditions in Eqs. (A.19) and(A.20).14 n dD ( n ) g (a) 10 n D ( n )TV (b) 10 n dD ( n ) b (c)10 n cPD ( n ) (d) 10 n T ( n ) (e) 10 n S ( n ) (f)10 n D ( n ) g (g) 10 n D ( n ) b (h) Figure B.12: Convergence metrics (a) d (cid:101) D ( n ) g , (b) (cid:101) D ( n )TV , (c) d (cid:101) D ( n ) b (d) (cid:103) cPD ( n ) , (e) (cid:101) T ( n ) , (f) (cid:101) S ( n ) , (g) (cid:101) D ( n ) g , and (h) (cid:101) D ( n ) b of the CPD algorithm, as functions of iteration n , from con-sistent data, obtained with single (dashed, red) and double(solid, dark) computer precision. Appendix B.2. The CPD algorithm inverting thelinear data model
While the CPD algorithm is demonstrated aboveto solve the convex optimization program in Eq.(A.2) above, it remains to be shown, however, asto whether it can accurately invert the linear datamodel in Eq. (A.1), and as to what data conditionsunder which it can. The verification study belowis thus to demonstrate that under appropriate datacondition, the CPD algorithm can invert the lineardata model in Eq. (A.1) in terms of achieving Eq.(A.25).In panel (h) of Fig. B.12, we display (cid:101) D ( n ) b (dash,red) computed from reconstructed images b ( n ) as afunction of iteration n . It can be observed that met-ric (cid:101) D ( n ) b decreases consistently until the four-bytefloating-point precision is reached, thus demon-strating numerically that under the full-scan datacondition, the CPD algorithm can invert the lineardata model in Eq. (A.1). Appendix B.3. Numerical demonstration ofachieving convergence conditions
It should be noted that in Fig. B.12, the con-vergence curves (dashed, red) decay until around ∼ iterations where we believe that four-bytefloating-point precision (i.e., single floating-pointprecision) is reached. To confirm this observation, we then perform the same study, but with eight-byte floating-point precision (i.e., double floating-point precision.) It can be observed that the conver-gence curves (solid, dark) continue their decayingtrends passing ∼ iterations until reaching thedouble floating-point precision around ∼ itera-tions. The convergence curves are plotted on log-logscales for revealing unambiguously not only theirvalues, but also equally importantly their decayingtrends. The result verifies that the convergence con-ditions are achieved up to the single floating-pointprecision around ∼ iterations. We have alsoconducted the studies in Section 3.1.1 for the NCPDalgorithm with the double floating-point precision,and results similar to those observed for the CPDalgorithm have been obtained. Acknowledgements
This work is supported in part by NIH R01 GrantNos. EB026282 and EB023968, and the Grayson-Jockey Club Research.
References
Alvarez, R.E., Macovski, A., 1976. Energy-selective recon-structions in X-ray computerised tomography. Phys. Med.Biol. 21, 733–744.Barber, R.F., Sidky, E.Y., 2016. Mocca: Mirrored con-vex/concave optimization for nonconvex composite func-tions. The Journal of Machine Learning Research 17,5006–5056.Barber, R.F., Sidky, E.Y., Schmidt, T.G., Pan, X., 2016. Analgorithm for constrained one-step inversion of spectralCT data. Phys. Med. Biol. 61, 3784–3818.Cai, C., Rodet, T., Legoupil, S., Mohammad-Djafari, A.,2013. A full-spectral Bayesian reconstruction approachbased on the material decomposition model applied indual-energy computed tomography. Med. Phys. 40,111916.Chambolle, A., Pock, T., 2010. A first-order primal-dual al-gorithm for convex problems with applications to imaging.J. Math. Imaging Vis. 40, 120–145.Chen, B., Zhang, Z., Sidky, E.Y., Xia, D., Pan, X., 2017.Image reconstruction and scan configurations enabled byoptimization-based algorithms in multispectral CT. Phys.Med. Biol. 62, 8763–8793.Chen, B., Zhang, Z., Xia, D., Sidky, E., Pan, X., 2018a. Non-convex Chambolle-Pock algorithm for multispectral CT,in: Proceedings of The Fifth International Conference onImage Formation in X-Ray Computed Tomography, pp.377–381.Chen, B., Zhang, Z., Xia, D., Sidky, E.Y., Pan, X., 2018b.Algorithm-enabled partial-angular-scan configurations fordual-energy CT. Med. Phys. 45, 1857–1870.Duchi, J., Shalev-Shwartz, S., Singer, Y., Chandra, T., 2008.Efficient projections onto the l1-ball for learning in highdimensions, in: Proceedings of the 25th international con-ference on machine learning, ACM. pp. 272–279. AMMEX, 2015. Dual energy CT phantom. URL: .Goldstein, T., Li, M., Yuan, X., Esser, E., Baraniuk, R.,2013. Adaptive primal-dual hybrid gradient methods forsaddle-point problems. arXiv preprint arXiv:1305.0546 .Hernandez, A.M., Boone, J.M., 2014. Tungsten anode spec-tral model using interpolating cubic splines: Unfiltered x-ray spectra from 20 kv to 640 kv. Med. Phys. 41, 042101.Hiriart-Urruty, J.B., Lemar´echal, C., 2013. Convex analy-sis and minimization algorithms I: Fundamentals. volume305. Springer science & business media.Hubbell, J., Seltzer, S., 2004. Tables of X-ray mass attenu-ation coefficients and mass energy-absorption coefficients(version 1.4). URL: http://physics.nist.gov/xaamdi .Liu, X., Chen, B., Zhang, Z., Xia, D., Sidky, E.Y., Pan, X.,2019. Optimization-based reconstruction for correctingnon-linear partial volume artifacts in CT, in: Proc. SPIEMed. Imag.: Phys. Med. Imag., p. 109482Q.Long, Y., Fessler, J., 2014. Multi-material decomposition us-ing statistical image reconstruction for spectral CT. IEEETrans. Med. Imag. 33, 1614–1626.Mory, C., Sixou, B., Si-Mohamed, S., Boussel, L., Rit, S.,2018. Comparison of five one-step reconstruction algo-rithms for spectral ct. Phys. Med. Biol. 63, 235001.Nakada, K., Taguchi, K., Fung, G.S.K., Amaya, K., 2015.Joint estimation of tissue types and linear attenuation co-efficients for photon counting CT. Med. Phys. 42, 5329–5341.Ochs, P., Chen, Y., Brox, T., Pock, T., 2014. iPiano: inertialproximal algorithm for nonconvex optimization. SIAM J.Imaging Sci. 7, 1388–1419.Pan, X., Chen, B., Zhang, Z., Rose, S., Sidky, E.Y., 2015.Optimization-based simultaneous determination of emis-sion activity and photon attenuation in PET, in: Proc.13th Int. Meet. Fully 3-D Image Reconstruct. Radiol.Nucl. Med., pp. 83–86.Rigie, D.S., LaRiviere, P.J., 2015. Joint reconstruction ofmulti-channel, spectral CT data via constrained total nu-clear variation minimization. Phys. Med. Biol. 60, 1741–1762.Roessl, E., Proksa, R., 2007. K-edge imaging in X-ray com-puted tomography using multi-bin photon counting detec-tors. Phys. Med. Biol. 52, 4679–4696.Sawatzky, A., Xu, Q., Schirra, C., Anastasio, M., 2014. Prox-imal ADMM for multi-channel image reconstruction inspectral X-ray CT. IEEE Trans. Med. Imag. 33, 1657–1668.Sidky, E.Y., Jorgensen, J.H., Pan, X., 2012. Convex op-timization problem prototyping for image reconstructionin computed tomography with the Chambolle-Pock algo-rithm. Phys. Med. Biol. 57, 3065–3091.Sidky, E.Y., Pan, X., 2008. Image reconstruction in circularcone-beam computed tomography by constrained, total-variation minimization. Phys. Med. Biol. 53, 4777.Valkonen, T., 2014. A primal–dual hybrid gradient methodfor nonlinear operators with applications to MRI. InverseProblems 30, 055012.Xia, D., Langan, D.A., Solomon, S.B., Zhang, Z., Chen,B., Lai, H., Sidky, E.Y., Pan, X., 2016. Optimization-based image reconstruction with artifact reduction in C-arm CBCT. Phys. Med. Biol. 61, 7300–7333.Zhang, R., Thibault, J.B., Bouman, C., Sauer, K., Hsieh,J., 2014. Model-based iterative reconstruction for dual-energy X-ray CT using a joint quadratic likelihood model. IEEE Trans. Med. Imag. 33, 117–134.Zhang, Z., Han, X., Pearson, E., Pelizzari, C., Sidky, E.Y.,Pan, X., 2016a. Artifact reduction in short-scan CBCTby use of optimization-based reconstruction. Phys. Med.Biol. 61, 3387–3406.Zhang, Z., Ye, J., Chen, B., Perkins, A.E., Rose, S., Sidky,E.Y., Kao, C.M., Xia, D., Tung, C.H., Pan, X., 2016b. In-vestigation of optimization-based reconstruction with animage-total-variation constraint in PET. Phys. Med. Biol.61, 6055–6084.Zhao, Y., Zhao, X., Zhang, P., 2015. An extended algebraicreconstruction technique (E-ART) for dual spectral CT.IEEE Trans. Med. Imag. 34, 761–768.Zou, Y., Silver, M.D., 2008. Analysis of fast kv-switching indual energy CT using a pre-reconstruction decompositiontechnique, in: Proc. SPIE Med. Imag.: Phys. Med. Imag.,p. 691313..Liu, X., Chen, B., Zhang, Z., Xia, D., Sidky, E.Y., Pan, X.,2019. Optimization-based reconstruction for correctingnon-linear partial volume artifacts in CT, in: Proc. SPIEMed. Imag.: Phys. Med. Imag., p. 109482Q.Long, Y., Fessler, J., 2014. Multi-material decomposition us-ing statistical image reconstruction for spectral CT. IEEETrans. Med. Imag. 33, 1614–1626.Mory, C., Sixou, B., Si-Mohamed, S., Boussel, L., Rit, S.,2018. Comparison of five one-step reconstruction algo-rithms for spectral ct. Phys. Med. Biol. 63, 235001.Nakada, K., Taguchi, K., Fung, G.S.K., Amaya, K., 2015.Joint estimation of tissue types and linear attenuation co-efficients for photon counting CT. Med. Phys. 42, 5329–5341.Ochs, P., Chen, Y., Brox, T., Pock, T., 2014. iPiano: inertialproximal algorithm for nonconvex optimization. SIAM J.Imaging Sci. 7, 1388–1419.Pan, X., Chen, B., Zhang, Z., Rose, S., Sidky, E.Y., 2015.Optimization-based simultaneous determination of emis-sion activity and photon attenuation in PET, in: Proc.13th Int. Meet. Fully 3-D Image Reconstruct. Radiol.Nucl. Med., pp. 83–86.Rigie, D.S., LaRiviere, P.J., 2015. Joint reconstruction ofmulti-channel, spectral CT data via constrained total nu-clear variation minimization. Phys. Med. Biol. 60, 1741–1762.Roessl, E., Proksa, R., 2007. K-edge imaging in X-ray com-puted tomography using multi-bin photon counting detec-tors. Phys. Med. Biol. 52, 4679–4696.Sawatzky, A., Xu, Q., Schirra, C., Anastasio, M., 2014. Prox-imal ADMM for multi-channel image reconstruction inspectral X-ray CT. IEEE Trans. Med. Imag. 33, 1657–1668.Sidky, E.Y., Jorgensen, J.H., Pan, X., 2012. Convex op-timization problem prototyping for image reconstructionin computed tomography with the Chambolle-Pock algo-rithm. Phys. Med. Biol. 57, 3065–3091.Sidky, E.Y., Pan, X., 2008. Image reconstruction in circularcone-beam computed tomography by constrained, total-variation minimization. Phys. Med. Biol. 53, 4777.Valkonen, T., 2014. A primal–dual hybrid gradient methodfor nonlinear operators with applications to MRI. InverseProblems 30, 055012.Xia, D., Langan, D.A., Solomon, S.B., Zhang, Z., Chen,B., Lai, H., Sidky, E.Y., Pan, X., 2016. Optimization-based image reconstruction with artifact reduction in C-arm CBCT. Phys. Med. Biol. 61, 7300–7333.Zhang, R., Thibault, J.B., Bouman, C., Sauer, K., Hsieh,J., 2014. Model-based iterative reconstruction for dual-energy X-ray CT using a joint quadratic likelihood model. IEEE Trans. Med. Imag. 33, 117–134.Zhang, Z., Han, X., Pearson, E., Pelizzari, C., Sidky, E.Y.,Pan, X., 2016a. Artifact reduction in short-scan CBCTby use of optimization-based reconstruction. Phys. Med.Biol. 61, 3387–3406.Zhang, Z., Ye, J., Chen, B., Perkins, A.E., Rose, S., Sidky,E.Y., Kao, C.M., Xia, D., Tung, C.H., Pan, X., 2016b. In-vestigation of optimization-based reconstruction with animage-total-variation constraint in PET. Phys. Med. Biol.61, 6055–6084.Zhao, Y., Zhao, X., Zhang, P., 2015. An extended algebraicreconstruction technique (E-ART) for dual spectral CT.IEEE Trans. Med. Imag. 34, 761–768.Zou, Y., Silver, M.D., 2008. Analysis of fast kv-switching indual energy CT using a pre-reconstruction decompositiontechnique, in: Proc. SPIE Med. Imag.: Phys. Med. Imag.,p. 691313.