On the Fusion of Compton Scatter and Attenuation Data for Limited-view X-ray Tomographic Applications

by   Hamideh Rezaee, et al.

In this paper we demonstrate the utility of fusing energy-resolved observations of Compton scattered photons with traditional attenuation data for the joint recovery of mass density and photoelectric absorption in the context of limited view tomographic imaging applications. We begin with the development of a physical and associated numerical model for the Compton scatter process. Using this model, we propose a variational approach recovering these two material properties. In addition to the typical data-fidelity terms, the optimization functional includes regularization for both the mass density and photoelectric coefficients. We consider a novel edge-preserving method in the case of mass density. To aid in the recovery of the photoelectric information, we draw on our recent method in r15 and employ a non-local regularization scheme that builds on the fact that mass density is more stably imaged. Simulation results demonstrate clear advantages associated with the use of both scattered photon data and energy resolved information in mapping the two material properties of interest. Specifically, comparing images obtained using only conventional attenuation data with those where we employ only Compton scatter photons and images formed from the combination of the two, shows that taking advantage of both types of data for reconstruction provides far more accurate results.


page 13

page 14

page 15

page 16

page 17

page 18

page 19

page 20


Material-separating regularizer for multi-energy X-ray tomography

Dual-energy X-ray tomography is considered in a context where the target...

Field model for complex ionic fluids

In this paper, we consider the field model for complex ionic fluids with...

Field model for complex ionic fluids: analytical properties and numerical investigation

In this paper, we consider the field model for complex ionic fluids with...

Adaptive Regularization in Convex Composite Optimization for Variational Imaging Problems

We propose an adaptive regularization scheme in a variational framework ...

Analysis of a fully discrete approximation for the classical Keller–Segel model: lower and a priori bounds

This paper is devoted to constructing approximate solutions for the clas...

Variational models for joint subsampling and reconstruction of turbulence-degraded images

Turbulence-degraded image frames are distorted by both turbulent deforma...

PGNets: Planet mass prediction using convolutional neural networks for radio continuum observations of protoplanetary disks

We developed Convolutional Neural Networks (CNNs) to rapidly and directl...

I Introduction

X-ray CT has been used widely in fields ranging from medical imaging [2] and non-destructive evaluation [3] to the investigation of the internal structures of geo-materials [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 multi-energy 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 tensor-based 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 tensor-based iterative algorithm which reconstructs spectral attenuation images is introduced in

[14]. There, a multi-linear image model and tensor-based 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 low-dose 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 energy-integrating projections. To reconstruct attenuation coefficient images in the different channels, a patch-based 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 edge-correlation 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 non-local mean (NLM) patch-based 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 photon-counting 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 limited-view 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 maximum-likelihood model employing patch-based 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 switching-based 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 low-rank 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 SVD-matrix minimization using ADMM. In [22]

an iterative algebraic reconstruction method is proposed for sparse-view 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 Radon-transform-like 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 ad-hoc pre-processing 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 X-ray 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 energy-resolved systems, only conventional attenuation data has been considered while the majority of Compton scatter-based 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 energy-resolved detectors and limited-view geometries to reconstruct both density and photoelectric coefficients. We consider a two-dimensional 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 multi-scale method is developed to provide an initial estimate of density. An edge-preserving method developed in [36] is employed to regularize the recovery of the density. In order to stabilize photoelectric reconstruction we apply a NLM batch-based 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 limited-view 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 gradient-based and edge-preserving regularization for density reconstruction and NLM patch-based 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 Levenberg-Marquardt 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 X-ray sources are collimated to produce pencil beams that illuminate the region of interest, and that these sources are rotated step-wise in angle to produce a set of X-ray beams. Two such beams are shown in Fig. 1. We refer to an X-ray 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 attenuation-based X-ray 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.

Fig. 1: Setup of the sources and detectors. A ray from source to primary detector is scattered with angle at the interaction point and is absorbed by secondary detector .
Fig. 2: X-ray energy spectrum of a pencil beam source. The axis shows the number of photons and the axis shows the energy levels varying from to within 1 energy bins.

Ii-a Attenuation Tomography Model

For a given primary raypath the total attenuated beam intensity is calculated at each detector as [17]


where is the intensity of the X-ray 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 energy-discriminating detectors, the energy integral in (1) is over the energy bandwidth of a particular energy channel of the detector, while for traditional energy-integrating 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]


where is the mass density, is the Avogadro number, and are the atomic number and atomic weight, is the photoelectric coefficient, and , the Klein-Nishina cross section is


where . The ratio can be approximated to for most of the elements [35]; therefore (2) can be summarized as


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



is the lexicographically ordered vector of attenuation coefficients at energy level


Referring to (2), the terms that depend on energy Klein-Nishina 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

Fig. 3: Energy dependent coefficients in mass attenuation. Comparison of Klein-Nishina cross section coefficient with . The vertical grid shows the bins over which the detectors in this study aggregate photons.

from which we obtain the following model which is linear in the unknowns of interest:


where . After substituting given by (2), a set of linear equations with respect to density and photoelectric coefficients is obtained as


where is the discretized attenuation-density system matrix obtained from the terms , is the discretized attenuation-photoelectric 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 .

Ii-B Scattering Tomography Model

Again referring to Fig. 1. in this paper, the Compton scattering model captures three physical processes [38]

  1. X-ray attenuation from the source to the interaction point along the line connecting and r.

  2. Compton scattering at the interaction point r.

  3. 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]



  • 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 II-A, attenuation of the beam intensity along the line connecting and r is a function of absorption coefficient and takes the form


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]:


