|Home | About | Journals | Submit | Contact Us | Français|
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() = 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:
where h(s) represents the ramp filtering kernel. Given the property that h(s) is analytic everywhere except at s = 0, f(s,) will be smooth and bounded for |s| < ρ0 because there is no singular operation in the integral. Actually, f(s,) is an analytic function of s for |s| < ρ0. In the backprojection step, the backprojection image g(ρ,θ) can be written as:
Eq. (2) shows that the reconstruction g(ρ,θ) only involves the data f(s,) with |s|≤ρ. That is, the reconstruction of g(ρ,θ) with ρ < ρ0 only involves the data f(s,) with |s|≤ ρ0. Hence, g(ρ,θ) with ρ < ρ0 is a smooth and bounded function.
For the reconstructed image g(ρ,θ) with |ρ|< ρ0, both and are smooth and bounded.
Perform a derivation operation on both side of Eq. (2), we have
Using the analytic property of f(s,) with the variable s for |s|< ρ0 mentioned in the proof of Lemma 2.1, we complete this proof.
Since Natterer has proved that g(ρ,θ) is a non-zero function, it is impossible that both and 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:
where μ(ρ,θ) represents a sparifying transform. For the commonly used gradient transform in the medical imaging field, we have
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:
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 λ = 0 for the given f(ρ,θ) and g(ρ,θ), the exactness of interior tomography in the CS framework should hold in this particular case.
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 for ρ < ρ0, the TV of f(ρ) inside the ROI can be computed as:
Clearly, tv reaches the minimal at λ = −1 instead of λ = 0 in this special case, which completes the proof.
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
where tm > 0 indicates a location of discontinuity, qm characterizes the corresponding jump defined as
and δ(•) represents Dirac delta function defined as:
For two bounded functions g1(ρ) and g2(ρ) with 0 < ρ < ρ0, there will be:
where qm and tm are the same as in Eq. (10).
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
Denoting the set of finite discontinuous points of f(ρ) at tm as ID, and the interval [0, ρ0] as I, we express the TV expansion Eq. (7) as
Eq. (20) can be simplified as
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 L1, L2,…Lm…,LM. While the length of Lm is denoted as tm > 0, the difference between the two neighboring sub-regions is denoted as qm, which is defined in the same fashion as in Eq. (11). An example is given in Figure 1, where the compact support Ω includes 6 sub-regions and the ROI Ω0 covers 8 arc-segments as boundaries. For a given point = (ρ,θ) on an arc-segment Lm, we set up a local coordinate system with the unit direction being the tangent direction of Lm on and unit the direction being the normal direction of Lm on (see Figure 2). By the assumption of piecewise constant, f() has the maximal directional derivation along the direction . What is more, the absolute value of the maximal directional derivative μ in a small neighborhood of can be written as
On the other hand, we have a global polar system with the unit vectors being and . If we denote the angle between and as γ, the gradient of f() along in a small neighborhood of can be written as qmδ(v)(cos γ+sin γ). When the image f() is overlapped with a smooth and bounded function g() as we mentioned above in Lemmas 2.1 and 2.2., the absolute value of maximal directional derivative along the direction for the composited function should be:
For two bounded functions g1() and g2(), there will be:
In the CS framework, an interior ROI of a general compactly supported function f() can be exactly determined by minimizing the total variation defined by Eq.(7) if f() can be decomposed into finitely many constant sub-regions.
For an ROI Ω0 of a general piecewise constant function f(), the TV formula Eq. (5) becomes
where the finite many undifferentiable points do not affect the integral result. After the introduction of g() as discussed in Lemmas 2.1 and 2.2, the TV expression Eq. (7) becomes
Using the same strategy for computing and Lemma 2.4, we immediately obtain that
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 fm,n = f(mΔ, nΔ) based on all the truncated local projections, where Δ represents the sampling interval, m and n are integers. In the second step, we minimize the 1 norm for a given sparsifying transform of the discrete image fm,n using the standard steepest descent method. These two steps were iteratively performed in an alternating manner. Specifically, the algorithm can be summarized in the following pseudo-code:
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 PTV for the minimization of TV, where α represents the maximal step for the steepest descent, αs the decreasing scale of α after each computation, and PTV the local loop time. Line 2 initializes the reconstructed image, the main loop count k and the number of subsets PART for OS-SART. Lines 5-13 perform the OS-SART reconstruction. Line 6 updates using the OS-SART technique, where only the projections in one subset are used. Lines 7-12 define the local loop to minimize the 1 norm. Line 8 computes the steepest decent direction dm,n, of fTV. For the discrete gradient transform, the magnitude of the gradient can be approximately expressed as:
Correspondingly, the total variation can be defined as . Then, we have the steepest descent direction defined by
For another sparsifying transform, we can deduce the corresponding dm,n, easily. Line 9 computes the normalized factor β which is the ratio between the maximum absolute values of and dm,n. Lines 10 and 11 update and α, respectively. Line 14 decides if the main loop should be stopped or not. The stopping criteria may include the iteration time, convergence speed, total variation, image quality, etc.
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 dm,n, using Eq. (35), we added a small constant to Eq.(34) when computing the gradient μm,n. That is,
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, (x0, y0) the center of the ellipse, ω denotes the rotation angle, f the relative attenuation coefficient. The units for a b, and (x0, y0) 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.
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:
where Ck 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
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:
where and 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
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:
In the CS framework, an interior ROI can be exactly determined by minimizing the 1 norm if we define a sparsifying transform as
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).