Limited-View X-Ray Tomography Combining Attenuation and Compton Scatter Data: Approach and Experimental Results

X-ray inspection systems are critical in medical, non-destructive testing, and security applications, with systems typically measuring attenuation along straight-line paths connecting sources and detectors. Computed tomography (CT) systems can provide higher-quality images than single- or dual-view systems, but the need to measure many projections leads to greater system cost and complexity. Typically, off-angle Compton scattered photons are treated as noise during tomographic inversion. We seek to maximize the image quality of limited-view systems by combining attenuation data with measurements of Compton-scattered photons, exploiting the fact that the broken-ray paths followed by scattered photons provide additional geometric sampling of the scene. We describe a single-scatter forward model for Compton-scatter data measured with energy-resolving detectors, and demonstrate a reconstruction algorithm for density that combines both attenuation and scatter measurements. The experimental results suggest that including Compton-scattered data in the reconstruction process can improve image quality for density reconstruction using limited-view systems.


I. INTRODUCTION
X-ray imaging is critical in medical [1], industrial [2] and airport security [3], [4] applications.While medical applications have attracted the greatest attention, security scanning applications such as luggage screening offer several unique challenges.Because a wide variety of materials are encountered in luggage screening, accurate material identification is important, leading to interest in using energy discriminating detectors to improve material identification [5]- [9].A second key difference is that while CT systems used in medical imaging are generally able to collect data projections at a large number of angles fully encircling the object, access to the object from multiple views is limited in many security The associate editor coordinating the review of this manuscript and approving it for publication was Yi Zhang .applications, including luggage screening and kVp spectral CT [10].There have been a variety of investigations into limited-view tomography methods, several of which exploit energy-resolved measurements and enforce similarity across energy channels [10]- [13].However, the limited-view CT problems remains quite challenging, even more so when coupled with the need for accurate material identification.
Our hypothesis is that image quality can be improved in limited-view tomography by processing not only straight-line attenuation projections through the object, but also brokenray data created through Compton scattering.Scattering occurs when photons travel from the source to a scattering object, deflect via Compton scattering to a new angle, and then are measured at a detector.Including these brokenray paths in the inversion process dramatically increases the number of geometric ray paths through the investigation domain, and thus holds the potential to reduce image artifacts.In addition, there are indications that Compton scatter data has the potential to improve materials identification as it provides a strong contrast mechanism comparing to total attenuation [14], [15].Compton scatter tomography has been explored previously and has been shown to have advantages over conventional CT systems in nondestructive evaluation applications [16] and materials characterization [17].However, processing this scattered data brings several challenges, notably the low number of counts associated with these raypaths and the computational burden of modeling the additional paths (the forward model for Compton scatter tomography models attenuation from the source to the scattering point, the scattering process itself, and attenuation of the scattered photon as it travels to the detector).In addition, the forward model could become computationally intractable in situations where multiple scattering dominates over singlescattering.
Analytic Compton scattering tomography reconstruction methods are available, but are limited to specific data acquisition geometries [18], [19].Numerical reconstruction algorithms, such as the one discussed in this paper, are applicable to more general geometries.Unlike the work presented here, many previous publications either assume that an X-ray attenuation map is known a priori from a traditional CT scan (resulting in a linear mapping from density to observations) [14], or do not fully model the energy dependence of the attenuation [20], [21].Our forward model is more closely related to that in [22], although that work includes fluorescence effects, which are not important in our application.
This work builds on and extends our recent simulationbased study that combines scattered photons with attenuation data to perform image reconstruction [23].While that work employed a similar forward model to the one considered here, the inversion approach was entirely different (as a minor note, the current model uses a different solid angle calculation which was proven to better match the experimental results).More specifically, the work in [23] focused on recovering spatial maps of both density and photoelectric coefficients.Because of the smaller problem sizes considered, a Newton type optimization method could be used.An edge preserving regularization method first developed in [24] was used to stabilize the density profile while a nonlocal means regularizer was used for photoelectric.While effective, these regularization methods are computationally intensive, requiring the solution of multiple least-square problems at every ''outer'' iteration of the Newton method.
Driven largely by the exigencies associated with processing a larger experimental data set, the effort here differs significantly from the simulation-based work in [23].Our goal is to provide an experimental validation of the benefit of combining attenuation and Compton scatter data.Thus, we focus on density inversion, deferring experimental validation of photoelectric coefficient inversion (which is more ill-posed) to future work.Specific contributions of the current work are as follows: a) we develop experimental hardware to support our studies, including design of a novel source collimator, and validate the accuracy of our forward model against experimental data; b) we employ an iterative reconstruction method that scales better to larger-scale problems than the methods from [23]; c) we reduce computational effort algorithmically, by linearizing the attenuation component of our inverse problem using a multi-energy sinogram decomposition method (see Appendix B) that exploits energy-discriminating detectors; and d) we address computational issues by demonstrating an implementation that exploits capabilities of multi-core processors.
The remainder of this paper is organized as follows.We first describe the proposed system geometry and describe the forward models for both Compton-scattered measurements and attenuation data.We then introduce a gradient descent algorithm that combines both attenuation and scatter measurements.In the Results section, we describe an experimental testbed and describe test scenarios.This testbed was used to validate our single-scatter forward model in a recent paper [25], indicating the multiple-scattering effects are negligible for this system.Finally, experimental results are shown for density reconstruction of several objects, demonstrating that incorporating scatter data into the inversion leads to noticeable improvements in image quality and accuracy of density estimates.Finally, we conclude and discuss future work.

II. FORMULATION
In this section, we describe our physical forward model, then outline the reconstruction algorithm.Fig ( 1) illustrates a two-dimensional cross section of the problem, where the slice plot is over the x − y plane and out-of-plane scattering is neglected (justified experimentally as our measurements are all in-plane).A pencil-beam X-ray source at location r S illuminates the investigation domain D inv at beam angle φ, which is the angle between r and r C that points toward center of D inv .During data collection, the angle φ is rotated to sweep across the domain.The source is moved in discrete angular steps, so in the discussion below, one beam is used to denote data collected for one source at one angle.Some photons travel on a straight-line path to a detector, providing a measurement of attenuation along the path, while others undergo non-coherent Compton scattering and are deflected to other detectors, providing broken-raypath data.For a particular beam, the detector receiving straight-line attenuation data is known as the primary detector D, while detectors receiving scattered data are known as secondary detectors D .Forward models for both types of data are now described.

A. FORWARD MODEL FOR COMPTON SCATTERING
When the incident beam crosses the investigation domain it interacts with material inside, causing non-coherent Compton scattering.Using notation similar to [22], the expected number of Compton-scattered photons collected at secondary detector r D at an energy E D is given by: where I(E S ) is the number of photons emitted by the source at energy E S .f (r, r S , E S ) computes the attenuation that occurs as the incident photons travel from source location r S to interaction point r.Similarly, h(r D , r, E ) computes the scattered beam attenuation that occurs as the scattered photons travel from interaction point r to secondary detector location r D .Notice that Compton scattering causes an energy shift, so h(r D , r, E ) depends on the scattered photon energy E .The function S(r, θ, E S ) computes the fraction of photons scattered toward angle θ at interaction point r (see Fig. 1).ρ (r) is the material density evaluated at the interaction point r.δ E D (E ) is a Dirac delta function that selects for scattered photons with energy E that matches the energy E D at the detector.Finally, δ r D ,r S is a Dirac delta function defined within D inv over a line connecting r S to the primary detector r D .We now describe the terms in Eq. ( 1) in more detail.The attenuation functions are computed as: where δ r,r S (r ) and δ r D ,r (r ) are Dirac delta functions defined along raypaths that connect r S to r and r to r D , respectively.The solid angle D for a small rectangular shape detector is approximated as [26]: where r T , r B , r L and r R are unit vectors pointing from the interaction point r to centers of top, bottom, left and right edges of the detector at r D .The attenuation function µ (r, E) in Eqs. ( 2) and ( 3) are the attenuation coefficient function, which is formulated in terms of Compton scatter and photoelectric effects as: where N A is the Avogadro number and Z (r) and A(r) are the atomic and mass numbers, respectively.ρ(r) and p(r) are the material density and photoelectric coefficients, while f p (E) = (E 0 /E) 3 is the photoelectric energy factor with E 0 as the referenced energy.An interested reader is referred to [23] for details concerning the computation of the Klein-Nishina cross section f KN (E) in Eq. ( 4), the scattered energy E and the scattering factor S(r, θ, E S ) in Eq. (1).Note that Eq. (1) computes the scattered photons that are incident on the detector at energy E D , and does not account for the fact that the detector provides an imperfect measurement of photon energy.The actual signal recorded by the detector at energy E q is found by integrating over the detector sensitivity function: where s q (E) is the sensitivity function for the q th energy channel, centered at E q .Idealized models for detector sensitivity functions are presented in [27], and model the sensitivity functions as being Gaussians whose bandwidth depends on detector properties.More realistic sensitivity functions can be computed for particular detectors and account for energydependence of the sensitivity function; we employ a model for the Multix ME100 detector found in [6].

B. DISCRETE FORM OF THE COMPTON SCATTERING MODEL
The spatial domain in Eq. 1 is discretized into equal rectangular cells, see Fig. 2, such that the dimension across x and y-axis has N x and N y segments, respectively, and the total number of cells is N Cells = N x N y .The source spectrum is discretized into N E energy levels, where the k th level is centered at E S k and has width of E S .Similarly, we discretize the photon energies incident on the detector into N D energy bins, each with bin width E D .Then, the discretized version of Eq. ( 1) can be written as: The indices i and j refer to a particular incident beam and secondary detector in the system, respectively, while the index m indexes the fine scale energy bins into which we aggregate the scattered photons prior to accounting for the finite energy resolution of the detectors in Eq. 10.The index l indicates a discretized cell (that contains the interaction point ri,l ) from the set of cells that the i th incident beam passes through.The interaction point ri,l is taken to be located at the middle of the discretized segment in the l th cell.The inner sum in Eq. ( 6) computes numerically the spatial integral in Eq. (1) (integral over r) using a Riemann sum.The two attenuation functions f and h are computed numerically using Riemann sums to approximate the integrals in Eqs. ( 2) and (3), giving where a T i,l is a row vector with elements ordered lexicographically from a matrix the size of the discretized grid.This vector stores zeros for elements corresponding to cells that do not intersect with the incident beam connecting r S,i to r i,l , and otherwise stores the length of the intersected beam segment across each cell.In the same manner, the vector b T i,j,l stores lengths along the scattered ray path from the interaction point r i,l to the secondary detector r D ,j .Moreover, µ(E S k ) and µ(E i,j,l,k ) are vectors with lexicographicallyordered elements that store the attenuation function (Eq.( 4)) evaluated over the discretized cells, assuming energies of E S k and E i,j,l,k , respectively.The subscripts on E i,j,l,k indicate that the energy shift is a function of the scattering geometry.
Referring back to Eq. ( 6), θ i,j,l in the scattering function S(r i,l, , θ i,j,l , E S k ) is the scattering angle that lies between i th incident beam and scattering beam connecting l th interaction point to j th secondary detector.On the other hand, the weighting function ω(i, j, l, k, m) is used to approximate the Dirac delta function δ E D (E ) in Eq. (1) as: The function δ i,l in Eq. ( 6) approximates δ r D ,r S (r)dr in Eq. ( 1) by storing the segment length of i th incident beam intersecting with l th interaction cell.ρ(r i,l, ) is the discretized value of the density for the cell located at ri,l, .Eq. ( 6) computes the photon counts before interaction with the detector.A discrete form of the detector sensitivity function can be found by computing a detector response matrix S (described in [6] for the Multix detectors used in our experiment).If the detector has Q output energy channels, then S is of size Q × N D and its entries are S q,m = s q (E m ).The signal g D (i, j, m) after the detector interaction occurs is computed as:

C. FORWARD MODEL FOR ATTENUATION DATA
The forward model for attenuation data is considerably simpler than the Compton model, and builds closely on the notation above.In subsequent discussion, we will also refer to attenuation data as transmission or 'Tx' data.For a raypath between a source-detector pair (r S , r D ), the continuous model for the number of counts collected on detector energy channel q is: where f (r D , r S , E S ) has the same form as defined in Eq. ( 2), but has a raypath that extends from r S to r D .All other terms are as defined above.
Similarly, the discretized version of the forward model follows closely from Eq. ( 7) and Eq.(10).Eq. ( 7) becomes for the attenuation case: where i denotes the raypath.a T i and µ(E S k ) are row vectors defined as in Eq. (7).Eq. ( 12) differs from Eq. ( 7) only in that there is no need to track the intersection point.The multiplication of a T i and µ implements the integral along the raypath as a Riemann sum.
When processing dual-or multi-energy X-ray attenuation data, a common approach is sinogram decomposition, in which the measured projections are decomposed into the two energy-independent terms shown in Eq. (4) (see [9] and references therein).In this representation, the spatial integral in Eq. ( 2) is rewritten as: where c p,i = a T i p is the integral of photoelectric absorption along beam i, where p is a column vector of size N Cells × 1 that contains the lexicographically unwrapped values of the photoelectric absorption in each cell.Similarly, c ρ,i = a T i ρ integrates a scaled version of the density along the raypath, with the scaled density defined as ρ(r) = N A 2 ρ(r), where we have taken advantage of the fact that Z (r) A(r) ≈ 1/2 for most elements [22].
This formulation of the forward model is convenient as data can be collected at various energies and used to estimate c ρ,i and c p,i , solving a separate nonlinear inverse problem for every raypath i.Details of the sinogram processing are described in Appendix B and as well as in [9].The result of this processing are estimates of the projected density and photoelectric coefficients along each raypath.These projections are related to the underlying data through a system matrix A of size (# raypaths x N Cells ) whose rows are the vectors a T i already defined: Thus, once sinogram decomposition is performed, recovery of the density and photoelectric images is a straightforward linear inverse problem.In the overall inversion method described below, we seek to solve this linear inverse problem as part of a joint inversion that combines attenuation and scattered data.

D. RECONSTRUCTION ALGORITHM
Eq. ( 10) above gives an expression for the scattered data for source i, receiver j and energy q.Collecting the results for all i, j and q into a vector g D , the forward model for the Compton scattering signal described in Sections II-A and II-B can be represented as: where f (ρ, p) is the nonlinear model of the Compton scattering signal as a function of both the density ρ and the photoelectric absorption coefficients p.To avoid confusion with the signal g D computed using the model above, the actual experimental data recorded from the instruments will be referred to as d.Similarly, in place of c ρ and c p from Eq. ( 14), s ρ and s p will refer to the sinogram decomposition of the density and the photoelectric coefficient, extracted from the actual transmission data.For density reconstruction, the optimization function can be written as: where γ sc ≤ 1 and γ tx ≤ 1 are weights used to control the relative importance of data misfit terms for scattering data (first term) and attenuation data (second term).The enforcement of the TV regularization is controlled through the λ TV parameter.Higher λ TV will encourage reconstructed profiles that are piecewise constant with clearly defined edges, a widely used regularization strategy in attenuation-only computed tomography [28]- [30].Equation.( 16) is subjected to an additional constraint ρ ≥ 0, i.e., negative densities are not allowed.
The formulation above lets us study the relative impact of scattering and attenuation data on our solutions.In this work, we will consider four scenarios to study the optimization problem in Eq. ( 16), namely i) inversion using Compton scatter data only (γ sc = 1, γ tx = 0), ii) inversion using attenuation data only (γ sc = 0, γ tx = 1), iii) inversion equally weighting both data types ( γ sc = γ tx = 1), and iv) and inversion that weights scattering data more heavily than attenuation data (γ sc = 1, γ tx < 1).The optimization function in Eq. ( 16) is written in term of ρ only, so an estimate of the photoelectric coefficients is required.Although the Compton scattering signal as described in Eq. ( 6) depends nonlinearly on both ρ and p, we can exploit the fact that the higher energy attenuation is dominated by ρ.Fig. 3 plots the two portions of the attenuation function (see Eq.( 4)) due to the Compton scattering term (N A (Z (r)/A(r))f KN (E)ρ(r) and due to the photoelectric absorption term (f p (E)p(r)) for Delrin, a plastic material used in our experimental work.The figure shows that the photoelectric absorption effect becomes increasingly negligible above 40KeV; similar results are seen for other materials as photoelectric absorption falls as the inverse cube of energy.Hence, in our density reconstructions, we neglect the photoelectric effect (set p = 0 in Eq. ( 16)) and measure data misfit for scattered data for higher energy photons only.While multiple higher-energy bins could be used, in the results below we sum all higher-energy scattered photons into a single energy channel.Note that this approach requires the use of X-ray detectors with the ability to resolve detected photon energy.
For the reconstruction process, the following algorithm shows how the optimization in Eq. ( 16) can be carried out using a steepest descent approach [25], [31], [32], [32]- [37] known as embedded nonlinear Landweber-Kaczmarz algorithm [25], [36], [37] Step 3.2 : While most of the notation above has already been introduced, we note that β is a small number used to remove the derivative singularity in the TV regularization term, N iter is the number of gradient descent iterations, and ∂ ρ i f ρ i , p 0 is the adjoint of the density differential, a matrix that is described in detail in Appendix A which relates each measurement to the change in density at every point.In the algorithm, Step 1.0 initializes the density profile ρ 0 to zero and sets the photoelectric profile p 0 ; note that the density profile is updated in Step 3, while photoelectric values remain at zero.In Step 2.0 the user sets the weighting parameters γ sc and γ tx of the Compton and transmission data misfits, respectively, and sets the TV regularization parameters.Step 3.0 is the iteration loop of the nonlinear Landweber algorithm with N iter iteration steps.Step 3.1 shows one step of the Landweber iteration over the density ρ, where γ sc ω sc i and γ tx ω tx i defines the step sizes per iteration for scattering and transmission updates, respectively.A sufficient condition for convergence of nonlinear Landweber-Kaczmarz states that γ sc ω sc i ≤ 1/(σ sc ) 2 and γ tx ω tx i ≤ 1/(σ tx ) 2 where σ sc i and σ tx are the maximum singular values of ∂ ρ i f ρ i , p 0 and (N A /2)A respectively, [36], [37].Since γ tx ≤ 1 and γ sc ≤ 1, to achieve maximum possible step sizes under the convergence condition, ω tx i = 1/(σ tx ) 2 and ω sc i = 1/(σ sc ) 2 are assumed.The upper script ( T ) is used to indicate the transpose operator.Finally, Step 3.2 enforces the nonnegativity constraint on density [38].

E. COMPUTATIONAL CONSIDERATIONS AND IMPLEMENTATION
The maximum step size for the Landweber algorithm is set by ω sc i and ω tx , which in turn depend on the maximum singular values σ sc i and σ tx for the sensitivity operators related to scattered data and attenuation data [33], [35], [37].These maximum singular values are computed using power iteration [39].While σ tx is a fixed quantity as the transmission model is a linear model, σ sc i changes each time the nonlinear Compton scattering model is updated.The computation of σ sc i is expensive as it requires evaluating the differential and the adjoint operators at least four times each.However, it is not necessary to re-evaluate σ sc i on each iteration, due to the fact that the nonlinear Landweber will have its largest changes in the first iteration steps, with changes in σ sc i becoming smaller as iteration proceeds.Hence for the results below, σ sc i will be only evaluated for selected iteration steps.For all the examples, σ sc i is computed at the first three iteration steps, then at every tenth step up to iteration 50, then again at iterations 75, 100 and 150.
The computation at each iteration is dominated by calculating the forward model Eq. ( 6), the differential ∂ ρ f (ρ, p) Eq. ( 26), and the adjoint ∂ ρ f (ρ, p) T Eq. ( 27) operators.Each of these terms is a matrix indexed by the incident source beam, i, and secondary detector, j.Computing each element of these matrices involves integrations over spatial cells and photon energies, rapidly increasing the dimensions of the problem and leading to a computational challenge.
Fortunately, these calculations (described in Appendix A) lend well to parallelization across beam/detector pairs.
Initially the reconstruction algorithm was fully implemented in Matlab 2016a and run on a High Performance Computing cluster at Tufts University.Calculation of the elements of the forward model, differential, and adjoint matrices were each parallelized by splitting the elements evenly across nodes using Matlab's Parallel Processing Toolbox.Reconstructions for the scan scenario and test images discussed in Section III.C typically required 3-5 days on up to 20 compute nodes (the number and type of cluster nodes available varied between runs).
To allow faster processing times, the forward model, differential, and adjoint calculations were translated to C and parallelized for a single node using OpenMP.The Intel icc compiler was used with O3-level optimization and AVX2 flags.Further acceleration was achieved by improving the load balancing across processors.Note from Figure 2 that different incident beams generally traverse over a different number of voxels in the investigation domain.This means that computing the integrations over voxels for different source/detector pairs require different amounts of computation.Rather than splitting the computation by assigning each processor a fixed number of operator elements, prior knowledge of the geometry allows the workload to be split more evenly thus achieving a better load balance.Running the overall reconstruction within Matlab, but escaping to the compiled routines for the three major operators, resulted in equivalent reconstructions that completed on the order of 9 minutes on a single dedicated server with two 20-core Xeon E5-2698v4 processors.

III. EXPERIMENTAL RESULTS
In this section, experimental results are illustrated to verify the forward model and to study reconstruction results.We describe the exeperimental testbed, discuss calibration and validation of the forward model, then show reconstruction results for several scenarios.

A. EXPERIMENTAL TESTBED
A multi-energy X-ray testbed system was constructed to experimentally measure both transmission (Tx) and Compton scattering data in a controlled environment using energydiscriminating detectors.This testbed, shown in Fig. 4, consisted of an L-shaped detector array mounted on an optical table (L-shaped detector arrays are the norm for conventional transmission X-ray baggage scanners, as they make efficient use of space for the rectangular tunnels used in these scanners).A roughly 60 cm × 40 cm section of the optical table was outfitted with regular 2'' × 2'' indentations that allow the repeatable positioning of material samples and other image targets.Both 2'' square and 2'' diameter circular test objects could be mounted directly to the testbed.Mechanical adaptors were built to allow for repeatable positioning of larger objects.8); markers show the origin of (x, y ) coordinate system.A slit opening built into the detectors ensures that only in-plane scatter is measured; b) zoom on coaxial target (reconstructions in Fig. 9); c) zoom on box target (reconstructions in Fig. 11).

