Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2759671

Formats

Article sections

- Abstract
- I. Introduction
- II. Main Results
- III. Uniqueness and Stability Analysis
- IV. Numerical Implementation and Simulation
- V. Discussions and Conclusion
- References

Authors

Related links

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.1206PMCID: PMC2759671

NIHMSID: NIHMS89644

See other articles in PMC that cite the published article.

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.* [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 [4–6] 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 [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

$${P}_{o}(\theta ,s)=\underset{-\infty}{\overset{\infty}{\int}}f(s\mathbf{\theta}+t{\mathbf{\theta}}^{\perp}){e}^{-\underset{t}{\overset{\infty}{\int}}\mu (s\mathbf{\theta}+\stackrel{\sim}{t}{\mathbf{\theta}}^{\perp})d\stackrel{\sim}{t}}dt,$$

(1)

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 *P _{o}* (

$$\mu (\mathbf{x})=\{\begin{array}{cc}\hfill {\mu}_{0}\hfill & \hfill \mathbf{x}\in \mathrm{\Omega}\hfill \\ \hfill 0\hfill & \hfill \mathbf{x}\notin \mathrm{\Omega}\hfill \end{array},$$

(2)

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 *P _{o}* (

$${P}_{o}(\theta ,s)={e}^{-{\mu}_{0}{t}_{max}(\theta ,s)}\underset{-\infty}{\overset{\infty}{\int}}f(s\mathbf{\theta}+t{\mathbf{\theta}}^{\perp}){e}^{{\mu}_{0}t}dt.$$

(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 *t*_{max} (*θ*, *s*) for all the whole projection dataset. By multiplying a weighting factor *e*^{μ0tmax (θ,s)}, the projection model of SPECT can be reduced to

$${P}_{w}(\theta ,s)={P}_{o}(\theta ,s){e}^{{\mu}_{0}{t}_{max}(\theta ,s)}=\underset{-\infty}{\overset{\infty}{\int}}f(s\mathbf{\theta}+t{\mathbf{\theta}}^{\perp}){e}^{{\mu}_{0}t}dt,$$

(4)

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

Denote a weighted backprojection of the differential projection data by

$$g(\mathbf{x})=\underset{0}{\overset{\pi}{\int}}{{e}^{-{\mu}_{0}\mathbf{x}\u2022{\mathbf{\theta}}^{\perp}}\frac{{\partial P}_{w}(\theta ,s)}{\partial s}\mid}_{s=\mathbf{x}\u2022\mathbf{\theta}}d\theta .$$

(5)

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

$$g(\mathbf{x})=-2\pi PV\underset{-\infty}{\overset{+\infty}{\int}}c{h}_{{\mu}_{0}}(y-\stackrel{\sim}{y})f(x,\stackrel{\sim}{y})d\stackrel{\sim}{y},$$

(6)

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

$$c{h}_{{\mu}_{0}}(y)=\frac{cosh({\mu}_{0}y)}{\pi y}=\frac{{e}^{{\mu}_{0}y}+{e}^{-{\mu}_{0}y}}{2\pi 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 (*c _{sb}*,

$$g(u)=-2\pi PV\underset{{c}_{sb}}{\overset{{c}_{se}}{\int}}c{h}_{{\mu}_{0}}(u-\stackrel{\sim}{u})f(\stackrel{\sim}{u})d\stackrel{\sim}{u},\phantom{\rule{0.38889em}{0ex}}u\in ({c}_{vb},{c}_{ve}).$$

(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 (*c _{sb}*,

Assume that *c _{sb}* <

$${m}_{{\mu}_{0}}=\underset{{c}_{sb}}{\overset{{c}_{se}}{\int}}f(\stackrel{\sim}{u})cosh({\mu}_{0}\stackrel{\sim}{u})d\stackrel{\sim}{u}.$$

(9)

Note that we have omitted a parameter *c _{kb}* =

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 (*c _{vb}*,

Without loss of generality, we assume *c _{sb}* = −1 and

$$h(u)={h}_{d}(u)+\frac{1}{\pi}{m}_{{\mu}_{0}}+Kh,$$

(10)

where

$$h(u)=f(u)\sqrt{1-{u}^{2}},$$

(11)

$${h}_{d}(u)=-PV\underset{-1}{\overset{1}{\int}}\frac{\sqrt{1-{\stackrel{\sim}{u}}^{2}}}{\pi (u-\stackrel{\sim}{u})}g(\stackrel{\sim}{u})d\stackrel{\sim}{u},$$

(12)

$$(Kh)(u)=\underset{-1}{\overset{1}{\int}}\frac{h(\stackrel{\sim}{u})}{\sqrt{1-{\stackrel{\sim}{u}}^{2}}}{\overline{k}}_{{\mu}_{0}}(u,\stackrel{\sim}{u})d\stackrel{\sim}{u},$$

(13)

with

$${\overline{k}}_{{\mu}_{0}}(u,\stackrel{\sim}{u})={\mu}_{0}PV\underset{-1}{\overset{1}{\int}}\frac{\sqrt{1-{v}^{2}}}{\pi (u-v)}\eta ({\mu}_{0}(v-\stackrel{\sim}{u}))dv+\frac{1}{\pi}(1-cosh{\mu}_{0}\stackrel{\sim}{u}),$$

(14)

$$\eta (u)=\{\begin{array}{cc}(coshu-1)/(\pi u)& u\ne 0\\ 0& u=0\end{array}.$$

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

Using the following identical equation

$$PV\underset{-1}{\overset{1}{\int}}\frac{\sqrt{1-{v}^{2}}}{\pi (u-v)}dv=u,\phantom{\rule{0.38889em}{0ex}}\mid u\mid \le 1,$$

(16)

_{μ0} (*u*, *ũ*) can be rewritten as

$${\overline{k}}_{{\mu}_{0}}(u,\stackrel{\sim}{u})={\mu}_{0}PV\underset{-1}{\overset{1}{\int}}\sqrt{1-{v}^{2}}\overline{\eta}(u,v,\stackrel{\sim}{u})dv+{\mu}_{0}u\eta ({\mu}_{0}(u-\stackrel{\sim}{u}))+\frac{1}{\pi}(1-cosh{\mu}_{0}\stackrel{\sim}{u}),$$

(17)

where

$$\overline{\eta}(u,v,\stackrel{\sim}{u})=\frac{\eta ({\mu}_{0}(v-\stackrel{\sim}{u}))-\eta ({\mu}_{0}(u-\stackrel{\sim}{u}))}{\pi (u-v)}.$$

(18)

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 = *c _{sb}* <

$$h(u)={h}_{1}(u)+{h}_{2}(u),$$

(19)

where

$${h}_{1}(u)=-PV\underset{{c}_{vb}}{\overset{{c}_{ve}}{\int}}\frac{\sqrt{1-{\stackrel{\sim}{u}}^{2}}}{\pi (u-\stackrel{\sim}{u})}g(\stackrel{\sim}{u})d\stackrel{\sim}{u}+\frac{1}{\pi}{m}_{{\mu}_{0}},$$

(20)

$${h}_{2}(u)=-PV\left(\underset{-1}{\overset{{c}_{vb}}{\int}}+\underset{{c}_{ve}}{\overset{1}{\int}}\phantom{\rule{0.16667em}{0ex}}\right)\left(\frac{\sqrt{1-{\stackrel{\sim}{u}}^{2}}}{\pi (u-\stackrel{\sim}{u})}g(\stackrel{\sim}{u})d\stackrel{\sim}{u}\right)+Kh(u),$$

(21)

By our assumption in Theorem 2.1, *h*_{1}(*u*) is known for *u* *R* and *h*(*u*) is known on the interval (*c _{vb}*,

$${h}_{2}(u)=h(u)-{h}_{1}(u)$$

(22)

is known on the interval (*c _{vb}*,

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 *g _{r}*(

- |
*h*(*u*) −*h*(_{r}*u*)|≤*ε*for*u*(*c*,_{vb}*c*) (_{ke}*f*(*u*) is known in (*c*,_{vb}*c*));_{ke} - |(
*h*)_{r}_{1}(*u*) −*h*_{1}(*u*) |≤*ε*for*u*(*c*,_{vb}*c*);_{ve} - $\mid f(u)\mid \le \frac{{M}_{1}}{2},\mid {f}_{r}(u)\mid \le \frac{{M}_{1}}{2}$, for
*u*(−1,1); - $\frac{1}{\pi}\sqrt{1-{u}^{2}}\mid Hf(u)\mid \le \frac{{M}_{2}}{2},\frac{1}{\pi}\sqrt{1-{u}^{2}}\mid {Hf}_{r}(u)\mid \le \frac{{M}_{2}}{2}$, for
*u*(−1,*c*)(_{vb}*c*, 1);_{ve} - ${m}_{{\mu}_{0}}=\underset{-1}{\overset{1}{\int}}f(\stackrel{\sim}{u})cosh({\mu}_{0}\stackrel{\sim}{u})d\stackrel{\sim}{u}=\underset{-1}{\overset{1}{\int}}{f}_{r}(\stackrel{\sim}{u})cosh({\mu}_{0}\stackrel{\sim}{u})d\stackrel{\sim}{u}$ (
*m*_{μ0}is known)

The above conditions lead to

- |
*h*(_{err}*u*) |≤*ε*for*u*(*c*,_{vb}*c*);_{ke} - |
*h*_{1,}(_{err}*u*))|≤*ε*for*u*(*c*,_{vb}*c*);_{ve} - |
*f*(_{err}*u*) |≤*M*_{1}for*u*(−1,1); - $\frac{1}{\pi}\sqrt{1-{u}^{2}}\mid H{f}_{\mathit{err}}(u)\mid \le {M}_{2}$ for
*u*(−1,*c*) (_{vb}*c*,1);_{ve}

where

$${h}_{\mathit{err}}(u)={h}_{r}(u)-h(u),{h}_{1,\mathit{err}}(u)={({h}_{r})}_{1}(u)-{h}_{1}(u),{f}_{\mathit{err}}(u)={f}_{r}(u)-f(u).$$

(23)

Similarly, we denote

$${h}_{2,\mathit{err}}(u)={({h}_{r})}_{2}(u)-{h}_{2}(u).$$

(24)

Choose two arbitrary constants * _{vb}*,

$${h}_{\mathit{err}}(u)={h}_{1,\mathit{err}}(u)+{h}_{2,\mathit{err}}(u).$$

(25)

Especially, on the interval (* _{vb}*,

$$\mid {h}_{2,\mathit{err}}(u)\mid \le \mid {h}_{1,\mathit{err}}(u)\mid +\mid {h}_{\mathit{err}}(u)\mid \le 2\epsilon .$$

(26)

In order to obtain a upper bound of |*h*_{2,}* _{err}* (

Let Ω be a domain and *D* be a segment of the boundary Ω. If *f* (*z*) is an analytical function in Ω such that

$$\mid f(z)\mid \le \{\begin{array}{l}\epsilon \phantom{\rule{0.16667em}{0ex}}\text{for}\phantom{\rule{0.16667em}{0ex}}z\in D\hfill \\ M\phantom{\rule{0.16667em}{0ex}}\text{for}\phantom{\rule{0.16667em}{0ex}}z\in \partial \mathrm{\Omega}\backslash D\hfill \end{array},$$

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

$$\omega (z)=\{\begin{array}{ll}0\hfill & \text{for}\phantom{\rule{0.16667em}{0ex}}z\in D\hfill \\ 1\hfill & \text{for}\phantom{\rule{0.16667em}{0ex}}z\in \partial \mathrm{\Omega}\backslash D\hfill \end{array},$$

then

$$\mid f(z)\mid \le {M}^{\omega (z)}{\epsilon}^{1-\omega (z)}.$$

Let *B* to be the disk with diameter (* _{vb}*,

$${h}_{2,\mathit{err}}(u)=-PV\left(\underset{-1}{\overset{{c}_{vb}}{\int}}+\underset{{c}_{ve}}{\overset{1}{\int}}\phantom{\rule{0.16667em}{0ex}}\right)\left(\frac{\sqrt{1-{\stackrel{\sim}{u}}^{2}}}{\pi (u-\stackrel{\sim}{u})}H{f}_{\mathit{err}}(\stackrel{\sim}{u})d\stackrel{\sim}{u}\right)+K{h}_{\mathit{err}}(u).$$

(27)

For *z* *B*,

$$\mid {h}_{2,\mathit{err}}(z)\mid \le \mid \left(\underset{-1}{\overset{{c}_{vb}}{\int}}+\underset{{c}_{ve}}{\overset{1}{\int}}\phantom{\rule{0.16667em}{0ex}}\right)\left(\frac{\sqrt{1-{\stackrel{\sim}{u}}^{2}}}{\pi (z-\stackrel{\sim}{u})}H{f}_{\mathit{err}}(\stackrel{\sim}{u})d\stackrel{\sim}{u}\right)\mid +\mid K{h}_{\mathit{err}}(z)\mid ,$$

(28)

where

$$\begin{array}{l}\mid \left(\underset{-1}{\overset{{c}_{vb}}{\int}}+\underset{{c}_{ve}}{\overset{1}{\int}}\phantom{\rule{0.16667em}{0ex}}\right)\left(\frac{\sqrt{1-{\stackrel{\sim}{u}}^{2}}}{\pi (z-\stackrel{\sim}{u})}H{f}_{\mathit{err}}(\stackrel{\sim}{u})d\stackrel{\sim}{u}\right)\mid \le \underset{-1}{\overset{{c}_{vb}}{\int}}\frac{\mid \sqrt{1-{\stackrel{\sim}{u}}^{2}}H{f}_{\mathit{err}}(\stackrel{\sim}{u})\mid}{\pi (Rez-\stackrel{\sim}{u})}d\stackrel{\sim}{u}+\underset{{c}_{ve}}{\overset{1}{\int}}\frac{\mid \sqrt{1-{\stackrel{\sim}{u}}^{2}}H{f}_{\mathit{err}}(\stackrel{\sim}{u})\mid}{\pi (\stackrel{\sim}{u}-Rez)}d\stackrel{\sim}{u}\\ \le \underset{-1}{\overset{{c}_{vb}}{\int}}\frac{{M}_{2}}{Rez-\stackrel{\sim}{u}}d\stackrel{\sim}{u}+\underset{{c}_{ve}}{\overset{1}{\int}}\frac{{M}_{2}}{\stackrel{\sim}{u}-Rez}d\stackrel{\sim}{u}={M}_{2}\left(ln\left(\frac{Rez+1}{Rez-{c}_{vb}}\right)-ln\left(\frac{1-Rez}{{c}_{ve}-Rez}\right)\right)\\ \le {M}_{2}\left(ln\left(\frac{{\stackrel{\sim}{c}}_{vb}+1}{{\stackrel{\sim}{c}}_{vb}-{c}_{vb}}\right)+ln\left(\frac{1-{\stackrel{\sim}{c}}_{ve}}{{c}_{ve}-{\stackrel{\sim}{c}}_{ve}}\right)\right),\end{array}$$

(29)

and

$$\mid K{h}_{\mathit{err}}(z)\mid \le \mid \underset{-1}{\overset{1}{\int}}{f}_{\mathit{err}}(\stackrel{\sim}{u}){\overline{k}}_{{\mu}_{0}}(z,\stackrel{\sim}{u})d\stackrel{\sim}{u}\mid \le C{M}_{1},$$

(30)

with $C=\underset{\mid z\mid \le 1,\mid \stackrel{\sim}{u}\mid \le 1}{max}\mid {\overline{k}}_{{\mu}_{0}}(z,\stackrel{\sim}{u})\mid $.

We conclude that

$$(\text{i})\phantom{\rule{0.16667em}{0ex}}\mid {h}_{2,\mathit{err}}(z)\mid \le {M}_{2}\left(ln\left(\frac{{\stackrel{\sim}{c}}_{vb}+1}{{\stackrel{\sim}{c}}_{vb}-{c}_{vb}}\right)+ln\left(\frac{1-{\stackrel{\sim}{c}}_{ve}}{{c}_{ve}-{\stackrel{\sim}{c}}_{ve}}\right)\right)+C{M}_{1}\phantom{\rule{0.16667em}{0ex}}\text{for}\phantom{\rule{0.16667em}{0ex}}z\in \partial \mathrm{\Omega}\backslash D,$$

(31)

$$(\text{ii})\phantom{\rule{0.16667em}{0ex}}\mid {h}_{2,\mathit{err}}(z)\mid \le 2\epsilon \phantom{\rule{0.16667em}{0ex}}\text{for}\phantom{\rule{0.16667em}{0ex}}z\in D.$$

(32)

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

$$\omega (z)=\{\begin{array}{ll}0\hfill & \text{for}\phantom{\rule{0.16667em}{0ex}}z\in D\hfill \\ 1\hfill & \text{for}\phantom{\rule{0.16667em}{0ex}}z\in \partial \mathrm{\Omega}\backslash D\hfill \end{array},$$

and for *u* (*c _{ke}*,

$$\omega (u)=\frac{4}{\pi}arctan\sqrt{\frac{2(u-{c}_{ke})({\stackrel{\sim}{c}}_{ve}-{\stackrel{\sim}{c}}_{vb})}{{({\stackrel{\sim}{c}}_{ve}-{\stackrel{\sim}{c}}_{vb})}^{2}-(2{c}_{ke}-{\stackrel{\sim}{c}}_{ve}-{\stackrel{\sim}{c}}_{vb})(2u-{\stackrel{\sim}{c}}_{ve}-{\stackrel{\sim}{c}}_{vb})}}.$$

(33)

Therefore, by Nevanlinna’s principle, for *u* (*c _{ke}*,

$$\phantom{\rule{0.16667em}{0ex}}\mid {h}_{2,\mathit{err}}(u)\mid \le {\left\{{M}_{2}\left(ln\left(\frac{{\stackrel{\sim}{c}}_{vb}+1}{{\stackrel{\sim}{c}}_{vb}-{c}_{vb}}\right)+ln\left(\frac{1-{\stackrel{\sim}{c}}_{ve}}{{c}_{ve}-{\stackrel{\sim}{c}}_{ve}}\right)\right)+C{M}_{1}\right\}}^{\omega (u)}{(2\epsilon )}^{1-\omega (u)},$$

(34)

and

$$\phantom{\rule{0.16667em}{0ex}}\mid {h}_{\mathit{err}}(u)\mid \le \epsilon +{\left\{{M}_{2}\left(ln\left(\frac{{\stackrel{\sim}{c}}_{vb}+1}{{\stackrel{\sim}{c}}_{vb}-{c}_{vb}}\right)+ln\left(\frac{1-{\stackrel{\sim}{c}}_{ve}}{{c}_{ve}-{\stackrel{\sim}{c}}_{ve}}\right)\right)+C{M}_{1}\right\}}^{\omega (u)}{(2\epsilon )}^{1-\omega (u)}.$$

(35)

One remark is in order that for *u* (*c _{ke}*,

Let us rewrite Eq. (8) as

$${g}_{w}(u)=-\frac{g(u)}{2\pi}=PV\underset{{c}_{sb}}{\overset{{c}_{se}}{\int}}c{h}_{{\mu}_{0}}(u-\stackrel{\sim}{u})f(\stackrel{\sim}{u})d\stackrel{\sim}{u},\phantom{\rule{0.38889em}{0ex}}u\in ({c}_{vb},{c}_{ve}),$$

(36)

and discretize the *u* -axis with a uniform sampling interval *δ* satisfying sampling theorem for both *f*(*u*) and *g _{w}*(

$$\mathbf{G}=\mathbf{AF},$$

(37)

where

$$\mathbf{F}={[{f}_{1}{f}_{2}\cdots {f}_{n}\cdots {f}_{N}]}^{T},$$

(38)

$$\mathbf{G}={[{g}_{1}{g}_{2}\cdots {g}_{m}\cdots {g}_{M}]}^{T},$$

(39)

$$\mathbf{A}=\left[\begin{array}{ccccc}{A}_{1,1}& \cdots & {A}_{1,n}& \cdots & {A}_{1,N}\\ \vdots & \ddots & \vdots & \ddots & \vdots \\ {A}_{m,1}& \cdots & {A}_{m,n}& \cdots & {A}_{m,N}\\ \vdots & \ddots & \vdots & \ddots & \vdots \\ {A}_{M,1}& \cdots & {A}_{M,n}& \cdots & {A}_{M,N}\end{array}\right].$$

(40)

In Eqs. (38) and (39), “*T*” represents the transpose operator, and in Eq. (40) *A _{m}*

$${A}_{m,n}=\{\begin{array}{cc}\frac{{e}^{{\mu}_{0}{u}_{m,n}^{d}}+{e}^{-{\mu}_{0}{u}_{m,n}^{d}}}{\pi {u}_{m,n}^{d}/\delta}& {u}_{m,n}^{d}/\delta \phantom{\rule{0.16667em}{0ex}}\text{is}\phantom{\rule{0.16667em}{0ex}}\text{odd}\\ 0& {u}_{m,n}^{d}/\delta \phantom{\rule{0.16667em}{0ex}}\text{is}\phantom{\rule{0.16667em}{0ex}}\text{even}\end{array},$$

(41)

with

$${u}_{m,n}^{d}=u({g}_{m})-u({f}_{n}).$$

(42)

With our sampling method,
${u}_{m,n}^{d}/\delta $ always produces integers, and the case
${u}_{m,n}^{d}/\delta =0$ in Eq. (42) exactly corresponds to the singular point in the Cauchy principal integral in Eq. (36). Because *f*(*u*) is known on (*c _{vb}*,

$$\mathbf{F}=\left(\begin{array}{l}{\mathbf{F}}^{k}\hfill \\ {\mathbf{F}}^{u}\hfill \end{array}\right),$$

(43)

where **F*** ^{k}* is the known part on the sampling points in (

$$\mathbf{A}=\left(\begin{array}{ll}{\mathbf{A}}^{k}\hfill & {\mathbf{A}}^{u}\hfill \end{array}\right).$$

(44)

Then, we immediately arrive at a linear equation system

$$\overline{\mathbf{A}}\overline{\mathbf{F}}=\overline{\mathbf{G}},$$

(45)

with = **G** − **A**^{k}**F*** ^{k}*,

Although Eq. (45) is generally ill-posed, we can still stably reconstruct the part of in the interval (*c _{ke}*,

$$\overline{\mathbf{A}}=\mathbf{U}\mathbf{\Lambda}{\mathbf{V}}^{T},$$

(46)

where **U** and **V** are orthogonal matrices of × and × respectively, and **Λ** is an × diagonal matrix whose diagonal elements *λ _{q}* satisfying

$$\widehat{\overline{\mathbf{F}}}={\widehat{\mathbf{F}}}^{u}=\mathbf{V}{\mathbf{\Lambda}}^{-1}{\mathbf{U}}^{T}\overline{\mathbf{G}},$$

(47)

where **Λ**^{−1} is a diagonal matrix of × whose diagonal elements
${\lambda}_{q}^{-1}$ 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 [7]. Hence, we can stably reconstruct *f* (*u*) at its finite sampling points on the interval (*c _{vb}*,

Recall that = **Ā** = **A**^{u}**F*** ^{u}*, we have

$${\widehat{\mathbf{F}}}^{u}-{\mathbf{F}}^{u}=\mathbf{V}{\mathbf{\Lambda}}^{-1}{\mathbf{U}}^{T}{\mathbf{A}}^{u}{\mathbf{F}}^{u}-{\mathbf{F}}^{u}=\left(\mathbf{V}{\mathbf{\Lambda}}^{-1}{\mathbf{U}}^{T}{\mathbf{A}}^{u}-{\mathbf{I}}^{u}\right){\mathbf{F}}^{u},$$

(49)

where **I*** ^{u}* represents an × unit diagonal matrix. Let

$${\mathbf{E}}^{u}={\mathbf{S}}^{u}{\mathbf{F}}^{u}.$$

(50)

For any 1 ≤ ≤ , the absolute reconstruction error can be expressed as

$$\mid {E}_{\overline{n}}^{u}\mid =\mid \sum _{{n}^{\prime}=1}^{\overline{N}}{S}_{\overline{n},{n}^{\prime}}^{u}{F}_{{n}^{\prime}}^{u}\mid \le \sum _{{n}^{\prime}=1}^{\overline{N}}\mid {S}_{\overline{n},{n}^{\prime}}^{u}\mid \mid {F}_{{n}^{\prime}}^{u}\mid \le \left(\sum _{{n}^{\prime}=1}^{\overline{N}}\mid {S}_{\overline{n},{n}^{\prime}}^{u}\mid \right){F}_{max}^{u}={C}_{\overline{n}}^{Eu}{F}_{max}^{u},$$

(51)

where
${F}_{max}^{u}$ represents the maximum absolute value of the elements in **F*** ^{u}*, and
${C}_{\overline{n}}^{Eu}={\displaystyle \sum _{{n}^{\prime}=1}^{\overline{N}}}\mid {S}_{\overline{n},{n}^{\prime}}^{u}\mid $. Because
${C}_{\overline{n}}^{Eu}$ reflects the accuracy of the reconstructed

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

$${\mathbf{G}}^{n}=\mathbf{G}+{\mathbf{W}}^{g},\mid {W}_{m}^{g}\mid <\epsilon ,$$

(52)

$${\mathbf{F}}^{nk}={\mathbf{F}}^{k}+{\mathbf{W}}^{fk},\mid {W}_{n}^{fk}\mid <\epsilon ,$$

(53)

where *ε* defines the maximum noise magnitude. Recall that = **G** − **A**^{k}**F*** ^{k}*, we have the noise expression of as

$${\overline{\mathbf{G}}}^{n}={\mathbf{W}}^{g}-{\mathbf{A}}^{k}{\mathbf{W}}^{fk}.$$

(54)

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

$${\widehat{\mathbf{F}}}^{un}=\mathbf{V}{\mathbf{\Lambda}}^{-1}{\mathbf{U}}^{T}({\mathbf{W}}^{g}-{\mathbf{A}}^{k}{\mathbf{W}}^{fk}),$$

(55)

which can be further reduced as

$${\widehat{\mathbf{F}}}^{un}={\mathbf{P}}^{un}{\mathbf{W}}^{un},$$

(56)

with **P*** ^{un}* = (

$$\mid {\widehat{F}}_{\overline{n}}^{un}\mid =\mid \sum _{\stackrel{\sim}{m}=1}^{\stackrel{\sim}{M}}{P}_{\overline{n},\stackrel{\sim}{m}}^{un}{W}_{\stackrel{\sim}{m}}^{un}\mid \le \sum _{\stackrel{\sim}{m}=1}^{\stackrel{\sim}{M}}\mid {P}_{\overline{n},\stackrel{\sim}{m}}^{un}\mid \mid {W}_{\stackrel{\sim}{m}}^{un}\mid \le \sum _{\stackrel{\sim}{M}=1}^{\stackrel{\sim}{m}}\mid {P}_{\overline{n},\stackrel{\sim}{m}}^{un}\mid \epsilon ={C}_{\overline{n}}^{Nu}\epsilon ,$$

(57)

where ${C}_{\overline{n}}^{Nu}={\displaystyle \sum _{\stackrel{\sim}{m}=1}^{\stackrel{\sim}{M}}}\mid {P}_{\overline{n},\stackrel{\sim}{m}}^{un}\mid $ and ${C}_{\overline{n}}^{Nu}$ is called the stability measure in this paper. By evaluating the values of ${C}_{\overline{n}}^{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 *c _{sb}* = −1,

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).

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.

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]

$${P}_{\text{o}}(\theta ,s)={e}^{-\underset{-\infty}{\overset{+\infty}{\int}}\mu (s\mathbf{\theta}+\stackrel{\sim}{t}{\mathbf{\theta}}^{\perp})d\stackrel{\sim}{t}}\underset{-\infty}{\overset{\infty}{\int}}f(s\mathbf{\theta}+t{\mathbf{\theta}}^{\perp})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.

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.

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.

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |