I Introduction
Xray CT has been used widely in fields ranging from medical imaging [2] and nondestructive evaluation [3] to the investigation of the internal structures of geomaterials [4] and luggage screening [5], the application of specific interest in this paper. Motivated by a desire to construct spatial maps of materials properties (in our case, mass density and photoelectric absorption) in these applications, dual and multienergy CT acquisition systems [6] have drawn much attention in recent years due to their ability to provide high quality images and enhanced material characterization. Specifically, the results in e.g., [7], [8], [9] suggest that energy resolving systems perform more robustly for material characterization. For example, in the context of medical imaging, a comparative evaluation performed by [10] between spectral CT and conventional CT shows that spectral CT is more reliable in clinical applications in terms of image noise, CT numbers [11] and quality of reconstruction.
Despite these efforts, simultaneous reconstruction of both photoelectric absorption coefficient and mass density (or the closely related property of Compton scatter attenuation [12]) is still challenging due to the lack of sensitivity in the data to variations in the photoelectric absorption coefficient as a function of space [1]. To address this problem, a number of approaches have been considered in recent years. In [13]
a tensorbased dictionary learning method is introduced for material characterization, taking advantage of high correlation of the attenuation map image between different energy channels. In that work filtered backprojection (FBP) reconstruction is applied to obtain the training dictionary and an alternating iterative optimization approach is used for reconstruction. Another tensorbased iterative algorithm which reconstructs spectral attenuation images is introduced in
[14]. There, a multilinear image model and tensorbased regularization combined with total variation regularization is proposed to enhance the reconstruction results of low energy channels. In [15] the fact that attenuation images are highly correlated in different energy channels has again been used to improve lowdose reconstruction CT. It is assumed that a high quality reference image (RI) of the same object is known. The RI is reconstructed either from a set of normal dose images or reconstructed using energyintegrating projections. To reconstruct attenuation coefficient images in the different channels, a patchbased cost function capturing the correlation between the reference and reconstructed images is introduced and is optimized using the simultaneous algebraic reconstruction technique (SART) [16].Other approaches to stabilize the photoelectric reconstruction were introduced in [1], [17] focusing on the use of structural regularizers. In [17]
high and low energy attenuation data were collected to estimate Compton and photoelectric attenuation coefficients. In that work an edgecorrelation regularization is proposed to aid in the recovery of the photoelectric coefficient. The same data collection scenario is considered in
[1] to characterize materials in luggage screening application. There a nonlocal mean (NLM) patchbased regularization scheme was employed to stabilize the recovery of the photoelectric coefficient and an alternating direction method of multipliers (ADMM) method was used to solve the resulting variational problem. In [18] both photoelectric and Compton attenuation coefficients are recovered from attenuation data for different energy bins by applying a linear mapping function between images of different energy bins to minimize the difference between those images. Also, total variation and the mean of the spectral images are combined to improve the performance of the algorithm. Instead of replacing the conventional integrating detectors with the photoncounting detectors in the hardware domain, a software solution is introduced in [19] to exploit the information embedded in attenuation data over different energy channels. This method provides the spectral attenuation information by a sparse representation of the reconstructed image at each iteration in a framelet system.The methods cited in the previous two paragraphs focus on cases in which either full view data are provided or, at worse, a limited number of sources (and associated detectors) which fully encircle the object are available for generating data. Reconstruction of photoelectric coefficients in applications with more severely limitedview geometries is more challenging. In many security applications including luggage screening and kVp spectral CT [20] access to the object from different views are limited while material characterization remains quite critical. In [20]
a maximumlikelihood model employing patchbased regularization is proposed to estimate attenuation coefficients and to exploit the similarity between images from different energy channels. An alternating optimization approach is applied to reconstruct attenuation coefficient images for a set of kVp switchingbased sparse spectral CT experiments. In another kVp switching spectral CT application
[21], attenuation coefficient images are transformed to the Fourier domain and presented in the form of a lowrank Hankel matrix with missing elements. Taking advantage of the high correlation of spectral attenuation images and sparsity in the Fourier domain, the missing elements are recovered by applying SVDmatrix minimization using ADMM. In [22]an iterative algebraic reconstruction method is proposed for sparseview CT in medical applications which uses discrete shearlet transformation (DST) for denoising the reconstructed attenuation image at each iteration. Also the effective number of views is increased by interpolating the existing angular views at each iteration.
In this paper, we consider an alternate approach to mapping mass density and photoelectric absorption using both attenuation and Compton scatter data. It has been shown that Compton scatter tomography has several advantages over conventional CT systems in e.g., nondestructive evaluation applications [23]. Compton tomography also provides a powerful tool for materials characterization [24]. More specifically, Compton scattering is sensitive to structural and density variation within the object [25] by providing a strong contrast mechanism compared to total attenuation [26]. Most of the Compton scattering tomography reconstruction methods can be divided into analytical and numerical approaches. A comprehensive review of the analytical solutions is provided in [27]. The ideas introduced in [28] are the basis of most of the research in the analytical domain. It has been shown in [28] that the scattered beams collected by detectors located on a circular arc connecting the source to the detector, called the ‘isogonic line’, allows for a closed form reconstruction algorithm not unlike conventional filtered backprojection. In a related study, a Radontransformlike model for a rotating single source/single detector system is introduced in [29] and provides a closed form solution for recovering the electron density on the arcs passing through the source and detector for each point inside the object. Further developments in [30] show that a Chebyshev integral transform is also applicable to the arcs passing through each point inside the object, which confirms the results provided by [29]. The same idea has been employed in [31] for luggage screening applications. There it was shown that a combination of the proposed method and conventional attenuation tomography can produce a map of atomic number. However the approach is not robust to noise, necessitating the use of an adhoc preprocessing step of smoothing of the data. In [32] an analytic approach is proposed for reconstruction of electron densities of tissues for medical applications.
Although the analytical methods provide efficient, closed form solutions, they can only be applied to very specific data acquisition geometries. Alternatively, numerical methods such as those considered here provide the flexibility to robustly process data for more general systems. In terms of the numerical methods for Compton scatter tomography, most of the work has focused on recovering either the electron density or the total attenuation. A generalized Compton scattering transform that falls in the first category was proposed in [33] to reconstruct the attenuation map of the object of interest. The energy dependency of the attenuation coefficient at the scattering point was not considered there. In [34], the energy dependency of the attenuation is taken into account by approximating the attenuation as a linear function of energy. The algorithm tried to recover the total attenuation coefficient with an iterative minimization method and performed robustly in the presence of noise. The linear approximation to the attenuation holds in the cases that the range of energy change is small. One of the few studies seeking to recover the electron density combines three different interactions, namely fluorescence, Compton scatter and absorption [35] to directly estimate the unknown fluorescence attenuation map using Compton scattering measurements. Another approach in Xray Compton tomography assumes the attenuation coefficient is known a priori, from a traditional CT scan, resulting in a linear mapping from density to observations [26]. Most of the research performed in Compton scattering tomography has focused on gathering the scatter data on energy integrating detectors. In recent years new detectors with good energy resolution have been developed, and a valuable contribution of this paper is exploring how those capabilities can be used in the context of Compton scatter tomography.
In most of the work performed in the context of energyresolved systems, only conventional attenuation data has been considered while the majority of Compton scatterbased imaging has focused on the recovery of attenuation coefficients. In this paper we propose an inversion scheme considering both Compton scattering and conventional attenuation data for applications with energyresolved detectors and limitedview geometries to reconstruct both density and photoelectric coefficients. We consider a twodimensional form of the problem in which scattered photons are collected along with conventional attenuation measurements. A cyclic descent approach is used where we alternate between estimating spatial maps of density and photoelectric attenuation. A multiscale method is developed to provide an initial estimate of density. An edgepreserving method developed in [36] is employed to regularize the recovery of the density. In order to stabilize photoelectric reconstruction we apply a NLM batchbased regularization [1]. To evaluate the performance of the proposed method we produce several synthetic phantoms. The simulation results suggests that including Compton scattering tomography as another source of information along with conventional attenuation data can significantly enhance materials characterization especially in challenging applications with limited view geometries.
The remainder of this paper is organized as follows. In Section II, we define a limitedview system and introduce the models we use for both energy resolved attenuation and Compton scatter data. In Section III we describe the optimization problem and the iterative reconstruction method for density and photoelectric coefficients. Also, the gradientbased and edgepreserving regularization for density reconstruction and NLM patchbased regularization for photoelectric reconstruction is described. In Section IV simulation results are presented and discussed. Section V provides concluding remarks and future directions. Finally in the Appendix, we elaborate on the derivative of the cost function and regularization terms required in LevenbergMarquardt optimization method.
Ii Problem Formulation
As illustrated in Fig. 1, here we consider the recovery of mass density and photoelectric absorption in a plane (i.e., a two dimensional problem) based on attenuation and Compton scatter data. While a 2D physical model for attenuation is commonly employed, Compton scattering is an inherently three dimensional process in that even for strictly “planar” objects, photons will be scattered into the third dimension. As discussed below, the model we develop accounts for this process and provides an accurate approach for modeling the 2D problem.
Shown in Fig. 1 are two types of raypaths and detectors which will be used repeatedly in the rest of the paper. We assume that Xray sources are collimated to produce pencil beams that illuminate the region of interest, and that these sources are rotated stepwise in angle to produce a set of Xray beams. Two such beams are shown in Fig. 1. We refer to an Xray pencil beam produced by a source traveling through the object on a straight line to a detector as a primary raypath and the associated detector a primary detector. The number of sources and detectors, and , determine , the total number of primary raypaths over which attenuation data will be collected. We note that the attenuation data collected along these primary raypaths constitute a typical data set for attenuationbased Xray imaging methods.
The Compton scatter data we use for reconstruction as generated by scattering of pencil beam photons at “interaction” points along the primary raypath passing through the object. At each interaction point along the primary raypath the beam scattered along the secondary raypath is observed by a secondary detector. As photons travel from the scattering point to the secondary detector, they are further attenuated. For each primary raypath the total attenuated beam intensity caused by scattering is calculated for each secondary detector , as shown in Fig. 1. Thus, Fig. 1 illustrates the secondary raypaths connecting two interaction points along the primary raypath to a detector at . At later beam positions, the same detector will measure (as a separate observation) scattering from interaction points along those beams, for example the path illustrated in Fig. 1. We assume that sources are capable of producing pencil beams which operates over a continuous range of energies. The source energy spectrum obtained from [17] is shown in Fig. 2. We consider detectors of finite energy resolution so that data are retained only in a band of width around each at detectors. Given this general system setup, we discuss the forward models associated with both absorption and scattering data in the following sections.
Iia Attenuation Tomography Model
For a given primary raypath the total attenuated beam intensity is calculated at each detector as [17]
(1) 
where is the intensity of the Xray source at energy , is a Dirac delta function supported along the primary raypath connecting the source position to the detector located at , and is the absorption coefficient at energy . In the case of energydiscriminating detectors, the energy integral in (1) is over the energy bandwidth of a particular energy channel of the detector, while for traditional energyintegrating detectors it is over all energy. As stated earlier the goal of this problem is material characterization which requires in our case recovery of mass density and photoelectric absorption coefficient which are related to according to [37]
(2) 
where is the mass density, is the Avogadro number, and are the atomic number and atomic weight, is the photoelectric coefficient, and , the KleinNishina cross section is
(3) 
where . The ratio can be approximated to for most of the elements [35]; therefore (2) can be summarized as
(4) 
In the event that detectors are perfectly energy resolving, the polychromatic projection can be replaced by a monochromatic projection so attenuated intensity given in (1) can be reduced to a collection of linear systems (one system per energy) relating data to the unknown density and photoelectric absorption coefficient [14]. For the problem of interest in this paper however, we consider detectors of finite energy resolution. For the imaging method considered in Section III, a linear model for attenuation is rather convenient. Toward that end, we consider the following discretized model for the attenuation data which exploits the fact that the energy dependence of the coefficients in (2) are well approximated as constant over the “bins” seen by the detectors even if varies more rapidly.
To discretize the attenuation model, we assume that the object area is discretized on a Cartesian grid with elements as shown in Fig. 1. The system matrix A is then defined where represents the length of that segment of primary raypath passing through pixel and is the th row of A. The size of A is given as , the product of the number of primary raypaths and number of pixels. For each primary raypath with detector energy bin and bandwidth of , the discrete equivalent to (1) is
Referring to (2), the terms that depend on energy KleinNishina cross section and are plotted as functions of energy in Fig. 3. Two characteristics of these graphs are important to us. First, is much smaller than which implies that the data are much less sensitive to photoelectric variations than those of density, a fact we shall exploit in Section III when we discuss the imaging algorithm. Second, both of the functions vary little over the windows (shown by the vertical lines in Fig. 3) over which the detectors in this study integrate energy. Thus we replace with so that the term can be factored out of the energy sum. Now, (5) simplifies to
(6) 
from which we obtain the following model which is linear in the unknowns of interest:
(7) 
where . After substituting given by (2), a set of linear equations with respect to density and photoelectric coefficients is obtained as
(8) 
where is the discretized attenuationdensity system matrix obtained from the terms , is the discretized attenuationphotoelectric system matrix defined by , for and , and and p are lexicographically ordered vectors of density and photoelectric images respectively. The vector consists of all of the observed attenuation data as a function of source location, primary detector location and energy. The number of elements in is equal to , the product of the number of primary raypaths and energy bins .
IiB Scattering Tomography Model
Again referring to Fig. 1. in this paper, the Compton scattering model captures three physical processes [38]

Xray attenuation from the source to the interaction point along the line connecting and r.

Compton scattering at the interaction point r.

Attenuation from the interaction point to the secondary detector along the line connecting r and .
Mathematically, we capture these three processes using the following model [35]
(9) 
where

is the attenuation of the beam intensity at energy along the line connecting and r.

is the attenuation of the scattered beam at energy . We describe below the relationship between , the energy of the photon emerging from the scattering event, and , the initial energy of the photon.

is the mass density at the interaction point.

is the scattering factor. We discuss below the relationship between the incident energy of the photon, , and the scattering angle, .
As in Section IIA, attenuation of the beam intensity along the line connecting and r is a function of absorption coefficient and takes the form
(10) 
The attenuation of the beam from the interaction point to the secondary detector is much the same except for the fact that the Compton scatter process is inherently three dimensional; i.e., photons are typically removed from the plane of scattering [39]. To capture this effect, we must be a bit more careful with our modeling of the detectors. Specifically, as shown in Fig. 4, we ascribe to each detector a height and width. Only those photons scattered within the solid angle subtended by the detector are in fact observed [40]. With this, the attenuation of the beam along the line connected r and is [35]:
(11) 
where is the solid angle subtended by detector . In the case of rectangular detectors, is given by [40]
(12) 
where , and are angles defined in Fig. 4 for two different secondary detectors and , and are height and width of a rectangular detector respectively and is the distance from the interaction point to the center of detector area.
To describe the scattering factor requires a bit of background regarding Compton scattering, an inelastic interaction in which an incident Xray photon transfers a portion of its energy to a bound electron of the material being probed and emerges at an angle with respect to the initial direction. As described in [41] the relationship among the energy of the incident photon, , the energy of the scattered photon, , and is
(13) 
where is the incident energy, and referring to Fig. 1, the scattering angle which can be calculated based on the position of sources and detectors via
(14) 
In (14), is the vector from the interaction point r to the detector located at and similarly for .
The scattering factor, , in (9) is [38]
(15) 
where is the electron density and is the differential KleinNishina cross section which gives the fraction of the Xray energy scattered at angle as
(16) 
where is the electron radius.
To discretize (9), we approximate the integral over energy using a Riemann sum and employ the same grid used in the case of attenuation data to approximate all spatial integrals to arrive at
(17) 
where is the location of th secondary detector , is the midpoint of the th pixel on the primary raypath , is the midpoint of the line segment along the primary raypath through this pixel, and is the length of the line segment along primary raypath crossing this pixel as illustrated in Fig. 5.
Because one of the goals in terms of the imaging is the recovery of the mass density, and different rays cross the same pixel in many ways, we assume that the mass density is constant within each pixel and we make the distinction between and so that only one unknown will be associated with each pixel but we still provide the correct geometry (specifically, the correct scattering angles) in the specification of the forward model. In this discrete model, attenuation due to absorption between two points is
(18) 
where is a row vector of length whose entries correspond to the length of the line segments crossing each pixel on the path from point to and is the lexicographically ordered vector of attenuation coefficients at energy . Finally, the discrete form of the scattering coefficient is related to the scattering angle and the initial energy such that
(19) 
where can be computed using (14).
The inelastic nature of Compton interactions imply that even monochromatic sources will give rise to observed scatter across a band of energies thereby significantly complicating the “bookkeeping” associated with this model. Because the scattering angle is a function of the energy after the Compton event, the finite bandwidth of our detectors requires that we introduce a window factor in the definition of the scattering coefficient defined in (19) as follows
(20) 
where, with defined by (13),
(21) 
Using standard linear algebra, (17) can be formulated as a set of equations nonlinear in the photoelectric coefficient and quasilinear in density resulting in a measurement model taking the form
(22) 
where is the discretized scattering system matrix obtained from the terms , and in (17). The vector is comprised of all of the observed scattered data as a function of source location, secondary detector location, and energy. is of size where is the number of secondary raypaths, computed from the number of sources and the number of detectors , and is the number of detector energy bins.
IiC Measurement Noise
While in principle a Poisson model is appropriate for describing both the attenuation and scattered data [42], we seek to focus initially on what can be learned from this new class of data in severely limited view geometries. Thus we assume here that the only uncertainty in the data arises from typical additive, white Gaussian noise [43], [44]. We leave it to future efforts to extend the ideas developed in this paper to the more complex, but very relevant and interesting, Poisson case. More specifically, the attenuation model after adding noise is defined by
(23) 
where
is a white Gaussian noise with zero mean and variance
. Similarly, the Compton scattering model is given by(24) 
where is a white Gaussian noise with zero mean and variance .
Iii Imaging Approach
We propose the following variational problem as the basis for the recovery of density and the photoelectric attenuation coefficient:
(25) 
where measures the mismatch between the scattering data and our prediction of the scattering data for a given and p, and measures the mismatch between the attenuation data and predicted data. The regularization terms and for density and photoelectric respectively stabilize the reconstruction by imposing prior information such as smoothness, and and are weighting factors. Following [45] we set and to basically normalize the impact of the two data sets in the reconstruction process.
We employ a cyclic coordinate descent method [46] for solving the optimization problem given in (25). At each iteration, density reconstruction is performed using the estimate of the photoelectric coefficient from the previous iteration. Density reconstruction itself is an iterative procedure detailed below in Section IIIA. Subsequently, we use the current estimated density image to recover photoelectric coefficient image in another iterative process described in Section IIIB.
Iiia Density Reconstruction
With representing our estimate of the photoelectric coefficient at iteration of the algorithm, from (25), we update the density estimate by solving
(26) 
where the term in (25) is not relevant as it does not depend on density.
In this paper we use an edgepreserving regularization method introduced in [36]. The approach is based on solving a series of traditional Tikhonovtype smoothness problems where at each iteration, an evolving set of weights is used to decrease the smoothness penalty in regions where edges are suspected. As this method has not appeared in the peerreviewed literature to date, we provide an overview here. To begin, recall the conventional Tikhonov smoothnessbased regularization approach defined as
(27) 
where is the regularization parameter which determines the balance between data mismatch and regularization terms, and L is a discrete gradient matrix including both vertical and horizontal derivatives computed as
(28) 
where I is an identity matrix (assuming we are reconstructing images containing pixels), is the Kronecker tensor product operator and is the first difference matrix with on the main diagonal and on the first upper diagonal.
As noted above, here we employ an approach based on a weighted Tikhonov regularizer for which (26) is solved repeatedly. From one iteration to the next the regularization is updated in a manner that deemphasizes the smoothing for locations in the image where edges are suspected. More specifically, at iteration the regularization term takes the form
(29) 
where is the regularization parameter, is a diagonal weighting matrix with elements between zero and one, , and we call the weighted gradient of . Those diagonal elements closer to one will enforce smoothness across the associated pixels while the values closer to zero indicate that those pixels belong to an edge and should be preserved.
To motivate our choice of , consider a problem like (26) where now we wish to estimate both and . As the elements of are nonnegative and we expect that most will be close to one and a few closer to zero (since edges are sparse), a reasonable approach for regularizing these quantities would be to employ a entropytype of functional [47, 48, 49]. In the event that the Boltzman entropy is used for regularizing and if one were to employ a Bregmantype of iteration for estimating then (26) takes the form
(30) 
where is the data fidelty terms in (26) and
is the generalized KullbackLeibler divergence defined for nonnegative vectors
x and y as where e.g. is the th element of x [50]. (See for example, [51] for the relationship between Boltmzman regularization and a KLbased Bregman problem.) Using an alternating minimization method for solving (30) gives a problem similar in structure to (26) for updating the density while a closed form solution for d is easily shown to be(31) 
That is, the new estimate for each weight is a scaled version of the old weight where the scale factor is a decreasing function of the strength of the edge at that location.
Though the ideas in the previous paragraph may be potentially useful in and of themselves for edgepreservation, the exponential dependence yields an approach which is not especially sensitive to edges of varying magnitude. To achieve such sensitivity, we propose the following iteration to replace (31):
(32) 
where is a monotonically decreasing function of its argument whose range is between zero and one. In this paper we specifically take . While we leave the detailed analysis of this method to future work, as partial justification note that in the idealized case where we know the true at every iteration, the update rule (32) gives (for , the number of nonzero elements in ),
(33) 
In other words, the vector d acts as an edge detector which, in this case, is zero wherever the gradient is nonzero and one otherwise as is required by an adaptive smoother. Moreover, as demonstrated and discussed in greater depth in [36], the evolution of in this case puts zeros at locations of large edges in the earlier stages of the iteration while smaller edges are better recovered as grows. In a sense then, the approach identifies somewhat coarser structure first and then evolves to recover finer scale details.
We incorporate this approach to regularization into our recovery of by replacing (26) with the following:
(34) 
which we write in the more convenient form
(35) 
with
(36) 
There remain two issues concerning this approach: how to solve (35) and how to terminate the iteration in . The quasilinear form of the cost function in (35) immediately suggests a fixed point iteration. Specifically, starting with an initial guess for the density, call it , we build so that the resulting problem, is a linear least squares problem for . Due to the size and sparsity of the matrices comprising , the iterative solver LSQR [52] is used to find the solution to this problem. That solution is then used to build a new and the process repeats. For the problems considered in Section IV, this “inner” iteration converges rather quickly with the norm of the difference between the density estimates below in roughly iterations.
The termination of the “outer” edgepreserving iteration over is required due to the monotonically decreasing nature of the diagonal weighting matrix implied by (32). Indeed, except in cases where the gradient is exactly zero, the weights will, as , go to zero resulting in an unregularized problem. In this paper, we choose to stop when the change in weighted gradient is small or we have exceeded some maximum number of iterations; i.e., when
(37) 
where is a small number and is the maximum number of iterations. For the cases in Section IV, typically we see convergence after approximately iterations.
Inputs: 
• , , p and L 
• , , and 
Initialize: 
• and % EPI = Edge Preserving Iteration 
• and 
• vector of to force at least one edgepreserving iteration 
• 
1: While true 
2: Set 
3: Set % FPI = Fixed Point Iteration 
4: While 
5: Build according to (36) 
6: Find by solving (35) with LSQR 
7: IF : % The inner, fixed point iteration has converged 
8: Update 
9: Set 
10: ELSE : 
11: 
12: end 
13: IF or : % The outer, edge preserving iteration has converged 
14: Update 
15: Set 
16: ELSE : 
17: Update d according to (32) 
18: Increase 
19: end 
The pseudocode in Table I summarizes the overall approach for determining . To begin the process, we require an initial estimate for the density, and assume at iteration . Starting with an appropriate initial guess for the density is crucial to the success of the approach. We note there are a number of ways this could be accomplished. For example, attenuation based CT images have been shown to be useful in this regard [29]. However for the limited view problems that interest most in this effort, reconstruction of the photoelectric and density from attenuation data is known to be a highly illposed problem. Thus to improve the convergence rate of the density reconstruction and reduce the overall time complexity, we are motivated to consider an alternate, multiscale approach which is used only at when we have essentially no prior information regarding the composition of the medium. Specifically, we begin with a coarse spatial representation of the density initialized to a constant value with the same constant used for all experiments in Section IV. The method in Table I is used to solve the problem at this spatial scale and the estimated density image at this level is “upscaled” employing nearest neighbor interpolation with the Matlab function ‘imresize()’ and used as an initial guess to build the system matrix at the next finer scale. This multiscale process continues until we reach the desired, finest scale.
IiiB Photoelectric Reconstruction
Given , the photoelectric subproblem takes the form
(38) 
where is the final estimate of density image at previous iteration as a solution to (26) and is the photoelectric regularization term. In contrast to the density problem, photoelectric recovery is a nonlinear least squares optimization problem which we solved using the LevenbergMarquardt method [53]. The approach requires the Jacobian matrix of the objective function which is given in Appendix A.
It is well known that the recovery of the photoelectric map is a challenging problem [1, 17] while density is, roughly speaking, far easier to obtain accurately. To stabilize the photoelectric problem, we have used patchbased nonlocal mean (NLM) regularization method [1] which benefits from the accuracy with which density can be recovered. In this approach the photoelectric reconstructed image is conditioned on a reference image which we take as , the density estimate obtained after the first iteration of the algorithm. Mathematically, the NLM regularization can be written in the form of quadratic regularization as
(39) 
where I is the identity matrix, W is the weight matrix which is calculated based on the reference image [54], [55] and is the regularization parameter. By stacking , and vectors (38) takes the form
(40) 
The reader is referred to the Appendix for further details of the solution procedure.
Iv Experiment
To evaluate our proposed method we consider a limited view system of the form provided in Fig. 1. The area to be imaged is taken to be . Three rotating pencil beam sources each with a spectrum shown in Fig. 2 are located exactly in the center of the left and bottom edges and leftbottom corner of the scanning area. Fortyone detectors with the width and height of are equally spaced along the top and right edges. All data are generated assuming a uniform grid of pixels covering the region. For the multiscale processing method described in Section IIIA, five uniform grids of , , …, are employed for the unknown mass density. We have generated synthetic data for two different phantoms consisting of different materials with moderate to high attenuation properties shown in Fig. 6. The first phantom in the shape of an elephant ^{1}^{1}1Tufts’ official mascot is Jumbo the elephant: https://www.tufts.edu/about/jumbo and with the material properties of plexiglass provides an interesting challenge in terms of recovering the intricate geometry of the object due to some rather challenging geometric details (e.g., the space between the legs, the trunk, etc). The second phantom is more complicated with three circular objects consisting of water, Delrin and graphite. The characteristics of the materials used in these phantoms are taken from the XCOM database [56] and described in detail in Table II.
Material  Density  Photoelectric 