Detectors:
The detector array consisted of ME-100 photoncounting detectors (Multix, Inc., Moirans, France), arranged in 14 modules each containing 128 detectors.Five modules (640 detectors) were arranged on one side of the L with 9 modules (1152 detectors) on the other, giving 1792 detectors in total.Each detector pixel was approximately 0.08 cm in width and height.A lead collimator is built into the detector package to suppress scatter from objects outside of the plane formed by the X-ray source and detector array.An additional radiation shield, shown in Fig. 4, further reduces scatter from outside the detector plane.During data collection, detectors were operated at the finest energy resolution, outputting 128 evenly spaced energy channels from 20-160 keV.This represents an over-sampling, as the actual energy resolution of the ME100 is closer to 5-7 keV, depending on energy.A detailed model of the ME100 energy resolution, based on [6], was used in reconstruction.
Source: The objects were interrogated using a single X-ray source (bremsstrahlung spectrum with peak energy 160 keV and average energy 62 keV, giving the spectrum in Fig. 5).A source collimator was used to produce a narrow pencil beam, as discussed in detail below.The source intensity (shown in Fig. 5) was computed pre-collimation, leading to a need to calibrate the emitted source flux, with calibration described in the next section.This source was mounted on a pivot at a radius of 76.2 cm from the origin of the coordinate system, which was taken to be in the middle of the interrogation region (see Fig. 4(b)).The testbed could be precisely configured to interrogate the scene from 6 different angles spaced apart by 22.5 • , ranging from −22.5 • to 90 • , with the 0 • source perpendicular to the longer detector array (source position (x = 0, y = −76.2cm) and the 90 • source perpendicular to the shorter detector array (source position (x = −76.2cm,y = 0).For each source position, the source was rotated so beams were swept across the interrogation domain in 0.4 • steps.Data from the various source positions were collected to mimic a multi-source system in which sources scan the region sequentially.
Initial experimental results showed that scatter from the pencil beam forming collimator itself was discovered to be a major limiting factor, which resulted in several iterations of hardware development.The original design employed a pencil beam forming system that is standard in commercial X-ray backscatter systems [40], shown in Fig. 6(a).The cone beam from the X-ray source is chopped into a pencil beam by a combination of a rotating disk with radial slot apertures and stationary pre-collimator slot in the plane of the detectors.As the disk rotates, the combination of disk aperture and pre-collimator form a rhombus shaped opening which sweeps across the field of view.Unfortunately, a rotating disk chopper requires a knife edge aperture, as shown in Figure 6(b),which created substantial scatter and fluorescence, both of which added noise to the overall data.Furthermore, because the point of scatter moves as the disk rotates, a complex set of systematic errors were embedded in the data, which included shifting transmission images of any objects in the tunnel.
To avoid these artifacts, the final design employed a snout which was fixed to the X-ray source, and the entire assembly of source and snout rotated to sweep the pencil beam across the tunnel, as shown in Fig. 6(c).The direct connection of snout to source results in containment of leakage radiation and a uniform pencil beam cross section.More importantly, this geometry controlled the scatter from the apertures themselves.Monte Carlo simulations predicted that simply replacing the knife edge aperture with a simple tunnel should reduce scatter and fluorescence by 10x-20x.The inclusion of a scatter shield further limited unwanted scatter in all directions except for a small region near the primary pencil beam.Collimation by the 'primary aperture' forms the pencil beam.The primary aperture itself is positioned only as far from the focal spot as deemed necessary to provide adequate spatial resolution (i.e. to adequately limit of the divergence angle of the pencil beam.)Most of the length of the snout serves as a scatter shield to contain scatter from the primary aperture.The exit opening of the snout defines a narrow angle of allowed scatter which is effectively a halo around the pencil beam.Direct beam from the halo illuminates a few dozen pixels on either side those in the path of the good pencil beam.Compton interactions from the narrow halo are orders of magnitude fewer than from the unwanted scatter in the original system, which was critical to obtaining a good match between forward model and experimental data.All data results shown in this paper were generated using this snout design.

B. CALIBRATION AND FORWARD MODEL VERIFICATION
The X-ray source used to generate the incident beam has intensity I S (E S ) shown in Fig. 5. Due to the source collimation, a calibration parameter α c is required to determine the effective photon counts on the incident beam I(E S ) in Eq. ( 1), such that I(E S ) = α c I S (E S ).The parameter α c was determined by matching the forward model to data collected using a reference object (a 2'' cylinder of known density) with the beam pointed at the center of the object.After collecting the data, high noise detectors were dropped and the predicted data for the remaining detectors was computed using the model described in Section II-A with I(E S ) = I S (E S ).The computed curve was then multiplied by the unknown α c and matched with the data using least squares minimization.
Example output from the resulting calibrated model is shown in Fig. 7.A Delrin (CH 2 O) target of 2.54 cm radius was placed at the center of the experimental fixture.The X-ray source was located at 0 • source (see description in Section III-A) with the beam directed toward the center of the object.Fig. 7(a) shows the sum of photon counts over all the energies at each detector for both experimental and simulation results.In the experimental results, the photons are collected over 0.1 msec time slots and are then averaged over 200 time slots for a total observation time of 20 sec.The jump at detector 1153 marks the transition from the longer to the shorter arms of the L-shaped detector array (shown in Fig. 4).Fig. 7(b) shows the number of photon counts per energy at a single detector (number 1550), for both simulation and experiment.The effect of statistical noise can be clearly seen in Fig. 7(b), as the number of photons per energy bin on one detector is relatively low.However, in general the forward model in Fig. 7 shows a close match to the experimental results.Note that spectrum shape in Fig. 7(b) is smoothed compared to the source spectrum (Fig. 5) due to the imperfect energy resolution of the detectors [6].
While the overall match in Fig. 7(a) is good, isolated detectors are seen that have anomolously high or low counts relative to their neighbors.Detectors are arranged in blocks of 32, and investigation showed that a majority of anomolous detectors were located at the edges of each detector block, while other anamolous detectors could be identified by measuring dark-scan current (i.e., detector output with the X-ray source turned off).Anomalous detectors were identified and removed before further processing.

C. RECONSTRUCTION RESULTS
To test image reconstruction, a series of phantoms objects were measured on the testbed.These phantoms range from simple geometric objects to objects more typical of airport luggage.In this section, we present reconstuction results for three scenarios: a) an arrangement of three plastic rods; b) a coaxial plastic object with an inner Delrin cylinder surrounded by an outer HDPE cylinder; and c) a cardboard box packed with cotton T-shirts that concealed Delrin and HDPE square rods.As described above, detector collimation was used to ensure that only in-plane scattering was measured, so all reconstructions are two-dimensional.
Data were collected from all 6 source locations.Because we seek to understand the benefit of Compton scattering data for few-view systems, we processed data for configurations of 3, 4 and 6 sources.The three-source configuration included the (90 • , 45 • , −22.5 • ) sources, while the 4-source configuration used the (90 • , 45 • , 0 • , −22.5 • ) sources.Reconstructed images are plotted below only for selected source configurations, but image quality metrics for a wider set of configurations are found in Table 1.
A common set of parameters were used to reconstruct all scenarios shown.Reconstructions were generated on an evenly spaced grid (4 mm grid spacing in both dimensions) for 150 iterations of the algorithm.
Total variation parameters were set as β = 10 −6 and λ TV = 0.01 for all examples.Selecting λ TV is in general problem dependent.For optimization problems with ''weighted'' misfit terms that are function of their misfits, λ TV should be updated to balance the optimization problem in accordance with the change in misfits weights [30].The later makes selecting λ TV a difficult task, so optimization schemes have been developed for such problems, e.g. the multiplicative constraint approach [38].However, if the misfit terms are multiplied with fixed weights, as in Eq. ( 16), then selecting λ TV becomes simpler and estimates can be found empirically [41], [42].For all the examples illustrated in this work, a fixed value of λ TV = 0.01 is used.Transmission sinogram decomposition was performed for all source-detector raypaths after elimination of bad detectors (as discussed in the previous section).Compton scattering measurements from good detector pixels were spatially averaged into 200 equivalent detectors.In effect, this simulates the use of larger-area scatter detectors (roughly 7x area increase), which helps to improve count statistics and to reduce computation time.

1) THREE TARGETS
In this example, the investigation domain contains three targets: a 2'' diameter PVC cylinder (C 2 H 3 Cl), a 2'' square Delrin rod (CH 2 O), and a 2'' diameter Delrin cylinder (CH 2 O) (see Fig. 4).The results of density reconstructions for a 4-source configuation are shown in Fig. 8. Fig. 8(a) shows the nominal (ground truth) density.To simplify the plotting, this and all subsequent density reconstructions are shown for a fixed range of 0 − 1.5 g/cm 3 (note that PVC and Delrin densities are quite similar).
The remaining plots show reconstructed density using Compton-scattered data only (Fig. 8(b)), Tx data only (Fig. 8(c)), and both Tx and scatter data (Fig. 8(d)).For results in Fig. 8(d), equal weighting is placed on both Tx and scatter data (i.e., λ tx = λ sc = 1).The use of Compton scatter data alone yields a reconstruction with reasonable geometric detail, but with greatly underestimated density values.As a note, this result is the solution of a nonlinear, non-convex problem without convergence guarantees.The Tx-only reconstruction in subfigure (b) shows improved density estimates, but density is still underestimated, and there is a noticeable increase in artifacts (see blue background region).The combination of Tx and Compton scatter data in (d) shows the most accurate (highest) reconstructed density values, and also shows a reduction in artifacts relative to Tx-only results.The geometric definition of the three shapes is also improved as compared to Tx-only results.This improvement is attributed to the additional sampling of the scene provided by broken-ray raypaths.
Because a reliable estimate of the actual density is availquantitative measures of the reconstructed image can be computed, using the nominal density image (subfigure (a)) as the reference image.Overall estimation error is measured using normalized mean-squared error (expressed as a percent), as well as the widely-used Structural Similarity Image Metric (SSIM) [43].SSIM is most typically used as an image quality metric that matches human perception, which has relevance for an X-ray imaging system where images would be reviewed by a human operator.Both metrics are shown in Table 1 for all phantoms and various numbers of sources.Note that the combination of Tx and scatter data gives the lowest error for all cases, and the best SSIM for all but one case.