where is the solid angle subtended by detector . In the case of rectangular detectors, is given by [40]


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.

Fig. 4: Two rectangular detectors placed in two different locations centered on the plane where we assume all scattering is taking place. A ray emitted by the source is scattered to different secondary detectors and . The height and width of the detector, the distance from the interaction point to detector and relative angles , and determine the solid angle for each interaction point-detector pair.

To describe the scattering factor requires a bit of background regarding Compton scattering, an inelastic interaction in which an incident X-ray 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


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


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]


where is the electron density and is the differential Klein-Nishina cross section which gives the fraction of the X-ray energy scattered at angle as


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


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.

Fig. 5: Along the primary raypath which is from source to detector , is the midpoint of the -th pixel, is the midpoint of the line segment along the primary raypath, and is the length of the line segment.

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


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


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


where, with defined by (13),


Using standard linear algebra, (17) can be formulated as a set of equations non-linear in the photoelectric coefficient and quasi-linear in density resulting in a measurement model taking the form


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.

Ii-C 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



is a white Gaussian noise with zero mean and variance

. Similarly, the Compton scattering model is given by


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:


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 III-A. Subsequently, we use the current estimated density image to recover photoelectric coefficient image in another iterative process described in Section III-B.

Iii-a Density Reconstruction

With representing our estimate of the photoelectric coefficient at iteration of the algorithm, from (25), we update the density estimate by solving


where the term in (25) is not relevant as it does not depend on density.

In this paper we use an edge-preserving regularization method introduced in [36]. The approach is based on solving a series of traditional Tikhonov-type 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 peer-reviewed literature to date, we provide an overview here. To begin, recall the conventional Tikhonov smoothness-based regularization approach defined as


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


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 de-emphasizes the smoothing for locations in the image where edges are suspected. More specifically, at iteration the regularization term takes the form


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 non-negative 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 entropy-type of functional [47, 48, 49]. In the event that the Boltzman entropy is used for regularizing and if one were to employ a Bregman-type of iteration for estimating then (26) takes the form


where is the data fidelty terms in (26) and

is the generalized Kullback-Leibler divergence defined for non-negative 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 KL-based 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


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 edge-preservation, 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):


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 ),


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:


which we write in the more convenient form




There remain two issues concerning this approach: how to solve (35) and how to terminate the iteration in . The quasi-linear 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” edge-preserving 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


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.

            •  , , p and L
            • , , and
            •  and % EPI = Edge Preserving Iteration
            •  and
            •  vector of to force at least one edge-preserving 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 :
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
TABLE I: Pseudo code for iterative quasi-linear solver

The pseudo-code 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 ill-posed problem. Thus to improve the convergence rate of the density reconstruction and reduce the overall time complexity, we are motivated to consider an alternate, multi-scale 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 multi-scale process continues until we reach the desired, finest scale.

Iii-B Photoelectric Reconstruction

Given , the photoelectric subproblem takes the form


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 non-linear least squares optimization problem which we solved using the Levenberg-Marquardt 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 patch-based non-local 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


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


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 left-bottom corner of the scanning area. Forty-one 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 multi-scale processing method described in Section III-A, 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 111Tufts’ 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.

(a) Phantom I - Density Image
(b) Phantom II- Density Image
(c) Phantom I- Photoelectric Image
(d) Phantom II- Photoelectric Image
Fig. 6: Simulated phantoms. Density and Photoelectric (at the energy level of ) ground truth images of different objects described in Table II. (a) and (b) are plotted in the range of and (c) and (d) are plotted in the range of .
Material Density Photoelectric
Delrin 1.4 .4134
Graphite 2.23 .2177
Plexiglass 1.18 .3263
Water 1 .5439
TABLE II: density and photoelectric coefficient of objects in simulated phantoms.

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 signal-to-noise (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 III-A, 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 Levenberg-Marquardt method is initialized with , where for is used.

The regularization parameters and discussed in Section III-A and Section III-B 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 two-dimensional 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 multi-scale 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


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]


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 fixed-point iteration and for the edge-preserving 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


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 multi-scale 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 attenuation-only, scattering-only 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 attenuation-only data are shown in Fig. 7 and Fig. 8 for the first and second phantoms respectively. These images indicate that attenuation-only 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.

(a) Scale 1: RMSE=0.6153
(b) Scale 2: RMSE=0.4775
(c) Scale 3: RMSE=0.4043
(d) Scale 4: RMSE=0.3962
(e) Scale 5: RMSE=0.3172
Fig. 7: Density reconstruction results obtained using attenuation data alone for Phantom-I at each scale of processing at the first iteration of the algorithm. While the amplitude of the reconstructions are reasonably accurate, geometric structure is less well resolved. Subplots (a)-(e) show the density reconstruction results for different grid sizes from to .
(a) Scale 1: RMSE=0.7319
(b) Scale 2: RMSE=0.4618
(c) Scale 3: RMSE=0.4220
(d) Scale 4: RMSE=0.3209
(e) Scale 5: RMSE=0.3137
Fig. 8: Density reconstruction results using attenuation data alone for the Phantom-II obtained at each scale of processing for the first iteration of the algorithm. Subplots (a)-(e) show the density reconstruction results for different grid sizes from to .

Density reconstruction images and RMSE with scattering-only data are demonstrated in Fig. 9 and Fig. 10. Scattering-only 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 attenuation-only 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 scattering-only data and RMSE at each scale is smaller relative to attenuation-only data case. From Fig. 8(a)-(e), for the second phantom the attenuation-only 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 attenuation-only 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.

(a) Scale 1: RMSE=0.7524