# Compressed Sensing Inspired Image Reconstruction from Overlapped Projections

Authors: Lin Yang; Yang Lu; Ge Wang
Journal: International Journal of Biomedical Imaging
 (2010)

---

## Abstract

The key idea discussed in this paper is to reconstruct an image from overlapped projections so that the data acquisition process can be shortened while the image quality remains essentially uncompromised. To perform image reconstruction from overlapped projections, the conventional reconstruction approach (e.g., filtered backprojection (FBP) algorithms) cannot be directly used because of two problems. First, overlapped projections represent an imaging system in terms of summed exponentials, which cannot be transformed into a linear form. Second, the overlapped measurement carries less information than the traditional line integrals. To meet these challenges, we propose a compressive sensing-(CS-) based iterative algorithm for reconstruction from overlapped data. This algorithm starts with a good initial guess, relies on adaptive linearization, and minimizes the total variation (TV). Then, we demonstrated the feasibility of this algorithm in numerical tests.


---

## Body

## 1. Introduction

The popular CT scheme takes projection data from an X-ray source being scanned along a trajectory and reconstructs an image from these data that are essentially line integrals through an object. In real-world applications, higher temporal resolution has been constantly pursued, such as for dynamic medical CT, micro-, and nano-CT.  The multisource scanning mode is well known to improve temporal resolution but the data acquisition and field of view are seriously restricted to avoid overlapped projections, such as in the case of the classic dynamic spatial reconstructor (DSR). As shown in Figure1, here we consider reconstructing an image from overlapped projections so that a new dimension of freedom can be offered to design novel CT architectures.Figure 1
Imaging geometry for collection of overlapped projections from twoX-ray sources.In the overlapped projection geometry, two (or more) sources, A and B, emit X-rays simultaneously through an object to be reconstructed from various orientations. As a result, the resultant X-ray projections are overlapped onto the same detector array. The overlapped projections use the same detector array at the same time but complicate the imaging model. To perform image reconstruction from overlapped projections, the conventional reconstruction approach (e.g., filtered backprojection (FBP) algorithms) cannot be directly used because of the two problems. First, overlapped projections represent an imaging system in terms of summed exponentials, which cannot be transformed into a linear form, since the X-ray intensity through an object follows an exponential decaying function. Second, overlapped measurement carries less information than the traditional line integrals, due to the additional uncertainty from mixing two ray sums, leading to an underdetermined imaging system.Compressive sensing (CS) is a new technique being rapidly developed over the past years [1, 2]. It has been shown that if a vector x contains at most S nonzero elements and there are K random measurements of x such that K≥C·S·log(N), where C is a constant and N is the dimension of x, then minimizing the L-1 norm of x reconstructs x perfectly with an overwhelming probability. Inspired by its success in signal recovery, we propose a compressive sensing- (CS-) inspired iterative algorithm for reconstruction from overlapped data. This algorithm starts with a good initial guess, relies on adaptive linearization, and minimizes the total variation (TV). Then, we demonstrated the feasibility of this algorithm in numerical tests. The rest of this paper is organized as follows: in Section 2 we describe our algorithms for data synthesis and image reconstruction, in Section 3 we report numerical results under different conditions, and in Section 4 we discuss relevant issues and conclude the paper.


## 2. Methodology

