|Home | About | Journals | Submit | Contact Us | Français|
Single photon emission computed tomography (SPECT) is an important biomedical imaging modality. However, since gamma cameras are expensive and bulky, truncated projection data are either preferred or unavoidable. Inspired by the recent results on interior tomography in the x-ray CT field, here we present the interior SPECT approach for exact and stable reconstruction of a region of interest (ROI) from uniformly attenuated local projection data, aided by prior knowledge of a sub-region in the ROI. First, by analytic continuation we prove that interior SPECT is both exact and stable, and by singular value decomposition (SVD) we establish the stability of interior SPECT. Then, given the constant attenuation coefficient and object boundary, our interior SPECT reconstruction is achieved by inverting a generalized truncated Hilbert transform using the SVD technique. Preliminary numerical simulation data demonstrate that our work has practical utilities. The theoretical generalization of our work to the variable attenuation case is underway, and the same numerical approach can be applied.
Facing the increasing radiation risk for CT examinations, a number of image reconstruction algorithms were developed to reduce the amount of necessary raw data. A recent milestone is the two-step Hilbert transform method developed by Noo et al. . In their framework, an object image on a PI-line/chord can be exactly reconstructed if the intersection between the chord and the object is completely covered by the field of view (FOV). In 2006, Defrise et al.  proposed an enhanced data completeness condition that the image on a chord in the FOV can be exactly reconstructed if one end of the chord segment in the object is covered by the FOV. Inspired by the tremendous biomedical implications including local cardiac CT at the minimum dose, local dental CT with high accuracy, CT guided procedures, nano-CT and so on , by analytic continuation we proved that the interior problem (exact reconstruction of an interior ROI only from data associated with lines through the ROI) can be exactly solved if a sub-region in an ROI/VOI in the FOV is known [4–6] from fan-beam/cone-beam projection datasets, while the conventional wisdom is that the interior problem does not have a unique solution . Similar results were also independently reported by others [8, 9].
While the CT reconstruction algorithms are being advanced rapidly, the single-photon emission computed tomography (SPECT) techniques are also experiencing remarkable improvements [10–14]. As a unique biomedical tomographic imaging technique, SPECT is to reconstruct an image of the radioactive source distribution. SPECT is performed with a gamma camera to acquire multiple 2D projections from multiple angles. Then, a tomographic reconstruction algorithm is applied to the fan-beam/cone-beam projections, yielding a 2D/3D image. Different from the line integral model for x-ray imaging, the SPECT projections can be mathematically modeled as a exponentially attenuated Radon transform [14, 15]. In this context, the CT reconstruction may be regarded as a special case of SPECT (all the attenuation coefficients are zeros). However, the reconstruction techniques of CT cannot be directly used for SPECT in general.
Inspired by the recent interior reconstruction results [4–6, 8, 9], here we report the corresponding results in terms of the attenuated Radon transform for SPECT. In the next section, we will summarize the major SPECT interior reconstruction results, assuming the constant attenuation coefficient and the object boundary are known. In the third section, we will analyze the uniqueness and stability. In the fourth section, using the generalized Hilbert transform and singular value decomposition (SVD) method we will describe a practical implementation of interior SPECT reconstruction, preliminary numerical analysis and simulation results. In the last section, we will discuss some relevant issues and conclude the paper.
Let f (x) be a 2D smooth distribution function defined on a convex compact support Ω with x=(x, y) Ω. Mathematically, in a parallel-beam geometry (Fig. 1), the SPECT projections of f (x) can be modeled as the attenuated Radon transform
where the subscript “o” indicates original projection data, θ = (cos θ, sin(θ)), θ = (−sin θ, cos(θ)), and μ(x) is the attenuation coefficient map on the whole compact support Ω. For practical applications Po (θ, s) can be acquired by a gamma camera. For most of the imaging subjects such as the brain, the attenuation map can be approximated as a uniform distribution
where μ0 is a constant. Since the object function is compactly supported, we can determine the maximum coordinate along the direction θ of the intersection between the support Ω and the integral line of the projection Po (θ, s). Without loss of generality, we denote this coordinate as tmax (θ, s) and Eq. (1) becomes
For practical applications, we can assume that the convex compact support Ω and the constant attenuation coefficient μ0 are known. Once the coordinate system is chosen, we can immediate determine tmax (θ, s) for all the whole projection dataset. By multiplying a weighting factor eμ0tmax (θ,s), the projection model of SPECT can be reduced to
where the subscript “w” indicates that the weighted projection data.
Denote a weighted backprojection of the differential projection data by
In 2004, Rullgard proved a relationship linking the object image f (x) and the weighted backprojection g(x) as 
where “PV” represents the Cauchy principle value integral, and chμ0 is defined by
When μ0 → 0, chμ0 becomes the Hilbert transform kernel, in consistence with the results in the CT field . Eq.(6) is called a generalized Hilbert transform, and the corresponding Hilbert filtering line a generalized PI-line. The backprojection operation for g(x) in Eq. (5) only involves local projection data whose integral paths intersect with the point x. Hence, inside a field of view (FOV) or region-of-interest (ROI) any point is irradiated at least from an angular range of 180 degrees, and its backprojection g(x) of differential data can be exactly computed.
Without loss of generality, let us denote a 2D f (x) on a PI-line/chord as f (u) and the backprojection g(x) as g (u), where u is a 1D coordinate along the PI-line. Denote the intersection of the convex compact support Ω and the PI-line as the interval (csb, cse) with csb < cse, which implies that f (u) = 0 for u (csb, cse) where the subscripts “sb” and “se” represent the starting and ending points of the support. Also, let us denote the intersection of the FOV and PI-line as an interval (cvb, cve) with cvb < cve where the subscripts “vb” and “ve” represent starting and ending points of the FOV. Using the above notations, Eq. (6) can be reduced as
In 2007, Noo et al. extended their two-step Hilbert CT method to the case of SPECT , and showed that an object image f (u) can be exactly reconstructed on the whole compact support interval (csb, cse) if cvb < csb < cse < cve (Fig. 2(a)). In this paper, we will further extend this result for SPECTin the case csb < cvb < cve < cse. Suppose that there exists a real number cke satisfying csb < cvb < cke < cve < cse with f (u) being known on the interval (cvb, cke), in the following we will show that f (u) can be exactly and stably reconstructed on the interval [cke, cve) (Fig. 2(b)), where the subscript “ke” represents the ending point of the known part. Our main results can be summarized as:
Assume that csb < cvb < cke < cve < cse and a smooth function f (u) is supported on (csb, cse). The function f (u) can be exactly reconstructed on (cvb, cve) if (i) f (u) is known on (cvb, cke), (ii) g(u) is known on (cvb, cve), and (iii) the constant μ0 and mμ0 are known, where
Note that we have omitted a parameter ckb = cvb in the above theorem, and in practical applications mμ0 can be directly computed from the projection datum along the path through the corresponding PI-line.
Similar to what was discussed in the CT case , we can assume that the known f (u) are on any subinterval of the line-segment (cvb, cve) or a union of such intervals. The corresponding results can be directly obtained by applying Theorem 2.1 repeatedly.
Without loss of generality, we assume csb = −1 and cse = 1. Otherwise, we can arrive at this standardization by a linear coordinate transform. For the case cvb < −1 and 1 < cve, Noo et al. reduced the reconstruction issue to a Fredholm integral equation of the second kind for h(u) on the interval [−1,1] :
What is more important, Noo et al. proved that the Cauchy principle value integral in Eq. (14) can be removed because Eq. (15) is a smooth function, leading to that μ0(u, ũ) is a smooth and continuous function over the region (u, ũ) [−1,1]×[−1,1].
Using the following identical equation
μ0 (u, ũ) can be rewritten as
Note that both (u, v, ũ) and η(μ0(u−ũ)) are analytical with respect to the variable u, and μ0(u, ũ) is analytical with respect to the variable u. Hence, (Kh)(u) can be analytically extended to the whole complex plane .
Now, let us consider our case −1 = csb < cvb < cke < cve < cse = 1. First, we rewrite Eq. (10) as:
By our assumption in Theorem 2.1, h1(u) is known for u R and h(u) is known on the interval (cvb, cke), we obtain that
is known on the interval (cvb, cke). Because the principal integral of the first term of Eq. (21) is defined on (−1, cvb) (cve, 1), h2 (u) is analytic on the interval (cvb, cve), and it can be analytically extended to the complex plane with cuts along the real axis from −∞ to cvb and from cve to +∞. Therefore, h2(u) can be determined on the interval (cvb, cve) by its value on the interval (cvb, cke). As a result, h(u) (or f(u)) can be uniquely reconstructed on the interval (cvb, cve).
Using techniques similar to that in [2, 16], a stability analysis can be performed for our theorem 2.1 as follows. In practice, we have the measurement gr(u) of g(u) and reconstruct fr (u) from gr (u). Assume that these functions satisfy the following conditions:
The above conditions lead to
Similarly, we denote
Choose two arbitrary constants vb, ve such that cvb < vb < cke < ve < cve. On the interval (vb, ve), we have
Especially, on the interval (vb, cke), we have
Let Ω be a domain and D be a segment of the boundary Ω. If f (z) is an analytical function in Ω such that
and ω(z) is a harmonic function in Ω such that
Let B to be the disk with diameter (vb, ve), Ω = B\(vb, cke), then Ω= D + B, where D = (vb, cke), and B is the circle with diameter (vb, ve) (Figure 3). Recall that
For z B,
We conclude that
On other hand, there exits a harmonic function ω(z) satisfying 
and for u (cke, ve),
Therefore, by Nevanlinna’s principle, for u (cke, ve)
One remark is in order that for u (cke, ve) and 0 < ω(u) < 1, ω(u) → 1 as u → ve). Eq. (35) implies that the reconstruction of f(u) for u (cke, cve) is stable near cke and probably less stable near cve.
Let us rewrite Eq. (8) as
and discretize the u -axis with a uniform sampling interval δ satisfying sampling theorem for both f(u) and gw(u). Denote f(u) at the discrete sampling points on the interval (csb, cse) as f1, f2, … fn …, fN. Also, denote gw (u) at sampling points on the interval (cvb, cve) as g1, g2, … gm, … gM. Eq. (36) can be then discretized as
In Eqs. (38) and (39), “T” represents the transpose operator, and in Eq. (40) Am, n is the weighting coefficient of fn. Let u(fn) and u(gm) be the coordinates on theu -axis for fn and gm. Based on the discrete Hilbert Kernel formulated in , we can write Am, n as
With our sampling method, always produces integers, and the case in Eq. (42) exactly corresponds to the singular point in the Cauchy principal integral in Eq. (36). Because f(u) is known on (cvb, cke), we may divide F into two parts
where Fk is the known part on the sampling points in (cvb, cke), while Fu is the unknown part on the sampling points in (csb, cse)/(cvb, cke). Correspondingly, the matrix A can be divided into two parts
Then, we immediately arrive at a linear equation system
with = G − AkFk, Ā = Au and = Fu.
Although Eq. (45) is generally ill-posed, we can still stably reconstruct the part of in the interval (cke, cve) using the singular value decomposition (SVD) method [19, 20]. Note that the dimension of is = M. If the dimension of is , the dimension of Ā will be × . According to the SVD theory [19, 20], the matrix Ā has a SVD decomposition in the following form
where U and V are orthogonal matrices of × and × respectively, and Λ is an × diagonal matrix whose diagonal elements λq satisfying λ1≥ … ≥ λq ≥ … ≥ λQ, Q = min(, ). In this way, we have a stable numerical solution for [19, 20],
where Λ−1 is a diagonal matrix of × whose diagonal elements are defined as
After a rigid stability analysis has been performed in Section 3.2., here we will numerically analyze the stability of Theorem 2.1 from a signal processing viewpoint. Recall that the unknown f (u) and known g (u) are linked by Eq. (8) and we have assumed that f (u) is smooth. Based on the sampling theorem in the classical signal processing theory, f (u) can be exactly recovered from its values at sampling points if (i) f (u) is band-limited and (ii) the sampling frequency is not smaller than the Nyquist frequency. In addition to the smoothness of f (u), we can further assume that f (u) is essentially band-limited, which is reasonable in most practical engineering applications . Hence, we can stably reconstruct f (u) at its finite sampling points on the interval (cvb, cve) using the aforementioned SVD method.
Recall that = Ā = AuFu, we have
where Iu represents an × unit diagonal matrix. Let Eu = u − Fu as the error vector of Fu, and Su = VΛ−1UT Au − Iu as an × matrix. Then, Eq. (49) can be reduced as
For any 1 ≤ ≤ , the absolute reconstruction error can be expressed as
where represents the maximum absolute value of the elements in Fu, and . Because reflects the accuracy of the reconstructed u, we call it the precision measure in this paper.
Next, let us consider the case of noisy data. Let us assume that
where ε defines the maximum noise magnitude. Recall that = G − AkFk, we have the noise expression of as
Finally, based on Eq. (47), the noise in the reconstructed image can be expressed as
which can be further reduced as
with Pun = (VΛ−1UT −VΛ−1UTAk) and . Clearly, the dimension of un is × 1. Assume that the dimension of Wun is × 1, the dimension of Pun must be × . For any 1 ≤ ≤ , the absolute magnitude of noise in the image can be further expressed as:
where and is called the stability measure in this paper. By evaluating the values of , we can directly reveal the stability of our algorithm for interior reconstruction of SPECT.
To demonstrate the exactness and stability of interior reconstruction of SPECT, we set csb = −1, cse = 1, cvb = −0.6, cve = 0.6 and cke = −0.2 for numerical computation. Assume that the number of sampling points on [−1,1] is N and the small constant in the TSVD method is ε. For different combinations of N, ε and μ0, we computed the precision and stability measures and as shown in Figures 4 and and5,5, respectively. Because f (u) was known on (cvb, cke), the corresponding portions of and in the interval (cvb, cke) were set to 0. From Figures 4 and and5,5, we can make the following three comments. First, there are little differences in terms of the precision and stability measures with respect to the number of sampling points on the whole support of f (u). If we continuously increase the number N, the sampling interval δ will become smaller and smaller. Eventually, we will arrive at the continuous case when N → ∞. It is expected that there would still be little difference in terms of the precision and stability measures of N → ∞ from what we have obtained (which is of course a numerical observation lack of mathematical rigor). Second, the precision measure in Figure 4 is very close to zero in the interval [cke, cve), especially the portions near cke. When we reduced the constant ε in the TSVD method, the precision measure would decrease in the interval [cke, cve), which means a higher accuracy. Undoubtedly, should ideally approach zero on the interval [cke, cve) when ε → 0. Meanwhile, a small constant attenuation parameter μ0 yields a better precision. Third, the stability measure in Figure 5 was smaller on the interval [cke, cve) than elsewhere (especially the portion near cke), which means a better stability. If we continuously decrease the constant ε in the TSVD method, the stability measure would correspondingly increase in an order of 1/ε, which means the stability would become worse and worse.
To demonstrate the usefulness of the proposed interior reconstruction for SPECT, we performed several numerical tests. In our simulation, a modified Shepp-Logan phantom of SPECT version was employed inside a circular support Ω of radius 10cm. As shown in Figure 6(a), the phantom consists of 10 ellipses whose parameters are the same as that given by Noo et al. . We assumed that a uniform attenuation coefficient μ0 throughout the whole support of the phantom. In parallel-beam geometry, 360 view-angles were uniformly sampled for an angle range [0, π]. For each view-angle, there were totally 600 detector elements uniformly distributed along a 20 cm linear detector array. For different attenuation coefficients μ0 (0.0cm−1, 0.15cm−1, 0.30cm−1), we analytically computed non-truncated projection datasets. Hence, the backprojection of differential data can be easily calculated at any point to simulate different fields of view (FOVs).
Without loss of generality, here we assumed a rectangular FOV and the distribution function f (x) was known in the stripe region in the FOV (Figure 6(b)). We implemented the TSVD method described in Section 3.2 to reconstruct the distribution image f (x) in the whole FOV. The constant parameter ε was set to be 0.02 and there were 512 uniform sampling points along each 20 cm side length in the image space. In the reconstructed procedure, the prior information of f (x) in the strip region was chosen as that reconstructed from the non-truncated data with μ0 = 0. The reconstructed results are presented in Figure 7. Because there was noise in the prior information, we had some errors in the reconstructed images, especially in the case μ0 = 0.30 cm−1. These results are consistent with our aforementioned stability analysis. To verify the stability of the proposed techniques, we added 1% Gaussian white noise into the raw projection data and repeated the above testing steps. The reconstructed images are also presented in Figure 7, which also support our analysis on the exactness and stability performed in Section III.
In this paper, we have proposed the interior reconstruction technique for SPECT in a parallel-beam geometry. Using a rebinning scheme, we can immediately extend Eq. (5) into a fan-beam geometry and other variants, such as the varying focal-length fan-beam geometry . A major difference lies in the involved Jacobian factor in Eq. (5) in the backprojection step, while the generalized Hilbert transform based reconstruction method remains the same as that for parallel-beam imaging. Using analytic techniques similar to our generalized backprojection filtration method [21, 22], our interior SPECT methodology can be further extended to more complex cone-beam configurations and more general focal spot trajectories.
When we derived the reconstruction formula and analyzed the exactness and stability, we have assumed a uniform attenuation map. Without giving a detailed proof and analysis, it should be pointed out that similar results can be obtained for a non-uniformly attenuated background. In that case, both analytic continuation and SVD techniques remain important tools to analyze the uniqueness, exactness and stability, as well as to reconstruct a distribution of radioactive sources. In addition to the SVD method proposed in this paper, other reconstruction methods, such as POCS , can be also developed for reconstruction from generalized attenuated Hilbert transform data.
As another important nuclear medicine imaging modality, positron emission tomography (PET) produces 3D images of functional and cellular features in the body. Different from SPECT, a radio-active source emits pairs of particles in opposite directions, and they are detected in the coincidence mode, i.e., only events with two particles arriving at opposite detectors within a narrow time window are counted. Thus, the prjoection model of PET can be written as 
After an attenuation correction, the both uniform and non-uniform PET reconstruction can be done as from CT projection data.
In conclusion, by analytic continuation we have extended the interior CT reconstruction technique for CT to that for SPECT. The uniqueness and stability of interior SPECT have been analyzed. Using the SVD technique, we have implemented an interior SPECT algorithm. Both our theoretical analysis and numerical simulations have demonstrated the feasibility of our proposed interior SPECT approach. The theoretical generalization of our work to the variable attenuation case is underway, and the same numerical approach can be applied.
This work is supported in part by NKBRSF (2003CB716101), NIH/NIBIB (EB002667, EB004287, EB007288) and NSFC (60325101, 60532080, 60628102), Ministry of Education of China (306017), Engineering Research Institute of Peking University, Key Laboratory of Machine Perception (Ministry of Education) of Peking University, and Microsoft Research of Asia.