2) COAXIAL GEOMETRY
As a second example, we consider a coaxial geometry in which a Delrin cylinder is placed inside an outer HDPE shell.An interesting feature of the Tx-only reconstruction in Fig. 9(c) is the noticeable variation in the reconstructed density of the outer ring.This is a result of the extremely limited number of angles used in reconstruction (the effect is reduced if all 6 sources are used).One potential gain from adding Compton scatter data is that the additional raypaths provide additional geometric information that can reduce artifacts in few-view reconstruction.A modest improvement is in fact visible, with the reconstructed HDPE ring being more homogenous in Fig. 9(d) than in Fig. 9(c).This leads to the question of whether Compton scatter data can be further exploited to reduce image artifacts.
As noted in Section II-D, convergence is guaranteed if the gradient descent steps based on Tx and scatter data are limited by choosing step size multipliers λ tx ≤ 1 and λ sc ≤ 1.In the results above, the largest possible steps were taken (λ tx = λ sc = 1) to accelerate convergence.However, given that the limited-view artifacts are primarily associated with transmission data, it is interesting to explore whether de-emphasizing Tx data can further reduce artifacts.To quantify this, the ground truth image in Fig. 9(a) was used to identify all pixels corresponding to the outer ring.Variability of reconstructed density for these pixels was computed using percentile measures, with results shown in Table 2.There appears to be a modest but measurable improvement in homogeneity, with the most homogeneous reconstruction given by the result plotted in Fig. 10(a).Fig. 10(b) compares the overall image error for various solutions as a function of iteration number, showing that varying the ratio λ tx /λ sc changes the shape of the error curve vs. ground truth values.While in this instance the lowest overall error is given by equal weighting, determining the optimum weighting for transmission vs. scattered data in the inversion is a potentially interesting avenue for future work.

3) PACKED BOX WITH COTTON AND TARGETS
As a final example, Fig. 11 shows a more 'real-world' scenario with more inhomogeneity in which polyethylene and Delrin objects are surrounded by several inches of cotton fabric (T-shirt material) and cellulose (the cardboard boxes used to contain the fabric.)This configuration approximates a cross section of a typical carry-on sized suitcase, packed with clothing and two masses of contraband organic material.The interior of the box is imaged, so the background density is that of cotton (taken to be 0.23 g/cm 3 ).In this scenario, the Compton scatter-only result shows very weak contrast at the location of the two objects.The Tx-only result does recover the objects, but underestimates their density.The combination of Tx and scatter data yields the most accurate density estimates for the hidden objects, as shown by the normalized error listed in Table 1.

IV. DISCUSSION AND FUTURE WORK
The work above has described physics-based models for both Compton scattering and transmission data, and used these models to reconstruct density of objects imaged by very fewview tomography systems.The experimental results showed that reconstructions that combined transmission (attenuation) data with Compton scatter measurements yield more accurate density estimates, and have reduced image artifacts, compared to transmission-only reconstructions.This suggests that current systems, which typically measure only transmission data and seek to eliminate scattered photons through collimation, are discarding potentially useful information that could improve image quality.A further result is that a first-order scattering model (neglecting potential multiple scattering) was sufficient for image recovery in our experiments.
The observed improvements in image quality can be attributed to underlying scattering physics.While Compton scattered data and attenuation measurements are both affected by the same underlying physical properties (density and photoelectric coefficient), the measurements are affected differently by these material properties.X-ray backscattering, which is a special case of Compton scattering, is known to highlight organic and low atomic number materials that have a low signature in transmission mode [40].Thus, the addition of scattering data provides new cues that can be exploited in material recovery.As discussed above, another key feature of Compton scattering is that the broken raypaths add geometric paths through the scene that can be important in few-view imaging scenarios, leading to reduction in image artifacts.
In this paper, we focused on density reconstruction only, deferring photoelectric reconstruction to future work.However, the photoelectric coefficients can be recovered using a similar Landweber iteration.Because the Compton scattering effect (proportional to density) is significant at all energies (see Fig. 3), reasonable density estimate is required for photoelectric recovery.The photoelectric profile could thus be recovered using a cyclic descent approach, in which Algorithm 1 above is used to estimate density, after which the estimated density is used to initialize the photolelectric recovery, and iteration between density and photoelectric reconstructions can be repeated as desired, yielding improved estimates on each iteration.
An important goal for future work is understanding the benefit of the proposed approach for a wider array of problems.Ideally future studies would include a wider variety of materials (textured materials, metal) as well as scenes with increased clutter.In addition, the geometric arrangement of X-ray sources and detectors could be optimized to ensure a more even coverage from all sides of the scanned object.
Multiple avenues exist to further develop X-ray systems that combine traditional attenuation-based imaging with Compton scatter data.Given the low count rates associated with Compton scattered data, larger area detectors or higherflux sources would be desirable.Algorithmically, there are open questions on how best to combine transmission and scattering data during inversion; Fig. 10 and the related discussion explored this idea briefly.In part for computational reasons, we primarily used the energy-discriminating detectors to eliminate low-energy data where photoelectric effects are significant, and summed high-energy data together into a single energy bin.However, multiple high-energy bins could be used during reconstruction, and the optimal energy binning for detector data remains to be fully explored.where N P = N E N I , N E is the number of incident energies and N I is the number of interaction points.Using Eq. ( 22) and Eq.27 gives the definition ∂ ρ f (ρ, p) The transposed operator is then computed as Computing ∂ p f (ρ, p) and ∂ p f (ρ, p) T follows similar steps as above.

APPENDIX B
In this appendix, we briefly review the sinogram processing for transmission data.Our approach is a multi-energy extension of the dual-energy sinogram decomposition method presented by Ying et al. [44].For more details, the reader is referred to [9].Based on the relation in Eq.13, we can recover c ρ,i and c p,i by minimizing the squared error between the measured projections and modeled projections.Here we define θ i = c ρ,i , c p,i as the unknown density and photoelectric projections along raypath i.Then, the log-normalized modeled projection along raypath i for source energy E s is K (θ i ) = − ln I (E s )exp(−c ρ,i f KN (E s ) −c p,i f p (E s ))dE s ] + ln I (E S )dE S (30) Ying et al. studied the dual-energy case where the object is scanned twice with low-and high-energy source spectra.They used low-and high-energy measured sinograms to estimate sinograms in Compton (related to scaled density) and photoelectric space by minimizing the squared error between the summed low-and high-energy scans projection and a modeled project computed using Eq. 30, subject to non-negativity constraints.We generalize their result to multi-energies by seeking a solution that minimizes with nonnegativity constraints: c ρ,i ≥ 0 and c i p,i ≥ 0 In the equation above, vector K(θ i ) is the N E dimensional column vector of modeled log-normalized mean values computed for different values of E s , using Eq.30 and m i is the N E dimensional column vector of log-normalized measured projections.
= diag{w} and w is the number of photon counts received in each energy bin.The matrix acts as a weighting matrix, making the problem a weighted least squares problem, and was derived by Bouman and Sauer by finding a quadratic approximation of the Poisson log-likelihood function for X-ray transmission problem using a Taylor series expansion, which can be described as [45]- [47].This weighting matrix gives more weight to the raypaths and energy bins whose measured photon count is larger.In practice, for rays passing through objects with high attenuation, this gives more weighting to the high energy bins.In the case where only two energies are present and is the identity, Eq. 31 reduces to the result in [44].For the interested reader, [9] explores cases where multi-energy sinogram decomposition has advantages over dual-energy sinogram decomposition.

FIGURE 2 .
FIGURE 2. Discretization of the investigation domain D inv .

FIGURE 4 .
FIGURE 4. Experimental setup for measurements; a) top view showing the X-ray source, source collimator, and L-shaped array, three-target scenario (reconstructions in Fig.8); markers show the origin of (x, y ) coordinate system.A slit opening built into the detectors ensures that only in-plane scatter is measured; b) zoom on coaxial target (reconstructions in Fig.9); c) zoom on box target (reconstructions in Fig.11).

FIGURE 6 .
FIGURE 6. Source collimation alternatives tested during experimentation.a) shows a conventional chopper system for creating a pencil beam; b)shows artifacts generated by scattering off the conventional chopper aperture; c) shows the rotating source and snout concept, with a shield to limit scatter from the primary aperture.

FIGURE 7 .
FIGURE 7. Forward model verification, showing (a) total photon counts per detector, and (b) photon counts per energy bin at detector 1550.Detectors 1-1152 are for the upper side of the detector (see Fig. 4), and the jump in figure (a) at detector 1153 marks the transition between sides of the L-shaped detector array.

FIGURE 8 .
FIGURE 8. Reconstructed result of three targets example, showing (a) actual density profile, (b) reconstruction using Compton scattering data only, (c) reconstruction using Tx data only, (d) reconstruction using both Compton and Tx data.All density reconstructions are shown for a dynamic range of 0 -1.5 g/cm 3 .

FIGURE 9 .
FIGURE 9. Reconstructed result of coaxial example, (a) actual density profile, (b) reconstruction using Compton scattering data only, (c) reconstruction using Tx data only, (d) reconstruction using both Compton and Tx data.
Ground truth and reconstruction results (again for the 4-source configuration) are shown in Fig.9.The results are generally similar to the three-target phantom, with the combination of Tx and Compton-scattered data showing best results.While all reconstructions underestimate density, Fig.9(d)shows more accurate values, particularly for the inner Delrin cylinder.In addition, the geometry of the outer cylinder is better defined, and artifacts outside the HDPE object are noticeably reduced.

FIGURE 10 .TABLE 2 .
FIGURE 10.Alternate reconstruction, which emphasizes the role of Compton data during reconstruction by setting γ tx = 0.5 and γ sc = 1.(a) reconstructed density (b) percent norm error vs. nominal density for all methods.While the alternate reconstruction does not improve SSIM or overall error, it slightly improves homogeneity in outer ring of the object.

Fig. 10 (
Fig.10(a)shows the result obtained by taking smaller steps based on transmission data, setting λ tx = 0.5 and λ sc = 1.Careful visual inspection suggests some improvement in the image, particularly in the homogeneity of the outer ring.To quantify this, the ground truth image in Fig.9(a) was used to identify all pixels corresponding to the outer ring.Variability of reconstructed density for these pixels was computed using percentile measures, with results shown in Table2.There appears to be a modest but measurable improvement in homogeneity, with the most homogeneous reconstruction given by the result plotted in Fig.10(a).Fig.10(b) compares the overall image error for various solutions as a function of iteration number, showing that varying the ratio λ tx /λ sc changes the shape of the error curve vs. ground truth values.While in this instance the lowest overall error is given by equal weighting, determining the optimum weighting for transmission vs. scattered data in the inversion is a potentially interesting avenue for future work.

FIGURE 11 .
FIGURE 11.Reconstructed result of clutter example, (a) actual density profile, (b) reconstruction using Compton scattering data only, (c) reconstruction using Tx data only, (d) reconstruction using both Compton and Tx data.

×
exp −a T l µ(E S k ) − b T l µ(E p )(26)The vector δg M can be expressed as δg M

TABLE 1 .
Quantitative comparison of methods.First number is SSIM, second number is normalized error (both computed vs. nominal density).High SSIM and low errors are desired; best result for each case is bolded.