An imagef can be discretized into a W by H matrix, which can be represented as a vector f of length n=W·H. Let Nsrc denote the total number of X-ray sources (for a dual-source system, Nsrc=2, which is the focused case in this paper), Nbin the total number of linear or area detector bins, and Nrot the total number of view angles. Then, the sampling process will yield Nbin·Nrot overlapped data. Since the X-ray attenuation is governed by an exponential decaying function, the overlapped projection data can be expressed as(1)p=[p1,1p2,1⋮pNbin,1⋮p1,Nrot⋮pNbin,Nrot]=exp(-M1f)+exp(-M2f)+⋯+exp(-MNsrcf)=[exp(-M1,1,1f)exp(-M1,2,1f)⋮exp(-M1,Nbin,1f)exp(-M1,1,2f)⋮exp(-M1,Nbin,2f)⋮exp(-M1,1,Nrotf)⋮exp(-M1,Nbin,Nrotf)]+[exp(-M2,1,1f)exp(-M2,2,1f)⋮exp(-M2,Nbin,1f)exp(-M2,1,2f)⋮exp(-M2,Nbin,2f)⋮exp(-M2,1,Nrotf)⋮exp(-M2,Nbin,Nrotf)]+⋯+[exp(-MNsrc,1,1f)exp(-MNsrc,2,1f)⋮exp(-MNsrc,Nbin,1f)exp(-MNsrc,1,2f)⋮exp(-MNsrc,Nbin,2f)⋮exp(-MNsrc,1,Nrotf)⋮exp(-MNsrc,Nbin,Nrotf)],
where p is an Nbin·Nrot by 1 vector whose element pm,r, m∈{1,2,…,Nbin} and r∈{1,2,…,Nrot}, is the overlapped projection datum detected by the mth detector bin at the rth view angle, Ml denotes the system matrix for the lth source, and l∈{1,2,…,Nsrc}. The 1 by n row vector Ml,m,r is the X-ray intersection length vector from the lth source to the mth detector bin at the rth view angle. The kth entry of Ml,m,r is obtained by calculating the intersection length of the involved X-ray through the kth pixel of f, which corresponds to the indices l, m, and r. The system matrix M and overlapped projection data can be readily computed in different ways. For example, in reference to [3], we have Algorithm 1 to generate projection data.Algorithm 1: Synthesis of overlapped projection data.
Initialization:(i)Initialize the coordinates for the two sources and detector bins:s1,s2,d;(ii)Define the rotation matrixQ=I2×2;(iii)Zero-out the data vectorpof sizeNbin·Nrotby 1;(iv)Specify the rotation angle incrementθ;Synthesis:Loop fork=1:Nrot:(i)With the current coordinates for the first source and the detector arraycompute the system matrixM1,k;(ii)With the current coordinates for the second source and the detector array;compute the system matrixM2,k;(iii)Update thekthblock of the data vectorp:pk=[p1,kp2,k⋮pNbin,k]=exp(-M1,kf)+exp(-M2,kf);(iv)UpdateQ:Q=[cos(θ)-sin(θ)sin(θ)cos(θ)];(v)Update source and detector coordinates:s1=Qs1;s2=Qs2;d=Qd.EndNow, the key issue is how to reconstruct an image from overlapped projection data. To alleviate the underdetermined measurement due to the overlapped nature of projection data, the compressed sensing (CS) principles are employed in our reconstruction process. To utilize the sparsity of an underlying image, it is first transformed into a gradient counterpart, and then the L-1 norm of the gradient, which is known as the total variation (TV), is minimized, subject to the overlapped projection data. The entire reconstruction process can therefore be casted into a constrained nonlinear optimization problem:Minimize: TV offsubject top=exp(M1f)+exp(M2f)and other constraints (such as intensity ranges and object features)Clearly, there are various ways to solve the above constrained TV minimization problem. In the CS field, a projection onto convex sets (POCS) and gradient descent search approach has been successfully used to solve this type of MRI and CT imaging problems [4–9]. POCS takes advantage of the fact that the linear constraints are hyper-planes in the n-dimensional space so that a closed form solution for the projection onto these hyper-planes can be derived. Such an algorithm works well in the single source geometry, because raw projection data can be processed into line integrals.However, in the case of overlapped projections from two sources, the constraint equations, which are the sums of two exponentials, cannot be transformed into a linear form. Therefore, a different approach is needed. Our solution is to make a good initial guess, such as a low-resolution CT image first. This blurry image will serve as a starting point, and the difference between this initial reference and the actual image will be iteratively updated, and at the same time the current guess will be also updated. Since the difference is assumed to be small, we can perform a Taylor series expansion to linearize the imaging system by omitting high-order terms. Then, we can apply the POCS-gradient algorithm on this linearly approximated system iteratively.Mathematically, let us denotef=g+df, where f is the original image, g the blurry image, and df the difference between f and g. Then, we have
(2)p=exp(-M1·f)+exp(-M2·f)=exp(-M1·g)exp(-M1·df)+exp(-M2·g)exp(-M2·df)=exp(-M1·g)(∑n=0∞(-M1·df)nn!)+exp(-M2·g)(∑n=0∞(-M2·df)nn!)≈exp(-M1·g)(1+(-M1·df))+exp(-M2·g)(1+(-M2·df))=exp(-M1·g)-exp(-M1·g)M1·df+exp(-M2·g)-exp(-M2·g)M2·df.
That is,(3)[exp(-M1·g)M1+exp(-M2·g)M2]·df=-p+exp(-M1·g)+exp(-M2·g).
Then, we have the approximate system(4)Mnewdf=pnew,
where(5)Mnew=exp(-M1·g)M1+exp(-M2·g)M2,pnew=-p+exp(-M1·g)+exp(-M2·g).
The above approximate system is linear with respect to df. This linearity allows us to perform POCS on df. To perform the gradient descent search on the TV of g+df, we compute the gradient of the TV explicitly in the image domain, for example, using the formulas described in [8]. After the linearization with respect to df and the formulation of the TV gradient, we can apply the POCS-gradient algorithm to estimate df. Note that such a reconstructed image g+df can be used as a new guess in the POCS-gradient process until a satisfactory reconstruction is achieved, as summarized in Algorithm 2.Algorithm 2:Image Reconstruction from overlapped projections (IROP).
Initialization:(i)  df=0;(ii)  g=A good initial guess, such as a blurry image from a low-resolution CT scan;(iii)  Nsrc=2;(iv)  Ndata=Nsrc·Nbin·Nrot;(v)  Nlinear=A number of the system linearization steps;(vi)Nitr=A number of the POCS-gradient or SART-gradient iterations for eachlinear approximation;(vii)   Ngrad=A number of the gradient descent search steps;(viii)  a=A step size for the gradient descent search, say 0.2;(ix) Obtain projection data using Algorithm1 or from a real scan;Reconstruction:Linearization updating loop forI=1:Nlinear(or until∥p-exp(M1g)-exp(M2g)∥pq≤ε)(i)ComputeMnew;POCS-gradient (or SART-gradient) reconstruction loop forJ=1:NitrForPOCSitr=1:Ndata(a)Projectdfonto the constraint equation:df=df-Mnewpnew-Mnewdf∥Mnew∥2;EndForGRADitr=1:Ngrad(b)Compute the TV with respect todf;(c) Perform the gradient descent search;EndEnd(ii)Update the current guessg=g+df;End


