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

**|**HHS Author Manuscripts**|**PMC2863332

Formats

Article sections

- Abstract
- I. Introduction
- II. Theoretical Analysis
- III. Algorithm Development
- IV. Discussions and Conclusion
- References

Authors

Related links

Phys Med Biol. Author manuscript; available in PMC 2010 May 4.

Published in final edited form as:

Published online 2009 April 15. doi: 10.1088/0031-9155/54/9/014

PMCID: PMC2863332

NIHMSID: NIHMS192583

Hengyong Yu and Ge Wang

Biomedical Imaging Division, VT-WFU School of Biomedical Engineering and Science, Virginia Tech, Blacksburg, VA 24061, USA

The publisher's final edited version of this article is available at Phys Med Biol

See "Interior region-of-interest reconstruction using a small, nearly piecewise constant subregion." in *Med Phys*, volume 38 on page 1307.

See other articles in PMC that cite the published article.

While the conventional wisdom is that the interior problem does not have a unique solution, by analytic continuation we recently showed that the interior problem can be uniquely and stably solved if we have a known sub-region inside a region-of-interest (ROI). However, such a known sub-region does not always readily available, and it is even impossible to find in some cases. Based on the compressed sensing theory, here we prove that if an object under reconstruction is essentially piecewise constant, a local ROI can be exactly and stably reconstructed via the total variation minimization. Because many objects in CT applications can be approximately modeled as piecewise constant, our approach is practically useful and suggests a new research direction of interior tomography. To illustrate the merits of our finding, we develop an iterative interior reconstruction algorithm that minimizes the total variation of a reconstructed image, and evaluate the performance in numerical simulation.

One conventional wisdom is that the interior problem (exact reconstruction of an interior ROI only from data associated with lines through the ROI) does not have a unique solution (Natterer, 2001). Nevertheless, it is highly desirable to perform interior reconstruction for radiation dose reduction and other important reasons. Hence, over past years a number of image reconstruction algorithms were developed that use an increasingly less amount of raw data (Parker, 1982; Noo *et al*., 2004; Defrise *et al*., 2006). Specifically, motivated by the major needs in cardiac CT, CT guided procedures, nano-CT and so on (Wang *et al*., 2007), by analytic continuation we recently proved that the interior problem can be exactly and stably solved if a sub-region in an ROI within a field-of–view (FOV) is known (Ye *et al*., 2007b; Ye *et al*., 2007a; Ye *et al*., 2008; Yu *et al*., 2008). Similar results were also independently reported by others (Kudo *et al*., 2008; Courdurier *et al*., 2008). Although the CT numbers of certain sub-regions such as air in a trachea and blood in an aorta can be indeed assumed, how to obtain precise knowledge of a sub-region generally can be difficult in some cases such as in studies on rare fossils or certain biomedical structures. Therefore, it would be very valuable to develop more powerful interior tomography techniques (Wang and Yu, 2008).

Another conventional wisdom is that data acquisition should be based on the Nyquist sampling theory, which states that to reconstruct a band-limited signal or image, the sampling rate must at least double the highest frequency of nonzero magnitude. Very interestingly, an alternative theory of compressed sensing (CS) has recently emerged which shows that high-quality signals and images can be reconstructed from far fewer data/measurements than what is usually considered necessary according to the Nyquist sampling theory (Donoho, 2006; Candes *et al*., 2006). The main idea of CS is that most signals are sparse in an appropriate orthonormal system, that is, a majority of their coefficients are close or equal to zero, when represented in some domain. Typically, CS starts with taking a limited amount of samples in a much less correlated basis, and the signal is exactly recovered with an overwhelming probability from the limited data via the _{1} norm minimization. Since samples are limited, the task of recovering the image would involve solving an underdetermined matrix equation, that is, there is a huge amount of candidate images that can all fit the limited measurements effectively. Thus, some additional constraint is needed to select the “best” candidate. While the classical solution to such problems is to minimize the _{2} norm, Donoho and Tao’s groups showed that finding the candidate with the minimum _{1} norm, also equivalent to the total variation (TV) minimization in some cases (Rudin *et al*., 1992), is a most reasonable choice, and can be expressed as a linear program and solved efficiently using existing methods (Donoho, 2006; Candes *et al*., 2006).

Inspired by the CS theory, in this paper we show the feasibility of exact interior tomography in the CS framework. In the next section, we perform a theoretical analysis. In the third section, we develop a CS-based interior tomography algorithm and report numerical simulation results. Finally, we discuss relevant issues.

It is well known that the CS theory depends on the principle of transform sparsity. For an object of interest such as a digital image, we can arrange it as a vector, and in numerous cases there exists an orthonormal basis to make the object sparse in terms of significant transform coefficients. In CS-based image reconstruction, frequently used sparsifying transforms are discrete gradient transforms and wavelet transforms. The discrete gradient sparsifying transform was recently utilized in CT reconstruction (Chen *et al*., 2008). This is because the x-ray attenuation coefficient often varies mildly within organs, and large image variations are usually confined to the borders of tissue structures. A sparse gradient image may also be a good image model in industrial or security applications.

Now, let us analyze the possible exactness of interior tomography in the CS framework subject to the TV minimization. Without loss of generality, let us consider a 2D smooth image *f*(*$\stackrel{\u20d7}{r}$*) = *f*(*x, y*) = *f*(*ρ, θ*), *ρ* [0,1], *θ* [0,2*π*) on the compact support unit disk Ω. Its Radon transform can be written as *R*(*s, *), *s* [−1,1], [0, *π*). Suppose that we are only interested in its interior part *f*(*ρ, θ*) with *ρ* < *ρ*_{0}, and we only know the corresponding local Radon transforms *R*(*s, *), |*s*| < *ρ*_{0}, which is also referred to as local parallel-beam projections. Based on the classic analysis (Natterer, 2001), in general there is no unique solution if we only know these local data. For any reconstructed image from the local data set, it can be viewed as an exact reconstruction from a complete dataset *R*(*s, *), *s* [−1,1], [0, *π*) and a global dataset (*s, *), *ρ*_{0} < |*s*|≤1, [0,*π*). Although (*s, *) = 0 for |*s*|<*ρ*_{0}, it can still produce a non-zero 2D local image *g*(*ρ,θ*), *ρ* [0, *ρ*_{0}), *θ* [0,2*π*) inside the ROI, which is the reason for the non-uniqueness (Natterer, 2001).

For any Radon transform (*s, *), *ρ*_{0} < |*s*|≤1, [0,*π*), the corresponding reconstructed image *g*(*ρ,θ*) with *ρ* < *ρ*_{0} is smooth and bounded if (*s,*) is continuous and bounded.

To prove Lemma 2.1, let us use the classical filtered backprojection (FBP) algorithm. In the filtration step, the filtered data can be written as:

$${\stackrel{\sim}{R}}_{f}(s,\varphi )=\underset{-1}{\overset{1}{\int}}\stackrel{\sim}{R}(t,\varphi )h(s-t)\mathit{dt}=(\underset{-1}{\overset{-{\rho}_{0}}{\int}}+\underset{{\rho}_{0}}{\overset{1}{\int}})\left(\stackrel{\sim}{R}(t,\varphi )h(s-t)\mathit{dt}\right),$$

(1)

where *h*(*s*) represents the ramp filtering kernel. Given the property that *h*(*s*) is analytic everywhere except at *s* = 0, * _{f}*(

$$g(\rho ,\theta )=\underset{0}{\overset{\pi}{\int}}{\stackrel{\sim}{R}}_{f}\phantom{\rule{0.1em}{0ex}}(\rho \phantom{\rule{0.1em}{0ex}}\mathrm{cos}(\varphi -\theta ),\varphi )d\varphi .$$

(2)

Eq. (2) shows that the reconstruction *g*(*ρ,θ*) only involves the data * _{f}*(

For the reconstructed image *g*(*ρ,θ*) with |*ρ*|< *ρ*_{0}, both
$\frac{\partial g(\rho ,\theta )}{\partial \rho}$ and
$\frac{\partial g(\rho ,\theta )}{\rho \partial \theta}$ are smooth and bounded.

Perform a derivation operation on both side of Eq. (2), we have

$$\frac{\partial g(\rho ,\theta )}{\partial \rho}=\underset{0}{\overset{\pi}{\int}}\mathrm{cos}(\varphi -\theta )\frac{\partial {\stackrel{\sim}{R}}_{f}(s,\varphi )}{\partial s}{|}_{s=\rho \mathrm{cos}(\varphi -\theta )}d\varphi ,$$

(3)

and

$$\frac{\partial g(\rho ,\theta )}{\rho \partial \theta}=\underset{0}{\overset{\pi}{\int}}\mathrm{sin}(\varphi -\theta )\frac{\partial {\stackrel{\sim}{R}}_{f}(s,\varphi )}{\partial s}{|}_{s=\rho \mathrm{cos}(\varphi -\theta )}d\varphi .$$

(4)

Using the analytic property of * _{f}*(

Since Natterer has proved that *g*(*ρ,θ*) is a non-zero function, it is impossible that both
$\frac{\partial g(\rho ,\theta )}{\partial \rho}$ and
$\frac{\partial g(\rho ,\theta )}{\rho \partial \theta}$ are zeros inside the ROI. The same results can be extended to higher order derivatives of *g*(*ρ,θ*).

Theoretically speaking, the _{1} norm of the image *f*(*ρ,θ*) inside the ROI can be expressed as:

$${f}_{\mathit{tv}}=\underset{0}{\overset{2\pi}{\int}}d\theta \underset{0}{\overset{{\rho}_{0}}{\int}}\left|\mu (p,\theta )\right|\rho d\rho ,$$

(5)

where *μ*(*ρ,θ*) represents a sparifying transform. For the commonly used gradient transform in the medical imaging field, we have

$$\mu (\rho ,\theta )=\sqrt{{\left(\frac{\partial f(\rho ,\theta )}{\partial \rho}\right)}^{2}+{\left(\frac{\partial f(\rho ,\theta )}{\rho \partial \theta}\right)}^{2}},$$

(6)

which is the gradient magnitude or absolute value of the maximum directional derivation at (*ρ,θ*). If there is no other statement in this paper, we always assume that *μ*(*ρ,θ*) defined by Eq. (6) and (5) represents the total variation. When there exists an artifact image *g*(*ρ,θ*) due to the data truncation, the total variation will be:

$$\begin{array}{ll}{\stackrel{\sim}{f}}_{\mathit{tv}}& =\underset{0}{\overset{2\pi}{\int}}d\theta \underset{0}{\overset{{\rho}_{0}}{\int}}\stackrel{\sim}{\mu}(\rho ,\theta )\rho d\rho \\ & =\underset{0}{\overset{2\pi}{\int}}d\theta \underset{0}{\overset{{\rho}_{0}}{\int}}\rho \sqrt{{\left(\frac{\partial f(\rho ,\theta )}{\partial \rho}+\lambda \frac{\partial g(\rho ,\theta )}{\partial \rho}\right)}^{2}+{\left(\frac{\partial f(\rho ,\theta )}{\rho \partial \theta}+\lambda \frac{\partial g(\rho ,\theta )}{\rho \partial \theta}\right)}^{2}}d\rho \end{array},$$

(7)

where (*ρ,θ*) represents the sparifying transform of a reconstructed image including an artifact image *g*(*ρ,θ*) and *λ* is a coefficient. If we can prove that * _{tv}* can be minimized at

In the CS framework, it is impossible to reconstruct exactly an interior ROI of a general 2D smooth function by minimizing the total variation Eq. (7).

First, let us consider a special circular symmetric case. Without loss of generality, we assume that (*s,*) = (−*s,*) = (*s*). Thus, the corresponding reconstructed image *g*(*ρ,θ*) will be independent of the variable *θ*, that is, *g*(*ρ,θ*) = *g*(*ρ*). Because *g*(*ρ,θ*) is bounded and smooth (Lemma 2.1), we can always construct a circular symmetric local reconstruction *f*(*ρ*) = *g*(*ρ*) + *C* > 0 for *ρ* < *ρ*_{0} with a constant *C*. Note that
$\frac{df(\rho )}{d\rho}=\frac{dg(\rho )}{d\rho}$ for *ρ* < *ρ*_{0}, the TV of *f*(*ρ*) inside the ROI can be computed as:

$${f}_{\mathit{tv}}=\underset{0}{\overset{2\pi}{\int}}d\theta \underset{0}{\overset{{\rho}_{0}}{\int}}\mid \frac{\partial f(\rho ,\theta )}{\partial \rho}\mid \rho d\rho =2\pi \underset{0}{\overset{{\rho}_{0}}{\int}}\mid \frac{df(\rho )}{d\rho}\mid \rho d\rho =2\pi \underset{0}{\overset{{\rho}_{0}}{\int}}\mid \frac{dg(\rho )}{d\rho}\mid \rho d\rho ,$$

(8)

and

$${\stackrel{\sim}{f}}_{\mathit{tv}}=2\pi \underset{0}{\overset{{\rho}_{0}}{\int}}\mid \frac{df(\rho )}{d\rho}+\lambda \frac{dg(\rho )}{d\rho}\mid \rho d\rho =2\pi \left|1+\lambda \right|\underset{0}{\overset{{\rho}_{0}}{\int}}\mid \frac{dg(\rho )}{d\rho}\mid \rho d\rho .$$

(9)

Clearly, * _{tv}* reaches the minimal at

Many biomedical images can be approximately modeled as being piecewise constants; for example in (Wang *et al*., 2004). Now, let us establish the exactness of interior tomography of a symmetric piecewise constant image *f*(*ρ*). In this case, we have

$$\frac{\partial f(\rho ,\theta )}{\partial \rho}=\frac{\partial f(\rho )}{\partial \rho}={f}^{\prime}(\rho )=\sum _{m=1}^{M}{q}_{m}\delta (\rho -{t}_{m}),\rho <{\rho}_{0},$$

(10)

where *t _{m}* > 0 indicates a location of discontinuity,

$${q}_{m}=\underset{\epsilon \to {0}^{+}}{\mathrm{lim}}(f({t}_{m}+\epsilon )-f({t}_{m}-\epsilon )),$$

(11)

and *δ*(•) represents Dirac delta function defined as:

$$\delta (\rho )=\{\begin{array}{ll}+\infty & \rho =0\\ 0& \rho \ne 0\end{array},\underset{-\infty}{\overset{\infty}{\int}}\delta (\rho )d\rho =1.$$

(12)

For two bounded functions *g*_{1}(*ρ*) and *g*_{2}(*ρ*) with 0 < *ρ* < *ρ*_{0}, there will be:

$$\underset{\epsilon \to {0}^{+}}{\mathrm{lim}}\underset{{t}_{m}-\epsilon}{\overset{{t}_{m}+\epsilon}{\int}}\rho \sqrt{{({q}_{m}\delta (\rho -{t}_{m})+{g}_{1}(\rho ))}^{2}+{({g}_{2}(\rho ))}^{2}}d\rho =\mid {q}_{m}\mid {t}_{m},$$

(13)

where *q _{m}* and

Since *g*_{1}(*ρ*) and *g*_{2}(*ρ*) are bounded, there will be
$\underset{\epsilon \to {0}^{+}}{\mathrm{lim}}\underset{{t}_{m}-\epsilon}{\overset{{t}_{m}+\epsilon}{\int}}\mid {g}_{1}(\rho )\mid d\rho =0$ and
$\underset{\epsilon \to {0}^{+}}{\mathrm{lim}}\underset{{t}_{m}-\epsilon}{\overset{{t}_{m}+\epsilon}{\int}}\mid {g}_{2}(\rho )\mid d\rho =0$. Denoting the left side of Eq. (13) as “*LS*”, we have

$$\begin{array}{ll}\mathit{LS}& \le \underset{\epsilon \to {0}^{+}}{\mathrm{lim}}\underset{{t}_{m}-\epsilon}{\overset{{t}_{m}+\epsilon}{\int}}\left(\mid {q}_{m}\delta (\rho -{t}_{m})\mid +\mid {g}_{1}(\rho )\mid +\mid {g}_{2}(\rho )\mid \right)\rho d\rho \\ & \le \underset{\epsilon \to {0}^{+}}{\mathrm{lim}}({t}_{m}+\epsilon )\underset{{t}_{m}-\epsilon}{\overset{{t}_{m}+\epsilon}{\int}}\left(\mid {q}_{m}\delta (\rho -{t}_{m})\mid +\mid {g}_{1}(\rho )\mid +\mid {g}_{2}(\rho )\mid \right)d\rho \\ & ={t}_{m}\underset{\epsilon \to {0}^{+}}{\mathrm{lim}}\underset{{t}_{m}-\epsilon}{\overset{{t}_{m}+\epsilon}{\int}}\left(\mid {q}_{m}\delta (\rho -{t}_{m})\mid +\mid {g}_{1}(\rho )\mid +\mid {g}_{2}(\rho )\mid \right)d\rho \\ & =\mid {q}_{m}\mid {t}_{m}\underset{\epsilon \to {0}^{+}}{\mathrm{lim}}\underset{{t}_{m}-\epsilon}{\overset{{t}_{m}+\epsilon}{\int}}\delta (\rho -{t}_{m})d\rho +{t}_{m}\underset{\epsilon \to {0}^{+}}{\mathrm{lim}}\underset{{t}_{m}-\epsilon}{\overset{{t}_{m}+\epsilon}{\int}}\mid {g}_{1}(\rho )\mid d\rho +{t}_{m}\underset{\epsilon \to {0}^{+}}{\mathrm{lim}}\underset{{t}_{m}-\epsilon}{\overset{{t}_{m}+\epsilon}{\int}}\mid {g}_{2}(\rho )\mid d\rho \\ & =\mid {q}_{m}\mid {t}_{m}\end{array}.$$

(14)

On the other hand, we have

$$\begin{array}{ll}\mathit{LS}& \ge \underset{\epsilon \to {0}^{+}}{\mathrm{lim}}\underset{{t}_{m}-\epsilon}{\overset{{t}_{m}+\epsilon}{\int}}\left(\mid {q}_{m}\delta (\rho -{t}_{m})\mid -\mid {g}_{1}(\rho )\mid -\mid {g}_{2}(\rho )\mid \right)\rho d\rho \\ & \ge \underset{\epsilon \to {0}^{+}}{\mathrm{lim}}({t}_{m}-\epsilon )\underset{{t}_{m}-\epsilon}{\overset{{t}_{m}+\epsilon}{\int}}\left(\mid {q}_{m}\delta (\rho -{t}_{m})\mid -\mid {g}_{1}(\rho )\mid -\mid {g}_{2}(\rho )\mid \right)d\rho \\ & ={t}_{m}\underset{\epsilon \to {0}^{+}}{\mathrm{lim}}\underset{{t}_{m}-\epsilon}{\overset{{t}_{m}+\epsilon}{\int}}\left(\mid {q}_{m}\delta (\rho -{t}_{m})\mid -\mid {g}_{1}(\rho )\mid -\mid {g}_{2}(\rho )\mid \right)d\rho \\ & =\mid {q}_{m}\mid {t}_{m}\underset{\epsilon \to {0}^{+}}{\mathrm{lim}}\underset{{t}_{m}-\epsilon}{\overset{{t}_{m}+\epsilon}{\int}}\delta (\rho -{t}_{m})d\rho -{t}_{m}\underset{\epsilon \to {0}^{+}}{\mathrm{lim}}\underset{{t}_{m}-\epsilon}{\overset{{t}_{m}+\epsilon}{\int}}\mid {g}_{1}(\rho )\mid d\rho -{t}_{m}\underset{\epsilon \to {0}^{+}}{\mathrm{lim}}\underset{{t}_{m}-\epsilon}{\overset{{t}_{m}+\epsilon}{\int}}\mid {g}_{2}(\rho )\mid d\rho \\ & =\mid {q}_{m}\mid {t}_{m}\end{array}.$$

(15)

Combining (14) and (15), we have

$$\mid {q}_{m}\mid {t}_{m}\le \mathit{LS}\le \mid {q}_{m}\mid {t}_{m},$$

(16)

which implies that *LS* = |*q _{m}*|

In the CS framework, an interior ROI of a circular symmetric piecewise constant function *f*(*ρ*) can be exactly determined by minimizing the total variation defined in Eq. (7).

For a symmetric and piecewise constant function *f*(*ρ*), the TV formula Eq. (5) becomes

$${f}_{\mathit{tv}}=2\pi \underset{0}{\overset{{\rho}_{0}}{\int}}\left|{f}^{\prime}(\rho )\right|\rho d\rho =2\pi \sum _{m=1}^{M}\mid {q}_{m}\mid {t}_{m}.$$

(17)

Denoting the set of finite discontinuous points of *f*(*ρ*) at *t _{m}* as

$${\stackrel{\sim}{f}}_{\mathit{tv}}=\left|\lambda \right|{\stackrel{\sim}{f}}_{\mathit{tvg}}+{\stackrel{\sim}{f}}_{\mathit{tvf}},$$

(18)

where

$${\stackrel{\sim}{f}}_{\mathit{tvg}}=\underset{0}{\overset{2\pi}{\int}}d\theta \underset{I-{I}_{D}}{\int}\rho \sqrt{{\left(\frac{\partial g(\rho ,\theta )}{\partial \rho}\right)}^{2}+{\left(\frac{\partial g(\rho ,\theta )}{\rho \partial \theta}\right)}^{2}}d\rho ,$$

(19)

and

$${\stackrel{\sim}{f}}_{\mathit{tvf}}=\underset{0}{\overset{2\pi}{\int}}d\theta \underset{{I}_{D}}{\int}\rho \sqrt{{\left(\sum _{m=1}^{M}{q}_{m}\delta (\rho -{t}_{m})+\lambda \frac{\partial g(\rho ,\theta )}{\partial \rho}\right)}^{2}+{\left(\lambda \frac{\partial g(\rho ,\theta )}{\rho \partial \theta}\right)}^{2}}d\rho .$$

(20)

Eq. (20) can be simplified as

$$\begin{array}{ll}{\stackrel{\sim}{f}}_{\mathit{tvf}}& =\underset{0}{\overset{2\pi}{\int}}d\theta \sum _{m=1}^{M}\underset{\epsilon \to {0}^{+}}{\mathrm{lim}}\underset{{t}_{m}-\epsilon}{\overset{{t}_{m}+\epsilon}{\int}}\rho \sqrt{{\left({q}_{m}\delta (\rho -{t}_{m})+\lambda \frac{\partial g(\rho ,\theta )}{\partial \rho}\right)}^{2}+{\left(\lambda \frac{\partial g(\rho ,\theta )}{\rho \partial \theta}\right)}^{2}}d\rho \\ & =\underset{0}{\overset{2\pi}{\int}}d\theta \sum _{m=1}^{M}\mid {q}_{m}\mid {t}_{m}=2\pi \sum _{m=1}^{M}\mid {q}_{m}\mid {t}_{m}={f}_{\mathit{tv}}\end{array},$$

(21)

where Lemma 2.3 has been used. Clearly, Eq. (18) reaches the minimal at *λ* = 0, which implies that exact reconstruction can be performed in terms of TV minimization only from purely local data.

Next, let us consider a general piecewise constant function *f*(*ρ,θ*) defined on the compact supported unit disk Ω. Without loss of generality, we assume that Ω can be divided into finite sub-regions Ω_{n}, *n* = 1,2,…, *N*, where each sub-region Ω*n* has a nonzero area measure, on which *f*(*ρ,θ*) is a constant. As a result, these sub-regions also define finitely many boundaries in terms of arc-segments, each of which is of a non-zero length and differentiable almost everywhere excluding at most finitely many points. We assume that an interior ROI Ω_{0} (defined by *ρ* < *ρ*_{0}) covers *M* boundaries and these line segments are denoted as *L*_{1}, *L*_{2},…*L _{m}*…,

$$\mu (\overrightarrow{r}+v\overrightarrow{v})={q}_{m}\delta (v).$$

(22)

On the other hand, we have a global polar system with the unit vectors being and . If we denote the angle between and *$\stackrel{\u20d7}{v}$* as *γ*, the gradient of *f*(*$\stackrel{\u20d7}{r}$*) along *$\stackrel{\u20d7}{v}$* in a small neighborhood of *$\stackrel{\u20d7}{r}$* can be written as *q _{m}δ*(

$$\stackrel{\sim}{\mu}(\overrightarrow{r}+v\overrightarrow{v})=\sqrt{{\left({q}_{m}\delta (v)\mathrm{cos}\gamma +\lambda \frac{\partial g(\overrightarrow{r}+v\overrightarrow{v})}{\partial \rho}\right)}^{2}+{\left({q}_{m}\delta (v)\mathrm{sin}\gamma +\lambda \frac{\partial g(\overrightarrow{r}+v\overrightarrow{v})}{\rho \partial \theta}\right)}^{2}}.$$

(23)

For two bounded functions *g*_{1}(*$\stackrel{\u20d7}{r}$*) and *g*_{2}(*$\stackrel{\u20d7}{r}$*), there will be:

$$\underset{\epsilon \to {0}^{+}}{\mathrm{lim}}\underset{-\epsilon}{\overset{\epsilon}{\int}}\sqrt{{\left({q}_{m}\delta (v)\mathrm{cos}\gamma +\lambda {g}_{1}(\overrightarrow{r}+v\overrightarrow{v})\right)}^{2}+{\left({q}_{m}\delta (v)\mathrm{sin}\gamma +\lambda {g}_{2}(\overrightarrow{r}+v\overrightarrow{v})\right)}^{2}}\mathit{dv}=\mid {q}_{m}\mid .$$

(24)

Denoting the left side of Eq. (24) as “*LS*”, we have

$$\begin{array}{lll}\mathit{LS}& =& \underset{\epsilon \to {0}^{+}}{\mathrm{lim}}\underset{-\epsilon}{\overset{\epsilon}{\int}}\sqrt{{\left({q}_{m}\delta (v)\right)}^{2}+{\left(\lambda {g}_{1}(\overrightarrow{r}+v\overrightarrow{v})\right)}^{2}+{\left(\lambda {g}_{2}(\overrightarrow{r}+v\overrightarrow{v})\right)}^{2}+2\lambda {q}_{m}\delta (v)\left(\mathrm{cos}\phantom{\rule{0.2em}{0ex}}\gamma \phantom{\rule{0.2em}{0ex}}{g}_{1}(\overrightarrow{r}+v\overrightarrow{v})+\mathrm{sin}\phantom{\rule{0.2em}{0ex}}\gamma \phantom{\rule{0.2em}{0ex}}{g}_{2}(\overrightarrow{r}+v\overrightarrow{v})\right)}\mathit{dv}\\ & \le & \underset{\epsilon \to {0}^{+}}{\mathrm{lim}}\underset{-\epsilon}{\overset{\epsilon}{\int}}\sqrt{\mid {q}_{m}\delta (v)\mid \left(\mid 2\lambda \left(\mathrm{cos}\phantom{\rule{0.2em}{0ex}}\gamma \phantom{\rule{0.2em}{0ex}}{g}_{1}(\overrightarrow{r}+v\overrightarrow{v})+\mathrm{sin}\phantom{\rule{0.2em}{0ex}}\gamma \phantom{\rule{0.2em}{0ex}}{g}_{2}(\overrightarrow{r}+v\overrightarrow{v})\right)\mid +\mid {q}_{m}\delta (v)\mid \right)+{\left(\lambda {g}_{1}(\overrightarrow{r}+v\overrightarrow{v})\right)}^{2}+{\left(\lambda {g}_{2}(\overrightarrow{r}+v\overrightarrow{v})\right)}^{2}}\mathit{dv}\\ & \le & \underset{\epsilon \to {0}^{+}}{\mathrm{lim}}\underset{-\epsilon}{\overset{\epsilon}{\int}}\sqrt{{\left(\mid 2\lambda (\mathrm{cos}\phantom{\rule{0.2em}{0ex}}\gamma \phantom{\rule{0.2em}{0ex}}{g}_{1}(\overrightarrow{r}+v\overrightarrow{v})+\mathrm{sin}\phantom{\rule{0.2em}{0ex}}\gamma \phantom{\rule{0.2em}{0ex}}{g}_{2}(\overrightarrow{r}+v\overrightarrow{v}))\mid +\mid {q}_{m}\delta (v)\mid \right)}^{2}+{\left(\lambda {g}_{1}(\overrightarrow{r}+v\overrightarrow{v})\right)}^{2}+{\left(\lambda {g}_{2}(\overrightarrow{r}+v\overrightarrow{v})\right)}^{2}}\mathit{dv}\\ & \le & \underset{\epsilon \to {0}^{+}}{\mathrm{lim}}\underset{-\epsilon}{\overset{\epsilon}{\int}}\left(\mid 2\lambda \left(\mathrm{cos}\phantom{\rule{0.2em}{0ex}}\gamma \phantom{\rule{0.2em}{0ex}}{g}_{1}(\overrightarrow{r}+v\overrightarrow{v})+\mathrm{sin}\phantom{\rule{0.2em}{0ex}}\gamma \phantom{\rule{0.2em}{0ex}}{g}_{2}(\overrightarrow{r}+v\overrightarrow{v})\right)\mid +\mid {q}_{m}\delta (v)\mid +\mid \lambda {g}_{1}(\overrightarrow{r}+v\overrightarrow{v})\mid +\mid \lambda {g}_{2}(\overrightarrow{r}+v\overrightarrow{v})\mid \right)\mathit{dv}\\ & =& \underset{\epsilon \to {0}^{+}}{\mathrm{lim}}\underset{-\epsilon}{\overset{\epsilon}{\int}}\mid {q}_{m}\delta (v)\mid \mathit{dv}+\underset{\epsilon \to {0}^{+}}{\mathrm{lim}}\underset{-\epsilon}{\overset{\epsilon}{\int}}\mid 2\lambda \left(\mathrm{cos}\phantom{\rule{0.2em}{0ex}}\gamma \phantom{\rule{0.2em}{0ex}}{g}_{1}(\overrightarrow{r}+v\overrightarrow{v})+\mathrm{sin}\phantom{\rule{0.2em}{0ex}}\gamma \phantom{\rule{0.2em}{0ex}}{g}_{2}(\overrightarrow{r}+v\overrightarrow{v})\right)\mid \mathit{dv}\\ & +& \underset{\epsilon \to {0}^{+}}{\mathrm{lim}}\underset{-\epsilon}{\overset{\epsilon}{\int}}\mid \lambda {g}_{1}(\overrightarrow{r}+v\overrightarrow{v})\mid \mathit{dv}+\underset{\epsilon \to {0}^{+}}{\mathrm{lim}}\underset{-\epsilon}{\overset{\epsilon}{\int}}\mid \lambda {g}_{2}(\overrightarrow{r}+v\overrightarrow{v})\mid \mathit{dv}\\ & =& \mid {q}_{m}\mid +0+0+0=\mid {q}_{m}\mid ,\end{array}$$

(25)

where we have used the facts that *λ*, *g*_{1}(*$\stackrel{\u20d7}{r}$*) and *g*_{2}(*$\stackrel{\u20d7}{r}$*) are bounded. On the other hand,

$$\begin{array}{lll}\mathit{LS}& =& \underset{\epsilon \to {0}^{+}}{\mathrm{lim}}\underset{-\epsilon}{\overset{\epsilon}{\int}}\sqrt{{\left({q}_{m}\delta (v)\right)}^{2}+{\left(\lambda {g}_{1}(\overrightarrow{r}+v\overrightarrow{v})\right)}^{2}+{\left(\lambda {g}_{2}(\overrightarrow{r}+v\overrightarrow{v})\right)}^{2}+2\lambda {q}_{m}\delta (v)\left(\mathrm{cos}\phantom{\rule{0.2em}{0ex}}\gamma \phantom{\rule{0.2em}{0ex}}{g}_{1}(\overrightarrow{r}+v\overrightarrow{v})+\mathrm{sin}\phantom{\rule{0.2em}{0ex}}\gamma \phantom{\rule{0.2em}{0ex}}{g}_{2}(\overrightarrow{r}+v\overrightarrow{v})\right)}\mathit{dv}\\ & \ge & \underset{\epsilon \to {0}^{+}}{\mathrm{lim}}\underset{-\epsilon}{\overset{\epsilon}{\int}}\sqrt{\mid {q}_{m}\delta (v)\mid (\mid {q}_{m}\delta (v)\mid -\mid 2\lambda (\mathrm{cos}\phantom{\rule{0.2em}{0ex}}\gamma \phantom{\rule{0.2em}{0ex}}{g}_{1}(\overrightarrow{r}+v\overrightarrow{v})+\mathrm{sin}\phantom{\rule{0.2em}{0ex}}\gamma \phantom{\rule{0.2em}{0ex}}{g}_{2}(\overrightarrow{r}+v\overrightarrow{v}))\mid )+{\left(\lambda {g}_{1}(\overrightarrow{r}+v\overrightarrow{v})\right)}^{2}+{\left(\lambda {g}_{2}(\overrightarrow{r}+v\overrightarrow{v})\right)}^{2}}\mathit{dv}\\ & \ge & \underset{\epsilon \to {0}^{+}}{\mathrm{lim}}\underset{-\epsilon}{\overset{\epsilon}{\int}}\sqrt{{\left(\mid {q}_{m}\delta (v)\mid -\mid 2\lambda (\mathrm{cos}\phantom{\rule{0.2em}{0ex}}\gamma \phantom{\rule{0.2em}{0ex}}{g}_{1}(\overrightarrow{r}+v\overrightarrow{v})+\mathrm{sin}\phantom{\rule{0.2em}{0ex}}\gamma \phantom{\rule{0.2em}{0ex}}{g}_{2}(\overrightarrow{r}+v\overrightarrow{v}))\mid \right)}^{2}+{\left(\lambda {g}_{1}(\overrightarrow{r}+v\overrightarrow{v})\right)}^{2}+{\left(\lambda {g}_{2}(\overrightarrow{r}+v\overrightarrow{v})\right)}^{2}}\mathit{dv}\\ & \ge & \underset{\epsilon \to {0}^{+}}{\mathrm{lim}}\underset{-\epsilon}{\overset{\epsilon}{\int}}\left(\mid {q}_{m}\delta (v)\mid -\mid 2\lambda \left(\mathrm{cos}\phantom{\rule{0.2em}{0ex}}\gamma \phantom{\rule{0.2em}{0ex}}{g}_{1}(\overrightarrow{r}+v\overrightarrow{v})+\mathrm{sin}\phantom{\rule{0.2em}{0ex}}\gamma \phantom{\rule{0.2em}{0ex}}{g}_{2}(\overrightarrow{r}+v\overrightarrow{v})\right)\mid -\sqrt{{\left(\lambda {g}_{1}(\overrightarrow{r}+v\overrightarrow{v})\right)}^{2}+{\left(\lambda {g}_{2}(\overrightarrow{r}+v\overrightarrow{v})\right)}^{2}}\right)\mathit{dv}\\ & =& \underset{\epsilon \to {0}^{+}}{\mathrm{lim}}\underset{-\epsilon}{\overset{\epsilon}{\int}}\mid {q}_{m}\delta (v)\mid \mathit{dv}-\underset{\epsilon \to {0}^{+}}{\mathrm{lim}}\underset{-\epsilon}{\overset{\epsilon}{\int}}\mid 2\lambda \left(\mathrm{cos}\phantom{\rule{0.2em}{0ex}}\gamma \phantom{\rule{0.2em}{0ex}}{g}_{1}(\overrightarrow{r}+v\overrightarrow{v})+\mathrm{sin}\phantom{\rule{0.2em}{0ex}}\gamma \phantom{\rule{0.2em}{0ex}}{g}_{2}(\overrightarrow{r}+v\overrightarrow{v})\right)\mid \mathit{dv}\\ & -& \underset{\epsilon \to {0}^{+}}{\mathrm{lim}}\underset{-\epsilon}{\overset{\epsilon}{\int}}\sqrt{{\left(\lambda {g}_{1}(\overrightarrow{r}+v\overrightarrow{v})\right)}^{2}+{\left(\lambda {g}_{2}(\overrightarrow{r}+v\overrightarrow{v})\right)}^{2}}\mathit{dv}\\ & =& \mid {q}_{m}\mid -0-0=\mid {q}_{m}\mid .\end{array}$$

(26)

In the CS framework, an interior ROI of a general compactly supported function *f*(*$\stackrel{\u20d7}{r}$*) can be exactly determined by minimizing the total variation defined by Eq.(7) if *f*(*$\stackrel{\u20d7}{r}$*) can be decomposed into finitely many constant sub-regions.

For an ROI Ω_{0} of a general piecewise constant function *f*(*$\stackrel{\u20d7}{r}$*), the TV formula Eq. (5) becomes

$${f}_{\mathit{tv}}=\underset{{\mathrm{\Omega}}_{0}}{\iint}\mu (\overrightarrow{r})d\overrightarrow{r}=\sum _{m=1}^{M}\underset{{L}_{m}}{\int}\mu (\overrightarrow{r})d\overrightarrow{r}.$$

(27)

To compute
$\underset{{L}_{m}}{\int}\mu (\overrightarrow{r})d\overrightarrow{r}$, we consider a curved strip Ω_{Lm} the central axis *L _{m}* and width 2

$$\underset{{L}_{m}}{\int}\mu (\overrightarrow{r})d\overrightarrow{r}=\underset{\epsilon \to {0}^{+}}{\mathrm{lim}}\underset{{\mathrm{\Omega}}_{{L}_{m}}}{\iint}\mu (\overrightarrow{r})d\overrightarrow{r}={t}_{m}\underset{\epsilon \to {0}^{+}}{\mathrm{lim}}\underset{-\epsilon}{\overset{\epsilon}{\int}}\mid {q}_{m}\delta (v)\mid dv={t}_{m}\mid {q}_{m}\mid ,$$

(28)

and

$${f}_{\mathit{tv}}=\underset{{\mathrm{\Omega}}_{0}}{\iint}\mu (\overrightarrow{r})d\overrightarrow{r}=\sum _{m=1}^{M}\mid {q}_{m}\mid {t}_{m},$$

(29)

where the finite many undifferentiable points do not affect the integral result. After the introduction of *g*(*$\stackrel{\u20d7}{r}$*) as discussed in Lemmas 2.1 and 2.2, the TV expression Eq. (7) becomes

$${\stackrel{\sim}{f}}_{\mathit{tv}}=\mid \lambda \mid {\stackrel{\sim}{f}}_{\mathit{tvg}}+{\stackrel{\sim}{f}}_{\mathit{tvf}},$$

(30)

where

$${\stackrel{\sim}{f}}_{\mathit{tvg}}=\underset{{\mathrm{\Omega}}_{0}-\mathrm{\Gamma}}{\iint}\sqrt{{\left(\frac{\partial g(\overrightarrow{r})}{\partial \rho}\right)}^{2}+{\left(\frac{\partial g(\overrightarrow{r})}{\rho \partial \theta}\right)}^{2}}d\overrightarrow{r},\phantom{\rule{0.2em}{0ex}}\mathrm{\Gamma}=\underset{m=1}{\overset{M}{\cup}}{L}_{m},$$

(31)

and

$${\stackrel{\sim}{f}}_{\mathit{tvf}}=\sum _{m=1}^{M}\underset{{L}_{m}}{\int}\stackrel{\sim}{\mu}(\overrightarrow{r})d\overrightarrow{r}.$$

(32)

Using the same strategy for computing $\underset{{L}_{m}}{\int}\mu (\overrightarrow{r})d\overrightarrow{r}$ and Lemma 2.4, we immediately obtain that

$${\stackrel{\sim}{f}}_{\mathit{tvf}}=\sum _{m=1}^{M}\mid {q}_{m}\mid {t}_{m}={f}_{\mathit{tv}}.$$

(33)

Clearly, Eq. (30) reaches the minimal at *λ* = 0, which implies Theorem 2.3.

Two comments are in order on the above theorems. First, the above results are independent of scanning geometries and coordinate systems although we have assumed the parallel-beam geometry and polar coordinate system. Second, if we change the parameter *λ* to a smooth and bounded function *λ*(*ρ,θ*), the same results can be proved.

To verify the theoretical results in Section II, we developed a numerical interior tomography algorithm in the CS framework. The algorithm consists of two major steps. In the first step, the ordered-subset simultaneous algebraic reconstruction technique (OS-SART) (Wang and Jiang, 2004) was used to reconstruct a digital image *f _{m,n}* =

*α*:= 0.005,*α*: 0.997,_{s}*P*: 5;_{TV}*k*: 0, ${f}_{m,n}^{0}$,*P*:= 20 ;_{ART}- Repeat the main loop (OS-SART Reconstruction and Minimizing TV)
*k*:=*k*+ 1; ${f}_{m,n}^{k}:={f}_{m,n}^{k-1}$- For
*p*= 1 to_{art}*P*do_{ART} - Update ${f}_{m,n}^{k}$ by OS-SART using the projections in the
*p*subset;_{art} - For
*p*= 1 to_{tv}*P*do_{TV} - Computing the steepest decent direction
*d*;_{m n} - $\beta :=\mathrm{max}(\mid {f}_{m,n}^{k}\mid )\xf7\mathrm{max}(\mid {d}_{m,n}\mid );$
- ${f}_{m,n}^{k}={f}_{m,n}^{k}-\alpha \times \beta \times {d}_{m,n};$
*α*=*α*×*α*;_{s}- End for (
*p*)_{tv} - End for (
*p*)_{art} - Until the stopping criteria are satisfied.

Because the OS-SART technique was discussed in our previous paper (Wang and Jiang, 2004), here we only explain the second step in the above pseudo-code. Line 1 initializes the control parameters *α, α _{s}* and

$${\mu}_{m,n}\cong \sqrt{\frac{{({f}_{m+1,n}-{f}_{m,n})}^{2}+{({f}_{m,n}-{f}_{m-1,n})}^{2}+{({f}_{m,n+1}-{f}_{m,n})}^{2}+{({f}_{m,n}-{f}_{m,n-1})}^{2}}{2{\nabla}^{2}}}.$$

(34)

Correspondingly, the total variation can be defined as ${f}_{\mathit{TV}}=\sum _{m}\sum _{n}{\mu}_{m,n}$. Then, we have the steepest descent direction defined by

$$\begin{array}{ll}{d}_{m,n}& =\frac{\partial {f}_{\mathit{TV}}}{\partial {f}_{m,n}}=\frac{4{f}_{m,n}-{f}_{m+1,n}-{f}_{m-1,n}-{f}_{m,n+1}-{f}_{m,n-1}}{{\mu}_{m,n}}\\ & +\frac{{f}_{m,n}-{f}_{m+1,n}}{{\mu}_{m+1,n}}+\frac{{f}_{m,n}-{f}_{m-1,n}}{{\mu}_{m-1,n}}+\frac{{f}_{m,n}-{f}_{m,n+1}}{{\mu}_{m,n+1}}+\frac{{f}_{m,n}-{f}_{m,n-1}}{{\mu}_{m,n-1}}\end{array}.$$

(35)

For another sparsifying transform, we can deduce the corresponding *d _{m,n}*, easily. Line 9 computes the normalized factor

Our algorithm was numerically implemented in MatLab on a PC (4.0 Gigabyte memory, 3.16G Hz CPU). While the basic structure was constructed in MatLab, all the computationally intensive parts were coded in C, which was linked via a MEX interface. A maximal iteration time was set to stop the main loop. Because many biomedical images can be approximately modeled as piecewise constants, we employed the commonly used discrete gradient transform in the numerical simulation and animal experiment. To avoid the singularity when computing *d _{m,n}*, using Eq. (35), we added a small constant to Eq.(34) when computing the gradient

$${\mu}_{m,n}\cong \sqrt{\frac{{({f}_{m+1,n}-{f}_{m,n})}^{2}+{({f}_{m,n}-{f}_{m-1,n})}^{2}+{({f}_{m,n+1}-{f}_{m,n})}^{2}+{({f}_{m,n}-{f}_{m,n-1})}^{2}}{2{\nabla}^{2}}+{\epsilon}^{2}}.$$

(36)

In our numerical simulation, we assumed a circular scanning locus of radius 57.0 cm and a fan-beam imaging geometry. We also assumed an equi-spatial virtual detector array of length 12.0 cm. The detector was centered at the system origin and made always perpendicular to the direction from the system origin to the x-ray source. The detector array included 360 elements, each of which is of aperture 0.033 cm. This scanning configuration covered a circular FOV of radius 5.967 cm. For a complete scanning turn, we equi-angularly collected 1300 projections. The reconstructed object was a 2D modified Shepp-logan phantom. This phantom is piecewise constant and includes a set of smooth ellipses whose parameters are listed in Table 1, where *a,b*, represent the *x, y* semi-axes, (*x*_{0}, *y*_{0}) the center of the ellipse, *ω* denotes the rotation angle, *f* the relative attenuation coefficient. The units for *a b*, and (*x*_{0}, *y*_{0}) are cm. The reconstructed images were in a 256×256 matrix covering an FOV of radius 10 cm. The 60 iterations took 60 minutes. For comparison, we also reconstructed an image using a local FBP method with smooth extrapolation from the truncated projections into missing data. Some typical reconstructed images were shown in Figure 4. The typical profiles were in Figure 5. As seen from Figures 4 and and5,5, the reconstructed images from the proposed algorithm are in a high precision inside the ROI.

Reconstructed images of a modified Shepp-logan phantom after 60 iterations. (a) The original phantom, (b) the reconstruction using the local FBP (after smooth data extrapolation), and (c) the reconstruction using the proposed CS-based interior tomography **...**

To demonstrate the real-world application of the proposed algorithm, we performed a CT scan of a living sheep, which was approved by Virginia Tech IACUC committee. The chest of a sheep was scanned in fan-beam geometry on a SIEMENS 64-Slice CT scanner (100kVp, 150mAs). The x-ray source trajectory of radius 57.0 cm was used. There were 1160 projections uniformly collected over a 360° range, and 672 detectors were equi-angularly distributed per projection. Thus, the FOV of radius 25.05 cm was formed. First, an entire 29.06 cm by 29.06 cm cross-section was reconstructed into 1024 ×1024 pixels using the popular FBP method from a complete dataset of projections. Second, a trachea was selected in reference to the reconstructed image. Around the trachea, a circular ROI of radius 120 pixels was specified. Then, only the projection data through the ROI were kept to simulate an interior scan. Third, the ROI was reconstructed by the local FBP with smooth data extrapolation and our proposed algorithms, respectively. The results were in Figures 6. Comparing the images in Figure 6, we observe that the proposed algorithm not only keeps image accuracy but also suppresses image noise.

There are still some differences between the phantom and reconstructed images (see Fig. 5 (b)). These artifacts in the CS-based reconstruction may be due to the following two sources. First, from the viewpoint of signal processing, the sampling process involved in the digitations can be interpreted as a truncation function in the frequency domain. The reconstructed images will not satisfy the piecewise constant condition due to the high-frequency truncation. Second, the reconstruction results will converge to the real values only after an infinite number of iterations. These two problems are also our opportunities to refine CS-based methods in the future. Since our paper is the first of its kind in terms of exact interior tomography, our numerical implementation of the iterative reconstruction algorithm is relatively preliminary. Neither the codes nor the control parameters have been optimized, needless to say the stopping criteria. Thus, we plan to perform optimizations over these issues in follow-up papers. Theoretically speaking, the gradient magnitude of a piecewise constant function will be infinite on the boundaries of different sub-regions. Using the Lambda tomography technique (Yu and Wang, 2006), we can identify these boundaries and handle them differently. Very likely, this type of special processing will give us superior reconstruction performance.

In this paper, we have established the feasibility of exact interior tomography via the total variation minimization in the CS framework. Our theoretical analysis has shown that exact ROI reconstruction is possible for a piecewise constant function. All the analyses in this paper have assumed that the image can be modeled in an appropriate functional form, which is a kind of prior information different from that in our previous papers (Ye *et al*., 2007b; Ye *et al*., 2008; Ye *et al*., 2007a; Yu *et al*., 2008). In the near future, we plan to introduce more prior information into the CS framework, such as low resolution image, a partially known image inside the ROI, the average intensity and standard deviation, *etc*. More generally, our overall intuition is that when we have strong knowledge on an interior ROI so that there is explicitly or implicitly a sparse representation of the true image on the ROI, an exact and stable interior reconstruction of the ROI should be achieved definitely or with a high probability. The key is always to define an appropriate sparsifying transform and an associated objective function. Then, among many solutions to the interior problem, the true one (or its good approximation) can be identified. For example, an image on the ROI is often compressible using the wavelet transform. That is, the desired interior reconstruction indeed has a sparse presentation, at least approximately. In this setting, it would be very interesting to develop a CS approach for interior tomography without any need for precise knowledge on a subregion in the ROI. We are actively working on this type of topics.

Furthermore, in the same spirit of the proofs in Section II by defining a different sparisfying transform *μ*, our results can be extended to some interesting families of functions. Particularly, we have the following conjectures 4.1, 4.2 and 4.3.

Consider a circular symmetric *f*(*ρ*) ≥ 0 defined on the support [0,1]. Suppose that [0,1] can be divided into finitely many sub-intervals on each of which *f*(*ρ*) can be expressed as:

$$f(\rho )=\sum _{k=0}^{K}{C}_{k}{\rho}^{k},K\ge 0,$$

(37)

where *C _{k}* are constants. In the CS framework, an interior ROI can be exactly determined by minimizing the

$$\mu (\rho ,\theta )=\sqrt{{\left(\frac{{\partial}^{K+1}f(\rho ,\theta )}{\partial {\rho}^{K+1}}\right)}^{2}+{\left(\frac{{\partial}^{K+1}f(\rho ,\theta )}{{\rho}^{K+1}\partial {\theta}^{K+1}}\right)}^{2}}.$$

(38)

Consider a general function *f*(*x, y*) ≥ 0 defined on the compact support Ω. Suppose that Ω can be divided into finitely many sub-regions on each of which *f*(*x, y*) can be expressed as:

$$f(x,y)=\sum _{k=0}^{K}({C}_{k}^{x}{x}^{k}+{C}_{k}^{y}{y}^{k}),K\ge 0,$$

(39)

where
${C}_{k}^{x}$ and
${C}_{k}^{y}$ are constants. In the CS framework, an interior ROI can be exactly determined by minimizing the _{1} norm if we define a sparsifying transform as

$$\mu (x,y)=\sqrt{{\left(\frac{{\partial}^{K+1}f(x,y)}{\partial {x}^{K+1}}\right)}^{2}+{\left(\frac{{\partial}^{K+1}f(x,y)}{\partial {y}^{K+1}}\right)}^{2}}.$$

(40)

Consider a general function *f*(*x, y*) ≥ 0 defined on the compact support Ω. Suppose that Ω can be divided into finitely many sub-regions on each of which *f*(*x, y*) has a constant Gaussian curvature *κ*(*x, y*) defined as:

$$\kappa (x,y)=\frac{\frac{{\partial}^{2}f(x,y)}{\partial {x}^{2}}\frac{{\partial}^{2}f(x,y)}{\partial {y}^{2}}-{\left(\frac{{\partial}^{2}f(x,y)}{\partial x\partial y}\right)}^{2}}{{\left(1+{\left(\frac{\partial f(x,y)}{\partial x}\right)}^{2}+{\left(\frac{\partial f(x,y)}{\partial y}\right)}^{2}\right)}^{2}}.$$

(41)

In the CS framework, an interior ROI can be exactly determined by minimizing the _{1} norm if we define a sparsifying transform as

$$\mu (x,y)=\sqrt{{\left(\frac{\partial \kappa (x,y)}{\partial x}\right)}^{2}+{\left(\frac{\partial \kappa (x,y)}{\partial y}\right)}^{2}}.$$

(42)

It should be pointed out that Eqs. (38) and (40) are the same thing when *K* = 0 because they represent the gradient magnitude. Eqs. (41) and (42) are independent of the coordinate system. While Conjecture 4.1 can be proved similar to that of Theorem 2.2, Conjectures 4.2 and 4.3 can be proved similar to that of Theorem 2.3. The major difference lies in that the operations of the higher order derivatives of Dirac delta function must be involved.

In conclusion, we have performed a theoretical analysis for interior tomography in the CS framework, and showed for the first time that exact knowledge on a sub-region of an ROI is not necessary for exact and stable interior reconstruction. Our numerical and experimental data have validated our theoretical analysis and algorithmic implementation, and demonstrated that our approach may have a great potential in medical, industrial and other applications. Furthermore, our work suggests a new direction of interior tomography, and more efforts and results are warranted.

This work is partially supported by NIH/NIBIB grants (EB002667, EB004287, EB007288).

- Candes EJ, Romberg J, Tao T. Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information. IEEE Transactions on Information Theory. 2006;52:489–509.
- Chen GH, Tang J, Leng S. Prior image constrained compressed sensing (PICCS): A method to accurately reconstruct dynamic CT images from highly undersampled projection data sets. Medical Physics. 2008;35:660–3. [PMC free article] [PubMed]
- Courdurier M, Noo F, Defrise M, Kudo H. Solving the interior problem of computed tomography using a priori knowledge. Inverse Problems. 2008;24 Article ID 065001, 27 pages. [PMC free article] [PubMed]
- Defrise M, Noo F, Clackdoyle R, Kudo H. Truncated Hilbert transform and image reconstruction from limited tomographic data. Inverse Problems. 2006;22:1037–53.
- Donoho DL. Compressed sensing. IEEE Transactions on Information Theory. 2006;52:1289–306.
- Kudo H, Courdurier M, Noo F, Defrise M. Tiny a priori knowledge solves the interior problem in computed tomography. Phys Med Biol. 2008;53:2207–31. [PubMed]
- Natterer F. The Mathematics of Computerized Tomography. Philadelphia: Society for Industrial and Applied Mathematics; 2001.
- Noo F, Clackdoyle R, Pack JD. A two-step Hilbert transform method for 2D image reconstruction. Physics in Medicine and Biology. 2004;49:3903–23. [PubMed]
- Parker DL. Optimal short scan convolution reconstruction for fan beam CT. Med Phys. 1982;9:254–7. [PubMed]
- Rudin LI, Osher S, Fatemi E. Nonlinear Total Variation Based Noise Removal Algorithms. Physica D. 1992;60:259–68.
- Wang G, Jiang M. Ordered-Subset Simultaneous Algebraic Reconstruction Techniques (OS-SART) Journal of X-ray Science and Technology. 2004;12:169–77.
- Wang G, Li Y, Jiang M. Uniqueness theorems in bioluminescence tomography. Medical Physics. 2004;31:2289–99. [PubMed]
- Wang G, Ye YB, Yu HY. 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.
- Wang G, Yu HY. Methods and Systems for Exact Local CT Based on Compressive Sampling. Patent disclosure submitted to Virginia Tech Intellectual Properties on Dec. 20 2008.
- Ye YB, Yu HY, Wang G. Exact Interior Reconstruction with Cone-beam CT. International Journal of Biomedical Imaging. 2007a;2007 Article ID: 10693, 5 pages. [PMC free article] [PubMed]
- Ye YB, Yu HY, Wei Y, Wang G. A general local reconstruction approach based on a truncated Hilbert transform. International Journal of Biomedical Imaging. 2007b;2007 Article ID: 63634, 8 pages. [PMC free article] [PubMed]
- Ye YB, Yu HY, Wang G. Exact interior reconstruction from truncated limited-angle projection data. International Journal of Biomedical Imaging. 2008;2008 Article ID: 427989, 6 pages. [PMC free article] [PubMed]
- Yu HY, Wang G. A general formula for fan-beam lambda tomography. International Journal of Biomedical Imaging. 2006;2006 Article ID: 10427, 9 pages.
- Yu HY, Ye YB, Wang G. Local Reconstruction Using the Truncated Hilbert Transform via Singular Value Decomposition. Journal of X-Ray Science and Technology. 2008;16:243–51. [PMC free article] [PubMed]

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