|Home | About | Journals | Submit | Contact Us | Français|
To maximize the time-integrated X-ray flux from multiple X-ray sources and shorten the data acquisition process, a promising way is to allow overlapped projections from multiple sources being simultaneously on without involving the source multiplexing technology. The most challenging task in this configuration is to perform image reconstruction effectively and efficiently from overlapped projections. Inspired by the single-source simultaneous algebraic reconstruction technique (SART), we hereby develop a multisource SART-type reconstruction algorithm regularized by a sparsity-oriented constraint in the soft-threshold filtering framework to reconstruct images from overlapped projections. Our numerical simulation results verify the correctness of the proposed algorithm and demonstrate the advantage of image reconstruction from overlapped projections.
Since the first computed tomography (CT) scanner was made , all the commercial scanners have been employing the X-ray source with a single small focal spot, which can be modeled as a point source. In micro-CT and even nano-CT applications, the reduced X-ray focal spot size becomes increasingly a limiting factor to achieve contrast and temporal resolution targets. To address this issue, our group recently proposed to use a line-shaped X-ray source so that more photons can be generated in a given data acquisition interval . In this context, the X-ray source can be modeled as a line-segment, which can be further discretely modeled as an array of points . In single source CT scanners, contrast resolution is limited by the finite focal-spot size necessary to generate a sufficient number of X-ray photons, and temporal resolution is limited by the time taken to acquire sufficiently many projections over a full-scan or half-scan angular range. Since a line source covers a wide angular range per view, the number of photons is increased to radiate an object to be reconstructed. Therefore, use of a line-shaped X-ray source or a multiple source array is a candidate scheme to balance spatial, contrast and temporal resolution.
Interestingly, the technology of field-emission X-ray sources based on carbon nanotubes (CNT) is a recent invention with several intrinsic advantages over conventional X-ray tubes. To maximize the time-integrated X-ray flux from multiple sources and improve the signal-to-noise ratio (SNR), multiplexing was used, where multi-source are excited with different temporal modulations . If two or more sources are simultaneously fired, X-ray photons reach the same detector together for any single measurement, and one would not be able to identify which photon comes from which source. To unmix the signals from various sources, several projections can be collected at different time instants for the same view . Due to the limited readout speed and mixed signals, there seems little advantage to use multiplexing configuration in terms of contrast and temporal resolution. In collaboration with Dr. Otto Zhou's group at University of North Carolina, we are developing a novel multi-source micro-CT system. In our system, we plan to fire a multi-source system simultaneously and acquire multiple projections at any viewing angle on nonoverlapped segments of a shared detector array . With recent developed compressive sampling (CS) techniques [7, 8], we are also working to improve temporal resolution and reduce radiation dose from a limited number of nonoverlapped projections .
Here, we consider how to reconstruct an image from overlapped projections. Previously, our group developed a generalized simultaneous algebraic reconstruction technique (SART) algorithm to reconstruct an image from data collected with an X-ray line source . Assuming that the differences between measured and predicted projections from various source points have equal weights, these differences are then backprojected to different X-ray source points in the SART framework. However, this algorithm converges to the least square solution that is not necessarily the true image. In our simulation, the reconstructed images sometimes suffered from blurring . Then, our group developed a CS-based algorithm to solve this problem . The algorithm is implemented in a projection onto convex sets (POCS) framework and employed a steepest gradient searching strategy. Although this algorithm often converges to the true image, its convergence speed is rather slow. Because the SART framework has an excellent convergence behavior especially when the ordered subset (OS) format is applied, in this paper we will develop a SART-type algorithm for image reconstruction from overlapped projections.
The paper is organized as follows. In the next section, we will formulate the overlapped projection model, which is not a sum of line integrals and quite different from that in the single-source case. In the third section, we will design a SART-type reconstruction algorithm for image reconstruction from overlapped projections. In the fourth section, we will report numerical results. In the last section, we will discuss related issues and conclude the paper.
For CT reconstruction, a two-dimensional digital image can be expressed as f = (fi,j) I × J, where the index 1 ≤ i ≤ I and 1 ≤ j ≤ J are integers. Define
with 1 ≤ n ≤ N and N = I × J, we have the image in a vector representation f = [f1,f2,...,fN]T N. In this paper, we will use both the signs fi,j and fn for convenience. Let g = [g1,g2,...,gM]T M be a measured vector with M being the product of the number of projections and the number of detector elements. They are linked by the following linear system:
where A = (am,n) M × N is a measurement matrix. In a typical fan-beam geometry, the nth pixel can be viewed as a rectangular region with a constant value fn, the mth measured datum gm as an integral of partially covered pixel areas by a narrow beam from an X-ray source to a detector element which are weighted by the corresponding X-ray linear attenuation coefficients. Thus, the component am,n in (2) denotes the intersection area between the nth pixel and the mth fan-beam ray (Figure 1). While the whole matrix A represents the forward projection, AT implements the backprojection.
While the imaging model (2) is valid for a single-source system, it cannot be used for multi-source-generated projections. The reason is that the measured data in (2) has been postprocessed by a logarithmic operation. In other words, we must model the raw data directly. For a multi-source system with Q sources, assuming that the qth source emits Iq photons towards each detector element. According to Beer's law, we have the following imaging model:
where Aq is the system matrix defined in Section 2.1 for the qth source, s a constant to normalize the area model (Figure 1) for the measurement matrix to the line integral model, and e−s Aqf is a vector whose element is the exponential function of the corresponding element of −s Aqf. In practice, we can approximate s as the reciprocal of the average width of an X-ray path in an object to be reconstructed. If we assume that all the X-ray sources have the same intensities, I1 = I2 = = Iq = I0, we have
where the constant s is absorbed by , and 1/I0 absorbed by . The key task in this paper is to reconstruct from and known measurement matrices Aq.
In the past decade, our group studied a block-iterative (BI) or ordered-subset (OS) version of a general Landweber scheme , of which the SART and OS-SART  are special cases, for minimization of a weighted least square functional in the real/complex space, and proved its convergence under quite general conditions [12, 13]. The SART or OS-SART technology has been widely used in the CT field. Particularly, for the system (2) the SART-type solution is given by 
where , , Am is the mth row of A, k the iteration index, and 0 < λk < 2 a free relaxation parameter. Let Λ+N N × N be a diagonal matrix with Λn,n+N = 1/a+n and ΛM+ M × M be a diagonal matrix with Λm,mM+ = 1/am+, (5) can be rewritten as
Note that for all the iteration index k, ΛM+, Λ+N and A remain unchanged.
Since (4) is a non-linear equation, there is no analytic solution to this problem. Here, we will find a solution in the SART framework summarized in the proceeding subsection. Our strategy is to linearize the equation and approximate the solution successively. Denote the approximate solution for (4) as after k iterations, we have
where is the error image, and “” represents the component-wise multiplication of vectors of the same size. Then, in (7) can be expanded in a Taylor series
where Substituting the 1st order approximation of into (7), we have
where Eq,k M × M be a diagonal matrix with . Equation (9) can be rewritten as
Since the projection errors from the involved sources are known, (10) can be understood as a linear system of with a measurement matrix Bk.
and ΛM+,k M × Mbe a diagonal matrix with
we have the SART-type solution for
where l is the iteration index and λk,l the relax parameter. Since
Equation (14) can be rewritten as
Once we have an approximation solution for after L iterations,we can update the reconstructed image by
For example, we can choose as the initial image and set L = 1 for one step iteration to approximate , which results in a simplified algorithm
For numerical implementation, our SART-type reconstruction algorithm for image reconstruction from overlapped projections can be summarized in the following pseudocode:
In the above pseudo-code, (S.1.) initializes the iteration index k, the relax parameter λ, and the initial image . In our numerical simulation, we always set λk = 1 and . The outer loop (S.2)–(S.7) solves for successively. (S.2) precomputes several important intermediate variables to update an reconstructed image in the iteration step k. (S.3)-(S.4) computes the current error image for one step iteration according to (16). Because the backprojection operation for different X-ray sources in the inner loop (S.4) has the same structure, it can be implemented by calling a common procedure. (S.5) updates an reconstructed image. (S.6) updates the iteration index. In (S.7), the convergence criteria are checked. The stopping criteria for (S.7) can be the maximum iteration number is reached and/or the relative reconstruction error comes under a predefined threshold . Finally, the reconstructed image f* is scaled by dividing with the constant s.
The conventional data acquisition is based on the Nyquist sampling theory, which states that for accurate reconstruction of a band-limited signal or image the sampling rate must be at least twice the highest frequency of the signal or image. However, the recently developed CS theory shows that a high-quality signal or image can be reconstructed from far fewer measurements than what is usually required by the Nyquist sampling theorem [7, 8]. In light of the work on solving the linear inverse problems with sparsity constraints by Daubechies et al. [14, 15], we recently adapted the single source SART to reconstruct an image from a limited number of projections subject to a sparsity constraint , and demonstrated that the sparsity constraint helped improve the quality of reconstructed images effectively and reduce the number of projections significantly. Using the same strategy described in our previous papers [16, 17], here we use the sparisty constraint to regularize the proposed multi-source SART algorithm. This can be done by adding a soft-threshold filtering step between (S.5) and (S.6) in the pseudo-code given in Section 3.2. Particularly, we have the following pseudo-code segment:
In the above pseudo-code, the sparse transform in (SS.5.1) can be either any invertible lossless compressible transform such as wavelet transform  and Fourier transform or uninvertible transforms such as discrete gradient transform (DGT) and discrete difference transform (DDT) . For an uninvertible transform, the inverse sparse transform in (SS.5.4.) is in terms of pseudoinversion as we performed for DGT and DDT . (SS.5.2.) determines an optimal threshold automatically using the projected gradient method for fast convergence . In fact, we can omit (SS.5.2.) and specify any fixed filtering threshold. However, both the convergence speed and final result depend on the choice of the threshold.
To verify the proposed SART-type algorithm for image reconstruction from overlapped projections, we implemented it in Matlab on a PC, with the computationally intensive segments coded in C and linked via the MEX mechanism. As illustrated in Figure 2, we simulated a triple-source fan-beam micro-CT system. In the system, the source XS0 is rotated on a circular scanning locus of radius 120mm. The object was a modified Shepp-Logan phantom in a compact support with a radius of 35mm. We used an equi-distance detector array of length 120mm. The detector was perpendicular to the direction from the origin to the X-ray source XS0 that is 40mm from the system origin. The detector array consisted of 500 elements. On the line through the X-ray source XS0 and parallel to the detector, we put two sources XS1 and XS2 being 25mm apart from XS0 on its right and left sides, respectively. With this triple-source configuration, we simulated single-source (only XS0 was fired), dual-source (XS1 and XS2 fired simultaneously), and triple-source cases (XS0, XS1 and XS2 fired simultaneously).
For each of our selected numbers of projections over a full-scan range, we first equiangularly acquired the corresponding projection dataset based on the aforementioned projection model in the single, dual, and triple source cases, respectively. Then, we reconstructed the images using our algorithm described in Section 3.2. In our simulation, the parameter λk in the SART formula (18) was set to 1.0, and the stopping criterion was defined as reaching the maximum iteration number 5000. Because the Shepp-Logan phantom is a piecewise constant function, its DGT and DDT are sparse. Hence, we also employed the sparsity regularization in terms of total difference minimization  and the threshold for filtering was automatically computed using the projected gradient method . Figure 3 shows the reconstructed 256 × 256 images from 9, 11, 13, and 15 projections, respectively. For real-world applications, measurement noise is unavoidable. To test the stability of the proposed algorithm against data noise, we repeated the above reconstructions from projections corrupted by Poisson noise, assuming I0 = 5 × 104 photons per detector element . The results are in Figure 4, which indicate the stability of the proposed algorithm.
In the CT field, the line integral model along an X-ray path has been widely used in consistency with Beer's law. However, it does not reflect the divergence due to the combination of the finite detector size and the source focal spot. As shown in Figure 1, we have assumed an area model for the X-ray path and normalized it for the multi-source imaging model. Our area model treats the X-ray path as a narrow fan-beam from the X-ray point source to the detector element, and we believe that the area model works better than the line model. Note that the proposed algorithm views the projection procedure as a matrix transform, both the area and line models can be handled by our algorithm. In other words, the proposed algorithm is independent of the imaging model as long as it is linear or can be transformed into a linear one. Additionally, to simplify the derivation and demonstrate our idea, we have assumed that the photon numbers emitting from all the sources to each detector are the same. In fact, these numbers may be different, and can be easily incorporated into our algorithm.
As far as the convergence of the proposed algorithm is concerned, it should converge to the least square solution in the cases of either noise-free and noisy projections. The reason is that the proposed method is in the framework of the general Landweber scheme, whose convergence has been well studied under quite general conditions [12, 13]. When only a small number of projections are available, we can use some sparse constraints to steer the solution to the truth. However, the convergence speed of the current soft-threshold filtering technology is still slow although it has been accelerated using the projected gradient method . In the future, we will employ more advanced techniques for a faster speed, which may include but not limited to optimizing the code, employing parallel computation, and developing new algorithms.
In conclusion, we have developed a SART-type algorithm for image reconstruction from overlapped projections. The algorithm has been verified and demonstrated in the numerical simulation. Our methodology has a potential to support more flexible designs of multi-source CT/micro-CT systems for better contrast and temporal resolution.
This work is partially supported by the NSF/MRI Program (CMMI-0923297) and NIH/NIBIB Grants (EB009275 and EB011785). The authors are grateful for constructive comments from the reviewers.