## 3. Numerical Experiments

To demonstrate the feasibility of our proposed algorithm for image reconstruction from overlapped projections, we developed a program in MATLAB, and implemented the traditional algebraic reconstruction technique (ART) for comparison. A Modified 2D Shepp-Logan phantom (Table2) was scaled into a 5 cm by 5 cm square and discretized into a 256 × 256 matrix. The phantom was centered at the origin of the reconstruction coordinate system.  A circular scanning trajectory of radius 121.66 cm was assumed with the two sources initially located at (-20 cm, 120 cm) and (20 cm, 120 cm), respectively. A 14 cm linear detector array was positioned opposite to the sources and 5 cm below the phantom with a distance of 131.53 cm from each of the sources. Gaussian white noise was drawn from the normal distribution N(0,0.005) and added to ideal projection data during the sampling stage. The scanning geometry is illustrated in Figure 1. The other parameters are listed in Table 1.Table 1
Parameters used in the numerical tests.
Test 1Test 2Test 3ART: Single SourceIROP: Two SourceART: Single SourceIROP: Two SourceART: Single SourceIROP: Two SourceNrot1515030150NA150Nbin50500100500NA500NlinearNA3NA3NA3Nitr50505050NA50NgradNA5NA5NA5Table 2
Parameters of the 2D modified Shepp-Logan phantom.
Axis length (a,b)Center (x,y)Angle (θ)Density(0.690,0.92)(0,0)01.0(0.6624,0.874)(0,-0.0184)0-0.8(0.11,0.31)(0.22,0)-18-0.2(0.16,0.41)(-0.22,0)18-0.2(0.21,0.25)(0,0.35)00.1(0.046,0.046)(0,0.1)00.1(0.046,0.046)(0,-0.01)00.1(0.046,0.023)(-0.08,-0.605)00.1(0.023,0.023)(0,-0.606)00.1(0.023,0.046)(0.06,-0.606)00.1We performed both ART and IROP reconstructions under these conditions, with blurry and constant initial guesses. Representative results are in Figures2, 3, and 4. It has been observed in our simulation that our IROP algorithm would work well if the initial guess resembles the ideal image through a moderate blurring process. Actually, in the first test the blurry images were obtained by blurring a low-quality ART image reconstructed under a severely under-sampling condition with only Nbin×Nrot = 15×50 = 750 measurements to reconstruct 256 × 256 = 65536 pixels. In the second test, more measurements were made in the single source scan, and the IROP reconstruction became better. Also, the IROP reconstruction turned to be smoother than the corresponding ART images, indicating that compressed sensing (CS) is more effective than ART in suppressing image noise. In the last test, we reconstructed an IROP image with a constant initial guess (a zero image).  The reconstructed image can be further improved if we use more iterations.Reconstructed images of the Shepp-Logan phantom in the first test. (a) A reconstruction using ART, (b) a reconstruction using IROP, and (c) the profiles along the central vertical line of the phantom, where the dotted and solid lines are for the phantom and the IROP reconstruction, respectively (the display window: [0,0.5]).
(a)(b)(c)Reconstructed images of the Shepp-Logan phantom in the second test. (a) A reconstruction using ART, (b) a reconstruction using IROP, and (c) the profiles along the central vertical line of the phantom, where the dotted and solid lines are for the phantom and the IROP reconstruction, respectively (the display window: [0,0.5]).
(a)(b)(c)Reconstructed images of the Shepp-Logan phantom in the third test. (a) A reconstruction using IROP with a constant initial guess, and (b) the profiles along the central vertical line of the phantom, where the dotted and solid lines are for the phantom and the IROP reconstruction, respectively (the display window: [0,0.5]).
(a)(b)To investigate the convergence of IROP, we first introduce an evaluation metricδ(n), which is defined as the sum of the component values in the error vector df at nth iteration:(6)δ(n)=∑i=12562dfi(n),
where the subscript i denotes the ith pixel component in the error vector df. We then plotted δ(n) for every iteration. The results for each of the three tests are shown in Figure 5. There are mainly three important observations from the convergence plots. First, the big jumps at multiples of Nitr indicate that linearization step played an important role in the convergence of IROP. Second, as IROP goes through more iterations, δ(n) approaches zero, showing that IROP effectively reduces the differences between the original and the reconstructed images.  Finally, the smaller values for δ(n) in test 2 show that a good initial guess can lead to better reconstruction quality. Figure 6 shows another 3 plots obtained with fewer iterations and more linearizations. It is observed that the error between the reconstructed image and the true image was reduced dramatically immediately after each new linearization. After 3 or 4 linearization processes, the convergence curve became stable and smooth. After that, if we increased the number of iterative steps in each linearization process, the image quality would be improved only slowly. Hence, to balance image quality and computational time, a good solution is for the method to use a limited number of iterative steps after each earlier linearization process and perform a sufficiently large number of iterative steps after the final linearization, for example, after 3 or 4  linearization processes.Convergence plots for (a) Test 1, (b) Test 2, and (c) Test 3, with 3 linearization steps and 50 iterations after each linearization.
(a)(b)(c)Convergence plots for (a) Test 1, (b) Test 2, and (c) Test 3, with 10 linearization steps and 3 iterations after each linearization step.
(a)(b)(c)