Delrin  1.4  .4134 
Graphite  2.23  .2177 
Plexiglass  1.18  .3263 
Water  1  .5439 
Attenuation data is collected in the range of on the energy resolution of for density and photoelectric coefficient reconstruction according to (23). Because of the size of the resulting data set ( primary ray paths scatter detectors per raypath energy bins observations), we have chosen to bin the scattered data into intervals so as to reduce the computational overhead of the processing. To consider measurement and discretization noise, a signaltonoise (SNR) ratio of dB is assumed for both attenuation and scattering measured data.
All the simulations are performed in MATLAB with the processing architecture of core Intel CPU and gigabytes of memory. The code used in these experiments is not optimized in terms of time and complexity efficiencies. The main computational load belongs to the LSQR solver and calculating forward model and Jacobin matrices, with , and on average per iteration respectively.
In the cyclic descent method described in Section IIIA, at each iteration, to reconstruct density the estimates of the photoelectric coefficient and density from previous iteration are required. At the initial iteration, , according to (26) the estimation of density requires photoelectric coefficient which we take as . The density is initialized with for both phantoms. For the photoelectric reconstruction at , is used in (38) and the LevenbergMarquardt method is initialized with , where for is used.
The regularization parameters and discussed in Section IIIA and Section IIIB are determined using the discrepancy principle [57] since the variance of the noise is assumed known. In theory these parameters should be selected by first discretizing the space of both and , calculating the reconstructions of density and photoelectric for all the points on this twodimensional discretized space, and then computing the value of the discrepancy function for each of these reconstructions. The optimal parameters and associated reconstructions output by the algorithm would be those associated with the minimum of the discrepancy function. Given the computational burden of the reconstruction process, we choose to employ the following suboptimal method. At iteration for density reconstruction where , each scale of the multiscale reconstruction process is repeated times for logarithmically spaced values of between and . At each scale, we choose that estimate of density which minimized the discrepancy function
(41) 
where is the regularized residual of the density reconstruction defined in (30), is the number of the elements of the data vector, is the regularization parameter indicator, corresponds to the scale level and is the noise variance. For , we use as the value of this parameter associated with the reconstruction selected at the finest scale of the iteration. Again at iteration , where the density estimation is used for reconstruction of photoelectric, an analagous approach is used to determine which is then used for the remainder of the iterations. Despite the suboptimal nature of this process the quantitative and qualitative measurements of the reconstruction results are highly acceptable given the limited nature of the source/detector geometry.
The stopping criteria for the overall algorithm is based on the density convergence. Thus, if the current estimation of density satisfies the convergence condition then the reconstruction of photoelectric using the final estimation of density will conclude the cyclic coordinate descent procedure. The stopping criteria is defined as [58]
(42) 
where is the density estimated vector at iteration and is a small, positive number defines the accuracy of the final results which is taken .The stopping criteria for the fixedpoint iteration and for the edgepreserving procedure defined in Table I are taken as and respectively, and .
To evaluate the performance of the proposed method quantitatively, we have calculated the relative mean square error (RMSE) for each of density and photoelectric images using
(43) 
where is the reconstruction of either the density or the photoelectric image and is the corresponding ground truth image.
There are a number of aspects of the reconstruction process we wish to explore with these examples. We first compare the recovered density and photoelectric maps after the first iteration (i.e., ) of the algorithm. This analysis allows us to explore the utility of the multiscale method for recovering density. After exploring these issues, we turn to the impact of iterating past and examine improvements seen in our ability to recover both parameters of interest. Finally, we compare our ability to quantify materials as a function of the data type used in the image formation process. We note that in all cases, the fusion of scatter data with traditional attenuation greatly improves both the quantitative as well as qualitative characteristics of the processing results.
We explore the effects of attenuationonly, scatteringonly and combination of both datasets in reconstructing density at first iteration by setting and ,then and and finally and respectively in (26).
Density reconstruction results and associated RMSE for attenuationonly data are shown in Fig. 7 and Fig. 8 for the first and second phantoms respectively. These images indicate that attenuationonly data, while providing reconstructions whose amplitudes are in the right range, suffer from significant artifacts making clear identification of the distinct regions in the scene virtually impossible.
Density reconstruction images and RMSE with scatteringonly data are demonstrated in Fig. 9 and Fig. 10. Scatteringonly data provides reconstructions where the structure of the objects are better recovered and the artifacts are reduced significantly, however the amplitudes are not completely in the right range relative to attenuationonly data. By comparing Fig. 7 (a)(e) and Fig. 9 (a)(e) for the first phantom, the shape of the elephant is better recovered in scatteringonly data and RMSE at each scale is smaller relative to attenuationonly data case. From Fig. 8(a)(e), for the second phantom the attenuationonly density reconstructions contain artifacts and noise around and along the objects so the structure of the objects are not well recovered. On the other hand, the reconstruction obtained using only scatter data contains fewer artifacts and the shape of the objects is generally better defined. We do note that the RMSE for the attenuationonly data is still smaller than that obtained using the scatter data as the absolute amplitudes of the objects are more accurate for the attenuation data even if their precise geometry is worse.