DLIMD: Dictionary Learning based Image-domain Material Decomposition for spectral CT

 Abstract — The potential huge a dvantage of spectr al computed tomography (CT ) is it s capability to provide accuracy m aterial identifica tion and qu antitative tissue information. This can benefit clinica l application s, such as brain angiography , early tumor recognition, etc . To achieve m ore accura te m aterial components with higher material image quality , we develop a dictionary learning based image-domain material deco mposition (DLIM D) for spectral CT in this paper. First , we reconstruct spect ral CT image fro m projections and ca lculate m aterial coefficients m atrix by selecting uniform regions of basis materials fro m image reconstructio n results. Second, we employ the direct inversion (DI) m ethod to obtain initial material decomposition resu lts , and a set o f image patches are extracted from the mode-1 unfolding of nor malized material image tensor to train a united dictionary by the K-SVD technique . Third, the trained di ctionary is employed t o explore the similarities from decomposed material images by constructing the DLIM D model. Fourth, more constraints ( i.e., volume conserva tion a nd the bounds of each pixel w ithin material maps) are furt her integrated into the model t o improve the accuracy of material decomposition. Finally, both physical phanto m and preclinical experi ments are employed to evaluate the performance of the proposed DLIMD in materia l decomposition a ccuracy, material im age edg e preserv ation and feature reco very. Index Terms — spectra l computed t omography ( CT ), material decomposition , image doma in, dictionary learning. I. I NTRODUCTI ON HE conventional co mputed tom ography (CT ) tr eats x-ray spectrum as single-energ y, and an energy integration detector (EID) is em ployed to collect photons fro m the This w ork w as suppor ted i n part by the N ational Natural Scie nce Foundation of China (No. 61471070 and No . 61501310) and NIH /NIBIB U01 grant (EB017 140). W. Wu, H.J. Yu, P. Chen, F. Liu and J. Feng are with Key Lab of Optoele ctronic T echnology and Systems, Ministry of Education, Chongqing University, Chongqing 4000 44, China. F. Luo is with State Key Laboratory of Information Engineering in Survey ing, Mapping and Remote Sensing, Wuhan University, Wuhan 4 30079, China. Q. Wang and H. Yu are with the Department of Electrical and Computer Engineering, University of Massachusetts Lowell , L owell, MA 01854, USA. Y. Zhu is with Schoo l of Mathematical Sciences, Capital Nor mal University, Beijing, 100048, China. He is also with Beijing Higher I nstitution Engineering Research Center of Testing and I maging, Beijing, 100048, China. Y. Zhang is w ith the PingAn Technolo gy, US Rese arch Lab, Palo Al to, CA 94306, USA . F. L iu and H. Y u serve as the correspondence authors (e -mail: liufl@cqu.edu.c n and hengy ong-yu@ieee.or g ). spectrum . As a resul t, d ifferent material co mponents may ha ve the same o r si milar linear attenuation coefficients i n the reconstructed CT images. T h is implies that it fails to discriminate and d ecompose material co mponents within t he imaging o bjects [1, 2], and it further loses the fu nctional information for CT images. Fortunately, due to t he advantages of dose redu ction, tissue contrast improvement, quantitative tissue analysis, beam hard ening artifacts reduction, material discrimination a nd deco mposition [3], the spectral CT has b een developed as an alter native solution. Compared with the traditional CT sc hemes, spectr al CT s ystems ca n obtain two or more ener gy bin proj ections for the sca nned obj ects. From t his view point , the dual energ y CT (DECT ) , as a simple commercialized architec ture, has achieved a great success in practical applications [4, 5] . Ho wever, the projectio ns from DECT o nly contain two energ y-spectru m measure ments. This limits t he m aterial decomposition capability of DECT in two basis materials without ad ditional constraint s [6, 7] . The development of photon- co unting detectors (P CDs) makes spec tral CT a hot topic in recent years [8, 9] . A typical spectral CT system equ ipped with P CD can co llect multi ple projections for d ifferent energy bins within one scan. Ho wever, different materials may have si milar attenuation coefficients in some of the energy b ins. Besides, how to perform a quantitative analysis for tis sue distribution is still a n op en problem. T o realize the aims of material discrimination and tissue q uantitative analysis, it is necessar y to develop advanced material deco mposition methods for spectr al CT. Regarding the material decomposition meth ods for spectral CT , they can be divided into two categories: direct and indirect methods [10] . Specifically, a direc t material decomposition method can directly obtain material components from pro jections [9, 11- 13] with known x-ray spectr um . However, it is difficult to ach ieve x-ra y transmission spectral m odel in practice. B esides, although some regularization prior ( for examples, total variation (TV) [9], nonlocal TV [11] and so on) have been considered in such a material decomposition model, the results fro m direct material deco mposition are still sensitive to noise. Indirect m aterial d ecomposit ion methods ca n be further divided into proj ection-based and image-based methods [14] . For the pr ojection -based methods, the proj ections are first decomposed to the speci fied material si nograms. Then , an image reco nstruction algorithm (for example, FBP [15] ) is followed. For such material decomposition methods, the err ors within the specified basis material si nograms ca n be magnif ied in the material maps, a nd this further comprise material decomposition accuracy. For the image-ba sed m ethods, the 1 st step is to r econstruct spectral CT i mages from proj ections, and the 2 nd step is to reconstruct material maps . In terms o f the 1 st step, there are a lot of iterative image reconstr uction methods DLIMD : Dictionary Learning based Image - domain Material Decomposition for spectral CT Weiwen Wu, Haijun Yu, Peijun Chen, Ful in Luo, Fenglin Liu, Member, IEEE , Qian Wang, Yining Zhu, Member, IEEE , Y anbo Zhang, Senior Member, IEEE, Jian Feng, Heng yong Y u, Senior Member, IEEE T available, such a s channel-in dependent total variatio n (TV) [16] , HighlY constrained backPRojectio n (HYPR) [17] , tight frame sparsity [18], p atch-based lo w-rank [3] , spectral p rior image co nstraint compressed sensing (spectral PICCS) [19] , prior rank and sparsit y [2 0], tensor dictionar y learning ( TDL) [21] and its advanced version (L0T DL) [22], nonlocal sparse matrix dec omposition [23], spatial-spectral cube match ing frame (SSCM F) [24], and non-local lo w-rank c ube-based tensor factorization (NLCTF) [25]. Regarding the 2 nd step , numer ous methods were p roposed for DE CT [26- 29] . However, the methods for multiple-bin spectral CT are r elative scarce . For examples, Tao et.al [30] develo ped a prior knowledge aware iterati ve denoising material d ecomposition (MD -PKAD) model by co nsidering the Pr ior Image Constrain ed Compressed Sensing (PICCS) [31]. Xie et.al proposed a multiple constraint model for image domain mater ial decomposition and validated it only by numerical simulation [32] . Besides, the deep learning bas ed multi-material decomposition models [33 , 34] ar e also studied. In this work , we will focus on the 2 nd step for image d omain material decomposition. The dictionary lear ning-base d m ethods have o btained a great achieve ment i n the CT field , such as lo w-dose reconstruction [35], sparse-view reconstruction [36] , m edical image denoising [37], spectral CT reconstruction [38] , etc . Indeed, the d ictionary lear ning-based m ethods can recov er image details and feature s b y training a lar ge number of image patches. Si milar to t he hyperspectral image[39], considering the correlatio n among di fferent ener gy channels (i.e., image structure similarities of the same o bject ), the tensor based dictionary lear ning was developed for spectral CT reconstruction a nd then obtain excellent per formance i n our previous studies [21, 2 2]. As for the multiple material decomposition o f spectral CT, multiple material i mages can be considered as a tensor, leading to a natural idea to estab lish tensor based dictionary learning for image-domain materi al decomposition. However, since it is difficult to ensure a rank- 1 tensor within a local region for multiple materials, it is inappropriate to ad opt the tensor dictio nary learning to formulate material decomposition model. Another idea is to develop a multiple-dictionary based dec omposition metho d, i.e., each material correspo nds to o ne dictio nary. However, the correlation of material images will be ignored. Besides, because some materials only have simple structures (see iodine contrast agent in Fig. 1(d ) in section III ), the trained dictionary cannot co de material information i n the decomposition process. To overcome the aforementioned limitations, we formulate a unified dictionary by training image patches from mod e-1 unfolding o f nor malized material image tensors. T hen, we construct a dictio nary learning based image-domain material decomposition (DLIMD) m o del to fully code the sp arsity and similarities of material i mages simultaneously. Fi nally, the volume co nservation a nd bound of each pixel are introduced into DLIMD m odel to fu rther improve material decomposition accurac y. T he advantages of DLIMD method include good cap abilities of noise suppression and edge preservation tha n other competing met hods. The rest of this paper is organized as follo ws. In section II, the material deco mposition model for spectral CT will be given. I n sectio n III, w e elab orate the m athematical model and solution for the DLIMD. In section IV, two real dataset s are employed to evaluate the proposed decomposition method. In section V, we d iscuss some related issues and make a conclusion. II. M AT ERIAL D ECOMPOSITI ON M ODEL FOR S PECT RAL CT A. Spectral CT imagin g model Assuming the detected photon number for the   x-ray path within the  th energy windo w   is   (      ). It can be modelled as       󰇛  󰇜       󰇛   󰇜         󰇛  󰇜 where     integrates o ver the range of  th energy channel and    indicates integral along the   x-ray p ath. In this study , M basis materials within the o bject is assumed,   󰇛    󰇜 represents the  th (      ) m aterial linear attenuation coefficient for energ y E at po sition  , and    󰇛    󰇜   correspond s to the summation of m aterial attenuation at position  .   󰇛  󰇜 represents the original x -ray photon intensity emit ting from the x -ray source for ener gy  . Considering t he basis material expansion , the summation o f x-ray attenuation coef ficient    󰇛    󰇜   in Eq. (1) can be furtherly written as a low-dimensio nal expansion. T hat is,    󰇛    󰇜       󰇛  󰇜   󰇛  󰇜    󰇛  󰇜 where   󰇛  󰇜 represents the m ass-attenuation coefficient for the m th material at energy  which ca n refer to the tab les in [40], and   󰇛  󰇜 is co mponent for  th material at location  . The task of m aterial decomposition is to recover all the material co mponent maps   󰇛  󰇜󰇛      󰇜 from Eq.(1). I f we direct ly reco nstruct m aterial m aps from Eq. (1), it can be considered as direc t material decomposition [10] . In this stud y, we only consider the i ndirect material deco mposition metho ds, i.e., image domain material decomposition, where the first step is to reconstruct the summation o f equivalent mass- attenuation co efficient f or    󰇛        󰇜   . Let   󰇛󰇜 be the original x-ray photon flux which ca n be expressed as   󰇛󰇜     󰇛  󰇜     󰇛  󰇜 Considering the equivalent mass-atten uation co efficient    󰇛        󰇜   and substituting Eqs. (2 ) and (3) into Eq. ( 1), we have      󰇛  󰇜       󰇛        󰇜      󰇛  󰇜 Considering t he multi-ener gy CT i maging model, t he averaged attenuation coefficients      󰇛        󰇜 is employed to replace    󰇛        󰇜   . Eq. (4) can be discretize d as following     󰇛  󰇜           󰇛       󰇜  󰇛  󰇜 where      󰇛        󰇜 is the projection matrix,   and   represent w idth and height o f the avera ged attenuation coefficients image,  is the number of total x-ray paths, and   presents the   row of  .      󰇛      󰇜    is the vectorization of averaged attenuation coefficient image at n th energy windo w . Then , a logarithm operation is operated on both sides of Eq . ( 5) . By ignor ing the index  , Eq. (5 ) can b e simplified as       󰇧     󰇛  󰇜 󰇨        󰇛  󰇜 For the image-do main materia l d ecomposition method, we can reconstruct multi -energy CT images      󰇛       󰇜 from the corresponding   󰇛       󰇜 . T here are numerous advanced multi -energy CT image reconstruction techniques, including TDL[21], L 0 TDL[22], NLCTF [41] and so on. In this stud y, we focus the image -do main material deco mposition, and FBP and NLCTF are e mployed to obtain 󰇝      󰇞   . B. Image- domain Material dec omposition model In or der to obtain the material component maps   󰇛  󰇜󰇛      󰇜 from 󰇝      󰇞   , Eq.(2) can be further r ead as a matrix form                                       󰇛  󰇜 where   represents a veraged at tenuation coe fficients of the  th material at  th energy window, which can b e deter min ed by calc ulating t he averaged mass coefficients[ 42]. Eq. (7) is equivalent to  󰇛  󰇜   󰇛  󰇜  󰇛  󰇜 where                     ,         and         are two tensors w hich resp ectively rep resent the reconstructed images and m aterial images.  󰇛  󰇜     and  󰇛  󰇜     are the mode-3 un folding of  and  . However, the reconstructed spectral CT i mages usually contain noise, which can comprise the d ecomposed m aterial images q uality. Considering noise in the reconstructed ima ges, Eq. (8) can be modified as  󰇛  󰇜   󰇛  󰇜   󰇛  󰇜 To recover spectral CT images from Eq. (9 ), w e have the following least square linear pro gram problem:       󰇛  󰇜   󰇛  󰇜     󰇛  󰇜 where     is Frob enius nor m. E q. (10) ca n be viewed as image domain basi s material decomposition model. The simple direct inversion (DI) method [26] can b e employed to solve Eq.(10) and obtain:  󰇛  󰇜  󰇛    󰇜     󰇛  󰇜  󰇛  󰇜 T o improve the decomposed i mage quality , using regularization is a feasib le strateg y to constrain t he solution during iteration pro cess, and E q. (10) can be re formulated as   󰇧    󰇛  󰇜   󰇛  󰇜      󰇛  󰇜 󰇨  󰇛  󰇜 where  is a regularization facto r to b alance the data fidelity term    󰇛  󰇜   󰇛  󰇜    and regularization ter m  󰇛  󰇜 . III. D I CTIONARY L EARNI NG BASED I MAGE - DO MAIN M ATERIAL D ECOMPOSITI ON (DLIMD) A. Dictionary learnin g A set of image patche s      ,       , are extracted from the traini ng datasets and unfold ed as vectors to train the dictionary     , where      and  represents the number of atoms. T he goal of dictionar y learn ing is to find sparse representation coef ficients      using the dictionary  . T hi s can b e establ ished b y sol ving the following optimization prob lem:                              󰇛  󰇜 where   is the sparsity level ,     is the quasi-   norm,       is sparse r epresentation coefficient s for i th image patch. Eq. (13) is a co nstrain ed li near programming pro blem, w hich can be converted into the follo wing unconstrain ed pr oblem:       󰇛                 󰇜     󰇛  󰇜 where   is a Lagrange multiplier. Eq. (14) can be solved by utilizing an alter nating minimization scheme. T he first step is to update   by fixing the dictionar y  ,       󰇛                 󰇜     󰇛  󰇜 Eq.(15) ca n be optimized using t he matching p ursuit (MP) [43] or orthogonal matching purs uit (OM P) algorith m [ 44].The second s tep is to update the dictionary with a fixed spars e representation coefficients  . There are many methods available to train the d ictionary, suc h as K -SVD [45] and online learning tech nique[46] . B. DLIMD Method Considering the general i mage do main material deco mposition model Eq . (12), s imilar to the dictio nary lear ning for lo w-dose CT reconstruction [35], th e proposed Dictionary Learning based Image-do main Material Decomposition (DLIMD) for spectral CT can be constructed as follow:    󰇝   󰇞         󰇛  󰇜   󰇛  󰇜          󰇧   󰇛   󰇜                󰇨          󰇛  󰇜 where   is m th channel of material images  ,       is sparse representation coefficients for i th image patch from m th material i mage and    󰇝   󰇞    .   󰇛   󰇜 is the i th i mage patch extractio n oper ator from   .   is the trained dictionar y. In this study, we use a united dictionary for different mater ial decomposition. High -quality dictionar y is b eneficial for sparse representation and accuracy of material decomposition . A natural idea is to train di fferent   for different material s . However, th is strateg y will be limited by three reaso ns. First , a specific material map may o nly contain a few image feat ures, i.e., bone and iodine contrast in Fig.1. I t is difficult for the trained   to encode the i mage features and r educe sp arse representation level during the p rocess of material decomposition, and this co mp rises t he material decomposition accuracy. Seco nd, tr aining different   is time consumi ng , which im plies higher co mputational cost is needed with various materials in practic e. T hird, if we train the dict ionary   for each material, the correlation between di fferent material maps will be lost. To overcome these limitations , we formulate a unified dictionar y       in this study. For a given spectral C T dataset, we can obtain the FBP or oth er reconstruction res ults and apply DI method to achieve material decomposition results  . T hen, before training the dictionary   , we should normalize m aterial images to avoid d ata inconsistency o f the pi xel values from differe nt materials . After that, a set of image patches are extracted from  󰇛  󰇜 󰆒      󰇛    󰇜 , the mode- 1 unfolding of normalized image  󰆒 , to train a global dictionary   . Alth ough the DI m aterial maps contain noise, our following experi ments de monstrate that the dictionary trainin g can preserve image in formation and against noise well . Finally, t he dictionary   can b e lear ned usin g t he K-SVD algorithm [45]. T he advantages of the for mulated unified dictionary   are threefold. First, it can fully encode similarities w ithin different materials images (for examples, the image structures within blue boxes from bone and soft tissue, th e i mage structures within red boxes from soft tissue and iodine contrast in Fig. 1 ) during the material decomposition proce ss. Second, it can enhance the redundanc y with the trained dictionary   . T hi s is good to encod e current material i mage structures b y employing the ato ms from ot her material image to enhance sparse representation. Third, we can also save training ti me to some extent in practice . Fig.1 Numerical simulation mo use phantom used in [21] consists of three materials, i.e. bone , soft tissue a nd iodine. (a) , (b), (c) and (d) repre sent the phantom, bone, soft tissue and iod ine, respectively . To further i mprove the ac curac y o f material d ecomposition, more constraint s sho uld b e considered. For example, it is reasonable to assume that bone and contrast agent don’t co - exist in a given pixel , and this assu mption is emplo yed in our experiments. I f the air is also treated as one b asis material, the summation o f pixel values at the same location in different material images should be equal to o ne [7 ], i.e.,                    󰇛                󰇜  󰇛  󰇜 where       represents the pixel value at 󰇛        󰇜  location.  is the air map and      represents the pixel value at 󰇛      󰇜  location. Besides, the pixel value within  should be in the ran ge of [0 1] , i.e.,           󰇛  󰇜 To further constrain the solution for optimization and improve the accurac y of material d ecomposition, Eqs. (1 7)-(18) can be treated as t wo conditions and are i ncorporated into the material deco mposition model Eq. (16) . Then, we have:    󰇝   󰇞         󰇛  󰇜   󰇛  󰇜          󰇡  󰇛   󰇜                 󰇢                                           󰇛  󰇜 To solve Eq. (19), we first introduce a tensor  to replace  , and E q. (19) can be convert ed in to the following constrain ed programmable prob lem:    󰇝   󰇞          󰇛  󰇜   󰇛  󰇜          󰇡  󰇛   󰇜                 󰇢                                             󰇛  󰇜 Eq. ( 20 ) is equal to the following unco nstrained o ptimization problem with certain conditio ns:    󰇝   󰇞          󰇛  󰇜   󰇛  󰇜                   󰇡  󰇛   󰇜                 󰇢                                          󰇛  󰇜 Here    is an optimization facto r to balance the estimation  for  . T he objection function E q. ( 21) can be divided into the following t wo s ub-problems:       󰇛  󰇜   󰇛  󰇜          󰇛  󰇜                                     󰇛   󰇜   󰇝   󰇞           󰇛  󰇜          󰇧   󰇛   󰇜                󰇨          󰇛  󰇜 i)  sub-prob lem: Regarding the op timization of Eq .(22a) , there ar e two strategies . A natural first strateg y is to treat the air as a basis material to  . In this case, the material attenuation matrix  can be modified as  󰆒                        󰇛   󰇜 , where  is a small positive value to represent air attenuation coefficie nt. The air can increase the number of b asis materials in this strategy, a nd it can result in instability of material deco mposition and comprise the deco mposition acc uracy, especiall y in the ca se o f    . The second strategy for so lving Eq. (22 a) can be divided into two step s . T he first step is to solve the objec tion function Eq. (23 ):       󰇛  󰇜   󰇛  󰇜          󰇛  󰇜                               󰇛  󰇜 Considering Eq . ( 23) with pixel level, it can be furtherl y written as       󰇼 󰇛  󰇜    󰇛  󰇜  󰇼                        󰇛  󰇜                                              󰇛  󰇜 where  󰇛  󰇜   󰇣                󰇤  ,  󰇛  󰇜   󰇣         󰇤  ,                         and       󰇛  󰇜        󰇛  󰇜          󰇛  󰇜   . Not e that        , Eq. ( 24 ) is equivalent to                                       󰇛  󰇜                                              󰇛  󰇜 We optimize Eq. ( 25 ) by minimizing the follo wing proble m:                                          󰇛  󰇜                                    󰇛  󰇜 Eq. ( 26 ) is a constrain ed convex pr ogrammable op timization problem . I t is equal to the follo wing optimization problem           󰇛      󰇜                       󰇛  󰇜                                   󰇛  󰇜 Eq. ( 27 ) is a constrain ed leas t square problem, which can be easily solved. The second step is to f ind the air region by operating a th reshold m ethod on the  󰇛   󰇜 . A gain, if the pixel value withi n  󰇛  󰇜 is greater t han a given t hreshold, the pixel v alue can be replaced by direct inversion. The given threshold is set as 0.99 in this work. To compare strategies 1 and 2, Fig .2 shows the material decomposition r esults using d irect inversion (DI) method from FBP images with noise -free projections . From Fig. 2, it can be seen t hat so ft tissue r esults using s trategy 1 contain severe strip artifacts. However, the artifacts are suppressed by u sing strategy 2. ii)  su b-problem: As for the optimization pr oblem of Eq.(22b) , it can be re-written as    󰇝   󰇞              󰇛   󰇜          󰇡  󰇛   󰇜                 󰇢         .( 28 ) Eq. ( 28 ) can b e divided into M sub-pro blems. Again,   󰇛       󰇜 can be updated independentl y. Here, Eq . (28) can be furtherly read as Fig.2 Material decompos ition results w ith soft tissue of the mouse phantom. (a) and (b) repre sent the resul ts using strategies 1 and 2, respectively .                 󰇛   󰇜            󰇛   󰇜      󰇛   󰇜                         󰇛  󰇜  Eq. ( 29 ) is a t ypical dictionar y learning b ased i mage denoisi ng model, and it can be solved by the method in [47] . During the process of d ictionary learni ng b ased image denoising, the parameters of spar sity level L and tolerance o f representation error  pla y i mportant roles in co ntrolling the d ictionary quality and material dec omposition accurac y. In this stud y, the sparsity le vel L is set as the sa me value for different materia ls. However, the tolerance of representation error     󰇛       󰇜 should be careful chosen. T he overall pseu do - co d es of the DLIMD algorithm is listed in Algorithm I. Al gorithm I : DLIMD Input :         and o ther parameter s; I nitialization of  󰇛󰇜    󰇛󰇜      . Output: Material decompositio n tensor  . Part I: Diction ary trainin g 1: Reconstruct ing spectral CT images us ing certain image reconstructio n method; 2: For mulating the material attenuation matrix using reconstruction results; 3: Decomposing t he reconstructe d images usi ng DI me thod; 4: Normalizing the D I results; 5: Extracting image patches to f orm a diction ary training da taset; 6: Training a dict ionary using the K-SVD tech nique. Part II : Material de composition 7: While not co nverge nce do 8: Upda ting  󰇛  󰇜 using Eq. (27); 9: Finding and pro cessing the p ixels in  󰇛  󰇜 using DI technique; 10: Updating  and 󰇝   󰇞   using Eq. (29) ; 11: k=k+1; 12: En d While C. Algorithm imp lementation d etails From the pseudo-codes of DLIMD method, it can be divided into two main parts : dictionar y learning and m aterial decomposition . To implement the dictio nary learning, we first reconstruct spectral CT i mages fro m pr ojections using analytic/iterative reconstruction methods . In thi s st udy, both FBP and NL CTF[2 5] methods are employed to i mplement image reconstruction. Cal culating the basis material attenuation matrix b y selectin g un iform regions is a key step for material d ecomposition. Then, w e adop t the DI method to decompose recon structed i mages into basi s material maps using the for mulated material attenuation matrix. Finally, 10 4 image patches with size o f    are extracted from normalized material images t o tr ain dictio nary    using t he K- SVD algorit hm . T he d ictionary   is o vercomplete to enforce the sparsity level . Agai n, the number of ato ms in a d ictionary should be much greater tha n t he size of an ato m, i.e .,      . Usually ,  should b e greater than      . Here, the number o f ato ms  within the dictionar y   is set as 512 . The sparsity level   in the dictionary training can be set empirically from 5 to 10, and it is uniformly set as 6. The iteration number to train dictionar y is uniformly se t as 20 0. Fig. 3 shows three trained dic tionaries and will be used in the following experi ments. Fig. 3 . The d ictionarie s used i n our experiments . 1 st -2 nd colum ns represent the dictionaries fo r physical phan tom and preclinic al experime nts. Especially, Eq. (27) is solved by adopting the lsq lin function in Matlab 201 4(b). The parameter selection s o f sparsity representation level L and     󰇛        󰇜 are dependent on specified application. All the m ethods are stopped after 30 iteratio ns . IV. E XPERIMENTS A ND R ESULTS T o evaluate the perfor mance of our proposed DLIMD method in material deco mposition, two real datasets ar e employed. T he advantages are demonstrated by co mparing with the DI and total variation-regularization material decomposition (TV MD ) methods. To quantitatively evalu ate the performance of all material decomposition tec hniques, the root m eans square error (RMSE), peak-signal- to - noise ratio (PSNR) and structural si milarity (SSIM) are employed. T he optimized parameters used in all exper iments ar e li sted in Table I. Table I . Parameters use d in the proposed DLI MD me thod for experiments Parameters  L      Physical phantom 0.003 10 (0.08, 0.03, 0.00 4) Preclinical (FBP) 0.001 12 ( 0. 004 ,0.012,0. 05 ) Preclinical (NL CTF) 0.001 12 (0.001,0.004,0.0 04) A . Physical Phanto m A p hysical phanto m consists o f three basis materials, i.e., water, iod ine a nd aluminu m, is employed to eval uate o ur method first. As shown in Fig. 4 , the phantom contains five cylinders, and ea ch o f them represents di fferent b asis material or dif ferent conce ntrations of iodine solution . A luminum is used to mimic bo ne. The spectral CT system i s eq uipped with a m icro-focus x-ray so urce (YXLON, 225Kv) and a flat-panel PCD (Xcounter, XC -Hydra FX20) . For the central slice, the PCD contains 2048 d etector cell s and each of them covers a length of 0. 1 mm. Here, e very 4 cells are combined to reduce noise . Because the P CD only has t wo energy c hannels, th e phantom i s scanned twice to obtain 4 energy bin projections . The full-scan proj ections are collected fro m 1080 unif ormly distributed vie ws with 137kV x-ray so urce . T he distances from the x-ray source to rotation al center and detector are 182.68 mm and 440.50 mm , respectively. The FOV dia meter is calculated as 82 .6 mm . Since the r econstructed material image matrix is 256 × 256, each pixel covers an area of 0.324 × 0.324 mm 2 . Fig. 4. Setups of p hysical phantom experiments . (a) i s the spectral CT system , (b) and (c) repre sent the physical phantom. To evaluate the propo sed material dec omposition methods in the ca se of analytic reconst ruction results, F ig. 5 shows t he reconstructed spectral CT images from four different energ y bins using FBP . Because the x-ray has beam hardening , the aluminum cylinder region is sl ightly damaged by met al artifacts. This can compris e the m aterial decomposition accuracy to so me extent. In this study, to validate t he advantages o f constraint E q. (1 7) in m aterial deco mposition model, the DI without equivalent transfor mation constraint (DIWET) , i.e. , Eq. (17), is also co mpared. The m aterial decomposition results w ith t hree basis material s (aluminum, water and iodine) are give n in Fi g. 6 . From Fig. 6 it can be observed that the constraint Eq. (17) can improve m aterial images q uality remarkably compar ed with the DIWE T and DI results. Besides, the material deco mposition with regularization terms can also provide higher image quality than that obtai ned by the DI method. In terms o f al uminum results, the p roposed DLIMD method can better protect image edges and avoid the blocky artifacts compared with the TVMD method. As for the d ecomposed water results, the proposed DLIMD can pro vide much clearer image ed ge with higher accuracy, whic h is confirmed b y the ma gnified ROI. Regarding the iodine results, three cylinders from the DLIMD are more co mplete a nd uni form tha n those o btained by other decomposition methods. Fig. 5 Phy sical phantom FBP reconstruction re sults. Fr om lef t to right, the images are fo r 1 st -4 th energy bins with a display window [0 1.3]cm -1 . Fig. 6 Material decompositio n results from Fig.5 . Fro m left to right, the columns represent the decompositio n results of aluminum, water and iodine , where t he display windows are [0.5 1], [0 .8 1] and [0 0.015 ]. From top to bottom, the rows represent the results decomposed by DIWET, D I, TVMD and DLI MD methods, respe ctively . To further evaluate the p erfor mance of DLIMD algorithm in im proving material decomposition , five region of interests (ROI) indicated with 1 -5 are extracted from different material maps a nd t heir quantitative evaluation results are li sted in Table II . Fro m Tab le II , we can see the performance of DIWET is alw ays the w orst . The proposed m ethod always performs the best follo wed by the TV -based and DI methods. Table II . Q uantitative eval uation results of ROI s 1-5. ROI-1 DIWET DI TVMD DLI MD RMSE 0.1890 0.0850 0.0827 0.0804 PSNR 14.471 21.406 21.649 21.893 SSIM 0.9278 0.9542 0.9823 0.9925 ROI-2 RMSE 0.2631 0.0324 0.0318 0.0273 PSNR 11.599 29.793 29.964 31.2 76 SSIM 0.7705 0.9732 0.9762 0.9978 ROI-3 RMSE( 10 -4 ) 48 .00 26 .09 9.754 6.719 PSNR 46.416 51.668 60.216 63.453 SSIM 0.3343 0.4071 0.7889 0.8972 ROI-4 RMSE( 10 -4 ) 66 . 50 34 .35 14 .43 10.281 PSNR 43.544 49.281 56.816 59.761 SSIM 0.4192 0.6105 0.8710 0.9474 ROI-5 RMSE( 10 -4 ) 87.83 62 .22 38 .06 33.640 PSNR 41.127 44.122 48.389 49.463 SSIM 0.4443 0.6114 0.8605 0.92 75 In this work, al l t he i nvolved methods are p rogrammed with Matlab (version 2014 b) on a PC (i7-6700, 8.0 GB mem ory). The p roposed DLIMD algorithm mainly incl udes t wo subroutines: dictionar y training and material decompositi on . As for the co mputatio nal cost with material decompositio n , the higher the number of the materials is, t he greater the computational cost is. Particularl y, t he DIWET , DI , TVMD and DLIMD consume 0.04 , 22 .79 24.16 and 44.17 seconds, respectively . B esides, the DLIMD needs more time (387.22 seconds) to train t he diction ar y co mpared with other regularization based methods. Obviously, the D LIMD algorithm r equires more computational co st than ot her algorithms. B . Preclinical exp eriments In this st udy, the PILATUS3 PCD w ith 4 ener gy-channels by DECTRIS is employed to acquire multi-energ y proj ections in preclinical application. Such PCD consists o f 515 cells a nd each has a length o f 0 .15 mm. Here, a full-sca n i s p erformed to obtain projections with 4 energy bi ns fro m 720 vie ws. Fig. 7 (a) sho ws the imaged objec t consisti ng of chicken foot and 5 mg/mL iodine sol ution cylinder . The distances fro m x-ra y source to rotational center a nd detector are 35.27 cm a nd 43.58 cm, resp ectively . The FBP reconstruction result s (see Fig. 7 (b)-(e) ) have    pixels and ea ch o f the m co vers an area of        . Fig. 7 Artificial mix ed tissues experime nt. (a) is the preclinical specime n fixed on the spectral CT sy stem . (b)-(e ) are FBP reconstruct ion results from 4 energy bins, where the dis play window is [0 0.5] cm -1 . To evaluate the advantages o f the DLIMD method , the material deco mposition res ults fro m FBP results with differ ent algorithms are shown in Fi g. 8 . For bon y components, the image edge of bon y structure is blurred in DI results. Although the TV based material dec omposition met hod can p rovide clear image edge, the block y artifacts ar e ap peared in component map. This point can be confirmed by the bony structure indicate d w ith red arrow within th e m agnified bony ROI “A” in Fig. 9 . As for so ft tissue deco mposition res ults, the DLIMD can obtain clear image str ucture and features. To clarify this point, three soft tis sue ROIs marked with “B”, “C” and “D” are extracted and magnified in Fig. 9 . From Fig. 9 , it can be easily seen that image ed ge a nd struct ures are remarkably clear tha n those o btained by T VMD and DI methods. Similar to the soft tissue components , as for iodine contrast, the DI can wrongly classif y the pixel s of bone component into iodine contrast f ollowed b y T VMD results . Obviously, t he DLIMD method can furtherl y i mprove the accuracy of decomposed materials. This conclusion can b e observed from the magnified ROI “E” in 3 rd row of Fig. 8 . According to Fig. 7( a) , the iodine concentration is 5.0 mg/mL and it ca n b e considered as gold standard to evaluate the decomposition acc uracy of different algorith ms. Here, the RMSEs a nd mean va lue o f extracted ROI 6 are calculated in Table III . From T able III, Here, the DLIMD can obtain th e largest mea n val ue of iodine and the s mallest RMSE, followed by TVMD and DI methods with FBP rec onstruction case. Fig. 8 Material decompositio n results from Fig .7 . From left to right, the columns re present the DI , TV MD and DL I MD methods. T he 1 st -3 rd ro ws represe nt the bone, soft tissue and iodine with the display windows [0. 25 0 .5], [0.85 0.95] and [0.0018 0. 005], respective ly. The advanced spectral CT imag e reconstruction methods can i mprove the accurac y of material decomposition by providing hi gh q uality recon structed i mages [25] . To further demonstrate the advanta ges of the pr oposed DLIMD algorit hm in advanced reconstruction result case , the NLCT F reconstruction method is emplo yed to implement ima ge reconstruction and corresp onding final material deco mposition results usi ng different algor ithms are furtherly given in Fig. 10 . Since the TVMD ca n obtain higher material decompositi on accuracy t han DI method, only the deco mposition results from TVMD and NLCTF m ethods are given to save space i n Fig. 10 . Although the TVMD method can provide the higher quality of bone component map, the block y ar tifacts st ill appear. T his point can be confirmed by the ex tracted ROIs “F” and “G” which are further mag nified in in Fig. 11 . From 1 st - 2 nd columns of Fig . 11 , it can be observed that the image structure fro m TVMD results indicated b y “3” conta ins blocky artifacts. Fo rtunately, thi s structure can be observed in DLIMD r esults. Besides, the gap bet ween t wo b one s indica ted by ar row “4” is still b lurred . Ho wever, it is clear in the proposed DLIMD method. As for so ft tissue decomposition results, the D LIMD ca n still obtain higher ima ge quality with clear i mage edges a nd much feature. Spec ifically, the image structure is blurred in TVMD in Fig . 10 , but DLIMD can provide clear im age ed ge . Especially, two ROI labelled by “H” and “K” are extracted and magn ifi ed in Fig s. 10 and 11 . The image features indicated by ar rows “2” and “6” ar e disappeared in TVMD. T hey can s till be seen i n DLIMD results. In add ition, the image str ucture is too blurr ing to observe shape of edge. However, w e can see t he details of edge in DLIMD results. Here, the MSEs and m ean value from location of ROI-6 are also calculated in Table III. From T able III, one can see the DLIMD can obtain the smallest MSE. In terms of mean val ue, the TVMD and DLIMD can o btain the same accuracy. Fig. 9 The magnified R OIs from Fig.8 . From left to right, th e columns represe nt the DI, TVMD and DLI MD methods. The 1 st -4 th rows represent the ROIs marked with “A ”, “B”, “C” and “D”, where t he di splay windows are [0. 29 0. 33 ], [0.85 0.95 ], [0.85 0.95] and [0.85 0.9 5], respective ly. Table I II . Quanti tative evaluation re sults of ROI -6. DI TVMD DLI MD FBP Mean (%) 0.1864 0.2282 0.2283 MSE ( 10 -6 ) 11.62 7.750 7.629 NLCTF Mean (%) 0.18 30 0.2353 0.2353 MSE ( 10 -6 ) 10.36 7.245 7.234 Fig. 10 Material decomposition results from t he NLCTF . From left to right , the columns represent the b one, so ft tissue and iodine with the display window s [ 0.25 0.5], [0.90 0.92] and [0.0018 0.005], respectively . The 1 st -2 nd rows re present TVMD and DLI MD methods. Fig. 11 The magnified ROI s in Fig.10 . 1 st -2 nd rows represe nt the TVMD and DLI MD me thods. The 1 st -3 rd columns represent the ROI s marked with “F”, “G” and “H”, where th e display window s are [0.3 0.45], [0.29 0.33] and [0.89 0.94], respectivel y. V. D ISCUSSI ONS AND C ONCLUSI ONS To improve the accurac y of material decomposition for spectral CT, we p ropose a DLIMD technique. Co mpared with the previo us i mage-do main material deco mposition meth od for spectral CT [12], the innovations of DLIMD are threefold . First, considering the similarities o f different material images, we co nstruct a unified dictionar y to code material image sparsity by training a set of image p atches. Here, t hose image patches are extracted from nor malized material ima ges. Second, w e formulate a DLIMD mathematical model by enhancing sparsity of materia l maps with the dictiona ry and analyzing the image-domain material decomposition. Third, additional constraints ( volume conservation [7] and the bound of each material pixel val ue) are incorporated into basis material deco mposition model to further improve the decomposition accuracy. T wo real datasets from physic al phantom and preclinical d ataset are employed to evaluate the proposed m ethod, and the results de monstrate the advantages of DLIMD i n improving image quality and material decomposition accurac y. As a m aterial deco mposition method in image d omain, the DLIMD ca n remove the artifacts to so me e xtent or impr ove the accuracy of material images from x-ray beam hardeni ng. For exa mple, although aluminum co mponent in the employ ed physical phantom ca n cause F BP image qualit y reduction a nd contain x-ray beam hard ening artifacts. Ho wever, the DI, TVMD and DLIMD methods with e nhanced constraints can still provide high accurac y of material deco mposition results. Because the reconstr ucted spectral CT images usually contain noise and artifacts, th ey can cause the instabil ities of image domain material deco mposition. A feasible strateg y is to for mulate regularizatio n- based decomposition mode ls. While TV is a co mmon regul arizer to be chosen, we i ntroduce an advanced sparsity representation (i.e., dictionary learning) into the material deco mposition model for better perfor mances in ter ms of image feature recovery and edge pr otection . Although t he D LIMD method can obtain satisfied decomposition results in image-do main, there ar e still some remaining i ssues that should be addressed. First, the parameters are only e mpirically optimized by comparing extensive experiment results. For example, r egarding t he selection of  , it is necessary to select an   for each material component indep endently. However, if  is great, it is difficult to select ap propriate  . In this case, it is feasible t o select a reaso nable  by est imating noise le vel within all material i mages in our follow -up work. Besides, we will also try some auto matic strategies to op timize other parameters in the DLIMD method. Seco nd, t wo real datasets o nly contain three different bas is materials. However, the i maging objects may co ntain multiple (greater than 3) materials. It is nece ssary to validate and eva luate the p erformance o f DLIMD in such cases. In conclusion, b ased on the d ictionary learning theory, we propose a DLIMD method for image-do main material decomposition of spectral CT. Real dataset exper iments demonstrate the outperfor mance of D LIMD technique. T his will be extremely useful for image d omain mate rial decomposition in practical ap plications of spectral CT. R EFERENCES [1] P. M . Shikhaliev and S. G. Fritz, "Photon counting spectral CT versus conventional CT: comparative evaluation for breast imaging application," Physi cs in Medici ne & Biology, vol. 56, p. 1905, 2011. [2] X. Z. Lin, Z. Y. Wu, R. Tao, Y. Guo, J. Y. Li, J. Zhang , et al. , "Dual energy spectral CT imaging of insulinoma — value in preoperative diagnosis compared w ith conventional mult i-detector CT," European journal of radiol ogy, vol. 81, pp. 2487-2494, 2012. [3] K. Kim, J. C. Ye, W . Worstell, J. Ouyang, Y. Ra kvongthai, G . El Fakhri , et al. , "Sparse-vie w spectral CT reconstr uction using spec tral patch- based low-rank penalty," IEEE t ransactio ns on medical imaging, vol. 34, pp. 748-760, 2015. [4] H. Saito, K. Noda, K. Ogasawara, S. Atsuji, H. Takaoka, H. Kajihara , et al. , "Reduced iodinated contrast media for abdominal imaging by dual - layer spectral dete ctor computed tomo graphy for patients with kidney disease," R adiology Case Reports, vol . 13, pp. 437-443, 2018/04/01/ 2018. [5] Y. Liu, A. Liu, L . Liu, S. Tian, J. Liu, R . Pu , et al. , "F easibil ity of spectral i maging with low ‐ concentration contrast medium in abdominal CT angiography of obese patients," International journal of clinical practi ce, vol. 70, 2016. [6] T. R. Johnson, B. K rauss, M. Se dlmair, M. Grasruck, H. Bruder, D. Morhard , et al. , "Material differentiation by dual energy CT: initial experie nce," European Radiolo gy, vol. 17, pp. 1 510-1517, 2007. [7] X. L iu, L . Yu, A . N. Primak, and C. H. McCollo ugh, "Qu antitative imaging of element composition and mass fraction u sing dual ‐ energy CT: Three ‐ material decomposition," Medical physics, vol . 36, pp. 1602-1609, 2009. [8] K. Taguchi and J. S. Iwanczy k, "Vision 20/20: Single photon counting x ‐ ray detectors in medical imaging," Medical Physics, vol. 40, p. 100901, 2013. [9 ] R. F. Ba rber, E. Y. Sidky, T. G. Schmidt, and X. Pan, "An algorithm for constrained one-step inversion of spectral CT d ata," Physics in Medicine & Biology, vol . 61, pp. 3784-3818 , 2016. [10] W. Wu, Q. Wang, F. Liu, Y. Zhu, and H. Yu, "Block Matching Fram e based Material Recons truction for Spe ctral CT," arXiv preprint arXiv:1810.1034 6, 2018. [11] J. Liu, H. Ding, S. Molloi, X. Zhang, and H. Gao, "TICMR: Total image constrained material reconstruction via n onlocal total variation regulariz ation for spectral CT," IEEE transactions on medical imaging, vol. 35, pp. 2578- 2586, 2016. [12] Y. Zhao, X. Zhao, and P. Zhang, "A n extended algebraic reconstruction technique (E-ART) for d ual spectral CT," IEEE Trans Med Imaging, vol. 34, pp. 761-8, Mar 2015. [13] Y. Long an d J. A . Fessler, "Multi-material decomposition using statistical image reconstruction for spectral CT," IEEE tra nsactions on medical imagin g, vol. 33, pp. 1614-1626, 201 4. [14] Q. Wang, Y. Zhu, a nd H. Yu, " Lo cally linear constraint based optimization model for material decomposition," Physics in Medic ine & Biology, vol. 62, p . 8314, 2017. [15] W. Wu, H. Yu, S. Wang, and F. Liu, "BPF-type region- of -in terest reconstructio n for parallel translational computed tomography ," Journal of X-ray science and technology, v ol. 25, pp. 487-50 4, 2017. [16] Q. Xu, H. Yu, J. Bennett, P. He, R. Zai non, R. Doesburg , et al. , "I mage Reconstruction for Hy brid True-Color Mic ro-CT," IEEE Transactions on Biomedical En gineering, vol. 59, pp. 1711-17 19, 2012. [17] S. Le ng, L. Yu, J. W ang, J. G. Fletcher, C. A. Mistretta, and C. H. McCollo ugh, "Noise reduction in spectral CT: R educing dose and breaking the trade ‐ off b etwee n image noise a nd energy bin sel ection," Medical physics, vol . 38, pp. 49 46-4957, 201 1. [18] B. Zhao, H. G ao, H. D ing, and S. Molloi, "Tight ‐ frame based iterative image reconstruction for spectral breast CT," Medical p hysics, vol. 40, 2013. [19] Z. Yu, S. Leng, Z. Li, and C. H. McCollo ugh, "Spectral prior image constrained compre ssed sensing (spectral PICCS) for phot on -counting computed tomography ," Physics in Medicine & Biology, vol. 61, p. 6707, 2016. [20] H. Gao, H. Yu, S. Osher, and G. Wang, "Multi -energy CT ba sed on a prior rank, intens ity and sparsity model (PRISM)," Inverse problems, vol. 27, p. 115012, 2011. [ 21] Y. Zhang, X. Mou, G . Wang, and H. Yu, "T ensor-based dictionary learning for spectral CT reconstruction," IEEE t ransactions on medical imaging, vol. 36, p p. 142-154, 201 7. [22] W. Wu, Y. Zhang, Q. Wang, F. Liu, P. Chen, and H. Yu, "L ow-dose spectral CT reconstructio n using image gradient ℓ0– norm and te nsor dictionary," Applied Mathematical Modelling, vol. 63, pp. 538-557, 2018/11/01/ 201 8. [23] S. Ni u, G . Yu , J. M a, and J. W ang, "No nlocal low-rank and sparse matrix decompositio n for spectral CT re construction," Inverse Problems, vol. 34, p. 02400 3, 2018. [24] W. Wu, Y. Z hang, Q. W ang, F . Liu, F. L uo, and H. Yu, "Spatial- Spectral Cube Matching Frame for Spectral CT Reconstructio n," inverse problems, vol . 34, p. 104003, 2018. [25] W. Wu, F. Liu, Y. Zhang, Q. Wang, and H. Yu, "Non-local Low-r ank Cube-based Tensor Factorization for Spectral CT Reconstruction," IEEE transactions on m edical im aging, vol. 38, pp. 1 079-1093, 2019. [26] T. Niu , X. Dong, M. Petrongolo , and L. Zhu, "I terative image ‐ domain decomposition fo r dual ‐ energy CT," Medical p hysics, vol. 41, 2014. [27] C. Maaß , M. Baer, and M. Kachelrie ß , "I mage ‐ based dual energy CT using optimized precorr ection functions: A practical n ew approach of material decomposition in image d omain," Medical physics, vol. 36, pp. 3818-3829, 2009. [28] W. Zhao, T. Niu, L. Xing, Y. Xie, G. Xiong, K. Elmore , et al. , "Using edge-pre serving algorithm with non-local mean for significantly improved image -domain mate rial decomposi tion in dual-energy CT," Physics in Me dicine & Biolog y, vol. 61, p. 1332, 2016. [29] Y. Xue, R . Ruan, X. Hu, Y. Kuang, J. Wang, Y. Long , et al. , "Statistical image ‐ domain multimaterial decompositio n for dual ‐ energy CT," Medical physics, vol . 44, pp. 88 6-901, 2017. [30] S. Tao, K. Rajendran, C. H. McCollough, and S. Le ng, "Material decomposition with prior k nowle dge aware iterative denois ing (MD- PKAID) ," Physics in Medicine & Biology, vol. 63, p. 195003, 2018. [31] G. H. Chen, J. Tang, and S . Leng, "Prio r Image Constrained Compre ssed Sensing (PI CCS)," Medical Phys ics, vol. 35, p. 66 0, 2008. [32] B. Xie, T. Su, V. Kaftandjian, P. Niu, F. Yang, M. Robini , et al. , "Material Decomposition in X-ray Spectral CT Using Multiple Constraints in I mage Domain," Journal of Nondes tructive Evaluatio n, vol. 38, p. 16, 20 19. [33] D. P. Clark, M. Holbrook, and C. T. Badea, "Multi-energy CT decomposition usi ng convolutional neural netwo rks," in Medical Imaging 2018: P hysics of Me dical Im aging , 2018, p. 105731O . [34] Z. Ch en and L . Li, "Robust multimaterial decomposition of spectral CT using convolutional neural networks," Optical Engineering, vol. 58, p. 013104, 2019. [35] Q. Xu, H. Yu, X. Mou, L. Zhang, J. Hsieh, and G. Wang, "Low-dose X- ray CT reconstructio n via dictionary lear ning," IEEE Tr ansactions o n Medical Imagin g, vol. 31, pp. 16 82-1697, 2012. [36] S. Li, Q. Cao, Y. Chen, Y. Hu, L. Luo, and C. Toumoulin, "Dictionary learning based sinog ram inpainting for CT sparse recons truction," Optik- International Journal for Light and Electron Optics, vol . 125, pp. 2862- 2867, 2014. [37] S. Li, H. Yin, and L. Fang, "Group-sparse representation with dictionary learning for medical image d enoisi ng and fusion," IEEE Transactions on biomedical engi neering, vol. 59, p p. 3450-3459, 2012. [38] Z. Bo, D. Huanjun, L. Yang, W. Ge, Z . Jun, and M . Sabee, "Dual - dictionary learning-based iterative image reconstruction for spectral computed tomography ap plication," Physics in Medicine & Biology, vol. 57, p. 8217, 20 12. [39] F. Luo, B. Du, L . Zhang, L . Zhang, and D. T ao, "Fe ature Learning Us ing Spatial-Spectral Hypergraph Discriminant Analysis for Hyperspe ctral Image," IEEE Trans actions on C ybernetics, 20 18. [40] J. H. H ubbell and S. M. Sel tzer, "Table s of X-ray mass attenuation coefficients and mass energy-absorption coefficients 1 keV to 20 MeV for e lements Z = 1 to 9 2 and 48 additional substances of dosimetric interest," National Inst. of Standards and Technology-PL , Gaithersburg, MD (United State s). Ionizing R adiation Div. 1995. [41] W. Wu, F. Liu, Y. Zhang, Q. Wang, and H. Yu, "Non-local Low- rank Cube-based Tensor Factorization for Spectral CT Reconstructio n," arXiv preprint arXiv:1 807.10610, 2018. [42] C. Badea, S. Johnston, Y. Qi, K. Ghaghada, and G. Johnson, "Dual- energy micro-CT imaging for differentiation of i odine-and gold-based nanoparticle s," in Medical Imaging 2011: Physics of Medical Imaging , 2011, p. 79611X. [43] S. G. Mallat a nd Z. Zhang, "Matching pursuits with time-frequency dictionaries," IEEE Transactions on sign al processing, vol. 41, pp. 3397-3415, 1993. [44] S. Chen, S. A. Billings, and W. Luo, "Orthogonal least squares methods and their application to no n-linear syste m identification," International Journal of contro l, vol. 50, pp. 1873-1896, 1989. [45] M. A haron, M. Elad , and A. Bruckstein, " K-SVD: A n algorithm for designing overcomple te dictionaries for sparse representation," IEEE Transactions on si gnal process ing, vol. 54, p. 4311, 20 06. [46] J. Mairal, F. Bach, J. Ponce , a nd G. Sapiro, "Online learning for matrix factorization and sparse coding," Journal of Machine Learning Research , vol. 11, pp. 19-60, 2010. [47] E. Michael and A . Michal, "Image d enoising via sparse and redundant re presentations ove r learned dictionaries," IEEE Tip, vol. 15, pp. 3736- 3745, 2006.