## 4. Discussions and Conclusion

The primary advantage of the IROP scheme is to improve the data acquisition speed. In one exemplary application, we can assume that the two sources are fairly close so that the detector collimation can work effectively for both the sources. If a good number of sources are used, scattering effects could be a concern. In that case, scattering correction may be needed using hardware (such as some degree of multiplexing) and/or software (such as model- or image-based compensation) methods [10–13].In Algorithm2, the key for the linearization to be successful is to have a good initial guess. It is underlined that it is practical to have such a good guess.  For example, in multiresolution CT studies, a low-resolution image serves as a guess naturally.  Also, in dynamic CT studies, an initial image represents a good guess to subsequent images. When we have a cluster of computers, we may use multiple random initial guesses to search for a more accurate and stable reconstruction.  Furthermore, Algorithm 2 may be adapted into an evolutionary scheme.The implementation of Algorithm2 can be improved in several ways. To reduce the smoothing artifact, one can reduce the number of gradient descent iterations or the step size. Other algorithmic parameters could also be tuned for a specific type of applications. Most importantly, the computational structure of Algorithm 2 is really based on simple heuristics and does not reflect all the constraints and requirements in a well-integrated and optimized fashion. It is possible and desirable to design brand new algorithms that involve less parameters and have better properties.A theoretical analysis on the convergence of the IROP scheme has not been performed yet but we hypothesize that the global convergence can be established if a guess is appropriately chosen, as numerically shown in the preceding section. Actually, the IROP problem is much better posed than many well-known inverse problems such as diffuse optical tomography (DOT) [14].  In IROP with two sources, each datum reflects information from two lines.  In DOT, each measure is related to a random zigzag trajectory.  Thus, it is not surprising to see better results with IROP than that with DOT. When we have infinitely many sources along a line, we have a line-source imaging geometry, which has been studied by Bharkhada et al. [15] and still yields better results than DOT reconstruction [15].Since the IROP scheme mixes line integrals pairwise, the IROP problem may lead to an underdetermined system of measurement equations, especially when the number of samples is not sufficiently large for ultrafast imaging performance.  To address this issue, we have implemented the CS principles in Algorithm2 by minimizing the TV. CS is a contemporary technique for solving an underdetermined system of linear equations, whose solution is known to be sparse. The main idea is to minimize cardinality, or equivalently to minimize the TV in many cases. In the context of IROP, an image itself is usually not sparse, but it can be sparsified in a transformed domain such as the gradient transform, and then we can apply the L-1 norm minimization in the transformed domain subject to the projection data constraints for good reconstructions, as numerically shown in the preceding section.It is emphasized that our IROP approach can be extended to multiple other imaging scenarios.  For example, in transmission ultrasound imaging, we can use multiple ultrasound sources and a single array of detectors (transducers). This may be also related to the area of signal unmixing. The common task would be to unravel an underlying signal or image from mixed measures. There seem good research opportunities along this direction.In conclusion, we have proposed the idea to perform image reconstruction from overlapped projection data and formulated a CS-based iterative algorithm for this new imaging problem. Our IROP algorithm starts with a good initial guess, relies on adaptive linearization, and minimizes the TV. Also, we have demonstrated the feasibility of this algorithm in numerical simulation. Further research is being performed to characterize and improve our IROP approach.



---

*Source: 284073-2010-06-22.xml*