PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Commun Numer Methods Eng. Author manuscript; available in PMC 2010 June 1.
Published in final edited form as:
Commun Numer Methods Eng. 2009 June 1; 25(6): 693–710.
doi:  10.1002/cnm.1206
PMCID: PMC2759671
NIHMSID: NIHMS89644

Interior SPECT- Exact and Stable ROI Reconstruction from Uniformly Attenuated Local Projections

Abstract

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.

Keywords: Single photon emission computed tomography (SPECT), exact interior reconstruction, generalized Hilbert transform, analytic continuation, singular value decomposition (SVD)

I. Introduction

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. [1]. 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. [2] 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 [3], 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 [46] from fan-beam/cone-beam projection datasets, while the conventional wisdom is that the interior problem does not have a unique solution [7]. 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 [1014]. 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 [46, 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.

II. Main Results

Let f (x) be a 2D smooth distribution function defined on a convex compact support Ω with x=(x, y) [set membership] Ω. Mathematically, in a parallel-beam geometry (Fig. 1), the SPECT projections of f (x) can be modeled as the attenuated Radon transform

Figure 1
Parallel beam coordinate system for the SPECT imaging.

Po(θ,s)=f(sθ+tθ)etμ(sθ+tθ)dtdt,
(1)

where the subscript “o” indicates original projection data, θ = (cos θ, sin(θ)), θ[perpendicular] = (−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

μ(x)={μ0xΩ0xΩ,
(2)

where μ0 is a constant. Since the object function is compactly supported, we can determine the maximum coordinate along the direction θ[perpendicular] 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

Po(θ,s)=eμ0tmax(θ,s)f(sθ+tθ)eμ0tdt.
(3)

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

Pw(θ,s)=Po(θ,s)eμ0tmax(θ,s)=f(sθ+tθ)eμ0tdt,
(4)

where the subscript “w” indicates that the weighted projection data.

Denote a weighted backprojection of the differential projection data by

g(x)=0πeμ0xθPw(θ,s)ss=xθdθ.
(5)

In 2004, Rullgard proved a relationship linking the object image f (x) and the weighted backprojection g(x) as [14]

g(x)=2πPV+chμ0(yy)f(x,y)dy,
(6)

where “PV” represents the Cauchy principle value integral, and chμ0 is defined by

chμ0(y)=cosh(μ0y)πy=eμ0y+eμ0y2πy.
(7)

When μ0 → 0, chμ0 becomes the Hilbert transform kernel, in consistence with the results in the CT field [1]. 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 [negated set membership] (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

g(u)=2πPVcsbcsechμ0(uu)f(u)du,u(cvb,cve).
(8)

In 2007, Noo et al. extended their two-step Hilbert CT method to the case of SPECT [13], 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:

Figure 2
Various exact reconstruction conditions based on the generalized Hilbert transform. (a) The exact reconstruction region by Noo et al. [15] and (b) that by ours in this paper.

Theorem 2.1

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

mμ0=csbcsef(u)cosh(μ0u)du.
(9)

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 [4], 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.

III. Uniqueness and Stability Analysis

3.1. Uniqueness Analysis

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

h(u)=hd(u)+1πmμ0+Kh,
(10)

where

h(u)=f(u)1u2,
(11)

hd(u)=PV111u2π(uu)g(u)du,
(12)

(Kh)(u)=11h(u)1u2k¯μ0(u,u)du,
(13)

with

k¯μ0(u,u)=μ0PV111v2π(uv)η(μ0(vu))dv+1π(1coshμ0u),
(14)
η(u)={(coshu1)/(πu)u00u=0.
(15)

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 kμ0(u, ũ) is a smooth and continuous function over the region (u, ũ) [set membership] [−1,1]×[−1,1].

Using the following identical equation

PV111v2π(uv)dv=u,u1,
(16)

kμ0 (u, ũ) can be rewritten as

k¯μ0(u,u)=μ0PV111v2η¯(u,v,u)dv+μ0uη(μ0(uu))+1π(1coshμ0u),
(17)

where

η¯(u,v,u)=η(μ0(vu))η(μ0(uu))π(uv).
(18)

Note that both [eta w/ macron](u, v, ũ) and η(μ0(uũ)) are analytical with respect to the variable u, and kμ0(u, ũ) is analytical with respect to the variable u. Hence, (Kh)(u) can be analytically extended to the whole complex plane C.

Now, let us consider our case −1 = csb < cvb < cke < cve < cse = 1. First, we rewrite Eq. (10) as:

h(u)=h1(u)+h2(u),
(19)

where

h1(u)=PVcvbcve1u2π(uu)g(u)du+1πmμ0,
(20)
h2(u)=PV(1cvb+cve1)(1u2π(uu)g(u)du)+Kh(u),
(21)

By our assumption in Theorem 2.1, h1(u) is known for u [set membership] R and h(u) is known on the interval (cvb, cke), we obtain that

h2(u)=h(u)h1(u)
(22)

is known on the interval (cvb, cke). Because the principal integral of the first term of Eq. (21) is defined on (−1, cvb) [union or logical sum] (cve, 1), h2 (u) is analytic on the interval (cvb, cve), and it can be analytically extended to the complex plane C 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).

3.2. Stability Analysis

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:

  1. |h(u) − hr (u)|≤ ε for u [set membership] (cvb, cke) (f(u) is known in (cvb, cke));
  2. |(hr)1(u) − h1(u) |≤ ε for u [set membership] (cvb, cve);
  3. f(u)M12,fr(u)M12, for u [set membership] (−1,1);
  4. 1π1u2Hf(u)M22,1π1u2Hfr(u)M22, for u [set membership] (−1, cvb)[union or logical sum](cve, 1);
  5. mμ0=11f(u)cosh(μ0u)du=11fr(u)cosh(μ0u)du (mμ0 is known)

The above conditions lead to

  1. |herr (u) |≤ ε for u [set membership] (cvb, cke);
  2. |h1,err (u))|≤ ε for u [set membership] (cvb,cve);
  3. |ferr (u) |≤ M1 for u [set membership] (−1,1);
  4. 1π1u2Hferr(u)M2 for u [set membership] (−1,cvb) [union or logical sum] (cve,1);

where

herr(u)=hr(u)h(u),h1,err(u)=(hr)1(u)h1(u),ferr(u)=fr(u)f(u).
(23)

Similarly, we denote

h2,err(u)=(hr)2(u)h2(u).
(24)

Choose two arbitrary constants cvb, cve such that cvb < cvb < cke < cve < cve. On the interval (cvb, cve), we have

herr(u)=h1,err(u)+h2,err(u).
(25)

Especially, on the interval (cvb, cke), we have

h2,err(u)h1,err(u)+herr(u)2ε.
(26)

In order to obtain a upper bound of |h2,err (u)| on the interval (cke, cve), we will use following Nevanlinna’s principle [2, 17].

Lemma (Nevanlinna’s principle)

Let Ω [subset or is implied by] C be a domain and D be a segment of the boundary [partial differential]Ω. If f (z) is an analytical function in Ω such that

f(z){εforzDMforzΩ\D,

and ω(z) is a harmonic function in Ω such that

ω(z)={0forzD1forzΩ\D,

then

f(z)Mω(z)ε1ω(z).

Let B to be the disk with diameter (Cvb, cve), Ω = B\(cvb, cke), then [partial differential]Ω= D + [partial differential]B, where D = (cvb, cke), and [partial differential]B is the circle with diameter (cvb, cve) (Figure 3). Recall that

Figure 3
Interior reconstruction configuration in terms of Nevanlinna principle.
h2,err(u)=PV(1cvb+cve1)(1u2π(uu)Hferr(u)du)+Kherr(u).
(27)

For z [set membership] [partial differential]B,

h2,err(z)(1cvb+cve1)(1u2π(zu)Hferr(u)du)+Kherr(z),
(28)

where

(1cvb+cve1)(1u2π(zu)Hferr(u)du)1cvb1u2Hferr(u)π(Rezu)du+cve11u2Hferr(u)π(uRez)du1cvbM2Rezudu+cve1M2uRezdu=M2(ln(Rez+1Rezcvb)ln(1RezcveRez))M2(ln(cvb+1cvbcvb)+ln(1cvecvecve)),
(29)

and

Kherr(z)11ferr(u)k¯μ0(z,u)duCM1,
(30)

with C=maxz1,u1k¯μ0(z,u).

We conclude that

(i)h2,err(z)M2(ln(cvb+1cvbcvb)+ln(1cvecvecve))+CM1forzΩ\D,
(31)
(ii)h2,err(z)2εforzD.
(32)

On other hand, there exits a harmonic function ω(z) satisfying [16]

ω(z)={0forzD1forzΩ\D,

and for u [set membership] (cke, cve),

ω(u)=4πarctan2(ucke)(cvecvb)(cvecvb)2(2ckecvecvb)(2ucvecvb).
(33)

Therefore, by Nevanlinna’s principle, for u [set membership] (cke, cve)

h2,err(u){M2(ln(cvb+1cvbcvb)+ln(1cvecvecve))+CM1}ω(u)(2ε)1ω(u),
(34)

and

herr(u)ε+{M2(ln(cvb+1cvbcvb)+ln(1cvecvecve))+CM1}ω(u)(2ε)1ω(u).
(35)

One remark is in order that for u [set membership] (cke, cve) and 0 < ω(u) < 1, ω(u) → 1 as ucve). Eq. (35) implies that the reconstruction of f(u) for u [set membership](cke, cve) is stable near cke and probably less stable near cve.

IV. Numerical Implementation and Simulation

4.1 Singular Value Decomposition

Let us rewrite Eq. (8) as

gw(u)=g(u)2π=PVcsbcsechμ0(uu)f(u)du,u(cvb,cve),
(36)

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

G=AF,
(37)

where

F=[f1f2fnfN]T,
(38)
G=[g1g2gmgM]T,
(39)
A=[A1,1A1,nA1,NAm,1Am,nAm,NAM,1AM,nAM,N].
(40)

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 [18], we can write Am, n as

Am,n={eμ0um,nd+eμ0um,ndπum,nd/δum,nd/δisodd0um,nd/δiseven,
(41)

with

um,nd=u(gm)u(fn).
(42)

With our sampling method, um,nd/δ always produces integers, and the case um,nd/δ=0 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

F=(FkFu),
(43)

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

A=(AkAu).
(44)

Then, we immediately arrive at a linear equation system

A¯F¯=G¯,
(45)

with G = GAkFk, Ā = Au and [F with macron] = Fu.

Although Eq. (45) is generally ill-posed, we can still stably reconstruct the part of [F with macron] in the interval (cke, cve) using the singular value decomposition (SVD) method [19, 20]. Note that the dimension of G is [M with macron] = M. If the dimension of [F with macron] is N, the dimension of Ā will be [M with macron] × N. According to the SVD theory [19, 20], the matrix Ā has a SVD decomposition in the following form

A¯=UΛVT,
(46)

where U and V are orthogonal matrices of [M with macron] × [M with macron] and N × N respectively, and Λ is an [M with macron] × N diagonal matrix whose diagonal elements λq satisfying λ1≥ … ≥ λq ≥ … ≥ λQ, Q = min([M with macron], N). In this way, we have a stable numerical solution for [F with macron] [19, 20],

F¯^=F^u=VΛ1UTG¯,
(47)

where Λ−1 is a diagonal matrix of N × [M with macron] whose diagonal elements λq1 are defined as

λq1={1/λqλq>ε0λqε,
(48)

and ε > 0 is a free small constant parameter. The system defined by Eqs. (47) and (48) is called a truncated SVD (TSVD) [19], which is a special case of regularization.

4.2 Numerical Analysis

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 [7]. Hence, we can stably reconstruct f (u) at its finite sampling points on the interval (cvb, cve) using the aforementioned SVD method.

Recall that G = Ā[F with macron] = AuFu, we have

F^uFu=VΛ1UTAuFuFu=(VΛ1UTAuIu)Fu,
(49)

where Iu represents an N × N unit diagonal matrix. Let Eu = FuFu as the error vector of Fu, and Su = −1UT AuIu as an N × N matrix. Then, Eq. (49) can be reduced as

Eu=SuFu.
(50)

For any 1 ≤ [n with macron]N, the absolute reconstruction error can be expressed as

En¯u=n=1N¯Sn¯,nuFnun=1N¯Sn¯,nuFnu(n=1N¯Sn¯,nu)Fmaxu=Cn¯EuFmaxu,
(51)

where Fmaxu represents the maximum absolute value of the elements in Fu, and Cn¯Eu=n=1N¯Sn¯,nu. Because Cn¯Eu reflects the accuracy of the reconstructed Fu, we call it the precision measure in this paper.

Next, let us consider the case of noisy data. Let us assume that

Gn=G+Wg,Wmg<ε,
(52)

Fnk=Fk+Wfk,Wnfk<ε,
(53)

where ε defines the maximum noise magnitude. Recall that G = GAkFk, we have the noise expression of G as

G¯n=WgAkWfk.
(54)

Finally, based on Eq. (47), the noise in the reconstructed image can be expressed as

F^un=VΛ1UT(WgAkWfk),
(55)

which can be further reduced as

F^un=PunWun,
(56)

with Pun = (−1UT−1UTAk) and Wun=(WgWfk). Clearly, the dimension of Fun is N × 1. Assume that the dimension of Wun is [M with tilde] × 1, the dimension of Pun must be N × [M with tilde]. For any 1 ≤ [n with macron]N, the absolute magnitude of noise in the image can be further expressed as:

F^n¯un=m=1MPn¯,munWmunm=1MPn¯,munWmunM=1mPn¯,munε=Cn¯Nuε,
(57)

where Cn¯Nu=m=1MPn¯,mun and Cn¯Nu is called the stability measure in this paper. By evaluating the values of Cn¯Nu, 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 Cn¯Eu and Cn¯Nu as shown in Figures 4 and and5,5, respectively. Because f (u) was known on (cvb, cke), the corresponding portions of Cn¯Eu and Cn¯Nu 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 Cn¯Eu 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 Cn¯Eu would decrease in the interval [cke, cve), which means a higher accuracy. Undoubtedly, Cn¯Eu 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.

Figure 4
Exactness measurement Cn¯Eu for different configurations of N, μ0 and ε.
Figure 5
Stability measurement Cn¯Nu for different configurations of N, μ0 and ε.

4.3 Numerical Simulation

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. [13]. 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).

Figure 6
Shepp-Logan phantom for interior reconstruction of SPECT. (a) Modified sheep-Logan phantom in a display window [0.15,0.45]. (b) Interior reconstruction configuration including a disk support, a rectangular field of view and a strip-shaped known prior ...

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.

Figure 7
SPECT interior reconstruction results. The left, middle and right columns are respectively for the case μ0 =0 cm−1, 0.15 cm−1 and 0.30 cm−1. The top row images are reconstructed from noise-free projection data, while the ...

V. Discussions and Conclusion

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 [10]. 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 [4], 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 [23]

Po(θ,s)=e+μ(sθ+tθ)dtf(sθ+tθ)dt.
(58)

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.

Acknowledgments

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.

References

1. Noo F, Clackdoyle R, Pack JD. A two-step Hilbert transform method for 2D image reconstruction. Physics In Medicine And Biology. 2004;49(17):3903–3923. [PubMed]
2. Defrise M, et al. Truncated Hilbert transform and image reconstruction from limited tomographic data. Inverse Problems. 2006;22(3):1037–1053.
3. Wang G, Ye Y, Yu H. General VOI/ROI Reconstruction Methods and Systems using a Truncated Hilbert Transform (Patent disclosure submitted to Virginia Tech. Intellectual Properties on May 15, 2007. 2007
4. Ye Y, et al. A General Local Reconstruction Approach Based on a Truncated Hilbert Transform. International Journal of Biomedical Imaging, 2007. 2007:8. Article ID: 63634. [PMC free article] [PubMed]
5. Ye Y, Yu H, Wang G. Exact Interior Reconstruction with Cone-beam CT. International Journal of Biomedical Imaging, 2007. 2007:5. Article ID: 10693. [PMC free article] [PubMed]
6. Ye Y, Yu HY, Wang G. Exact Interior Reconstruction from truncated limited-angle projection data. International Journal of Biomedical Imaging, 2008. 2008:6. Article ID: 427989. [PMC free article] [PubMed]
7. Natterer F. Classics in applied mathematics. Vol. 32. Philadelphia: Society for Industrial and Applied Mathematics; 2001. The Mathematics of Computerized Tomography.
8. Kudo H, et al. Tiny a priori knowledge solves the interior problem. 2007 IEEE Nuclear Science Symposium and Medical Imaging; 2007. p. Paper No.: M21-1.
9. Kudo H, et al. Tiny a priori knowledge solves the interior problem in computed tomography. Phys Med Biol. 2008;53(9):2207–2231. [PubMed]
10. Li T, et al. An efficient reconstruction method for nonuniform attenuation compensation in nonparallel beam geometries based on Novikov’s explicit inversion formula. IEEE Transactions on Medical Imaging. 2005;24(10):1357–1368. [PMC free article] [PubMed]
11. Tang QL, et al. Exact fan-beam and 4 pi-acquisition cone-beam SPECT algorithms with uniform attenuation correction. Medical Physics. 2005;32(11):3440–3447. [PubMed]
12. Tang QL, Zeng GSL, Gullberg GT. Analytical fan-beam and cone-beam reconstruction algorithms with uniform attenuation correction for SPECT. Physics in Medicine and Biology. 2005;50(13):3153–3170. [PubMed]
13. Noo F, et al. Image reconstruction from truncated data in single-photon emission computed tomography with uniform attenuation. Inverse Problems. 2007;23(2):645–667.
14. Rullgard H. An explicit inversion formula for the exponential Radon transform using data from 180 ~. Ark Mat. 2004;42:353–362.
15. Rullgard H. Stability of the inverse problem for the attenuated Radon transform with 180 degrees data. Inverse Problems. 2004;20(3):781–797.
16. Courdurier M. Restricted Measurements for the X-ray Transform. University of Washington-Seattle; 2007.
17. Nevanlinna R, editor. Grundlehren der mathematischen Wissenschaften. 2. Vol. 46. 1953. Eindentige Analytische Functionen.
18. Yu H, Wang G. Studies on implementation of the Katsevich algorithm for spiral cone-beam CT. Journal of X-ray Science and Technology. 2004;12(2):96–117.
19. Hansen PC. Truncated Singular Value Decomposition Solutions To Discrete Ill-Posed Problems With Ill-Determined Numerical Rank. Siam Journal On Scientific And Statistical Computing. 1990;11(3):503–518.
20. Hansen PC. Regularization, Gsvd And Truncated Gsvd. Bit. 1989;29(3):491–504.
21. Ye Y, et al. Exact reconstruction for cone-beam scanning along nonstandard spirals and other curves. Developments in X-Ray Tomography IV, Proceedings of SPIE VOL; Aug. 4–6. 2004; Denver, CO, United States. pp. 293–300.
22. Ye YB, et al. A General Exact Reconstruction for Cone-Beam CT via Backprojection-Filtration. IEEE Transactions on Medical Imaging. 2005;24(9):1190–1198. [PubMed]
23. Natterer F, Wubbeling F, editors. Mathematical Methods in Image Reconstruction. Society for Industrial and Applied Mathematics; Philadelphia: 2001.