PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Phys Med Biol. Author manuscript; available in PMC 2010 April 22.
Published in final edited form as:
PMCID: PMC2858336
NIHMSID: NIHMS192195

Supplemental Analysis on Compressed Sensing Based Interior Tomography

Abstract

Recently, in the compressed sensing framework we proved that an interior ROI can be exactly reconstructed via the total variation minimization if the ROI is piecewise constant. In the proofs, we implicitly utilized the property that if an artifact image assumes a constant value within the ROI then this constant must be zero. Here we prove this property in the space of square integrable functions.

Keywords: Computed tomography (CT), interior tomography, compressed sensing, total variation minimization, artifact image

I. Background

It is well known that an interior region-of-interest (ROI) cannot be uniquely reconstructed from projection data only associated with lines through the ROI. Interestingly, our group and others independently established that this interior problem can be exactly and stably solved if a sub-region in the ROI is precisely known (Ye et al., 2007; Yu et al., 2008; Courdurier et al., 2008; Kudo et al., 2008). Recently, in the compressed sensing framework we proved that an interior ROI can be exactly reconstructed via the total variation minimization if the imaging object is piecewise constant (Yu and Wang, 2009; Han et al., 2009). In our proofs of Theorems 2.2 and 2.3, after we made the conclusion that the minimum was achieved at λ = 0 in (Yu and Wang, 2009) we implicitly utilized the property that if an artifact image assumes a constant value within the ROI then this constant must be zero. Here we provide a proof for this property so that our compressed sensing based interior tomography theory is expressed rigorously.

Let us consider a 2D smooth image f(ρ,θ), ρ[set membership][0,1], θ[set membership][0,2π) on a compact support within a unit disk Ω . The Radon transform of f(ρ,θ) can be written as R(s,[var phi]), s [set membership][–1,1], [var phi][set membership][0,π) . Suppose that we are only interested in reconstructing its interior part defined as ρ<a (0<a<1) from the corresponding local Radon transform data R(s,[var phi]),|s|<a , which is also referred to as purely local projections. Based on the classic analysis, given these local data there is in general no unique solution (Natterer, 2001). Thus, any reconstructed image from such a local dataset can be viewed as a superposition of the exact reconstruction from a complete dataset R(s,[var phi]),s[set membership][–1,1],[var phi][set membership][0,π) and an artificial reconstruction, which is an artifact image corresponding to another projection dataset R(s,[var phi]), a <|s|≤1,[var phi][set membership][0,π) . Although R(s,[var phi])=0 for |s|<a, it can still produce a non-zero 2D local image g(ρ,θ), ρ[set membership][0,a), θ[set membership][0,2π) , inside the ROI. That is why we have the non-uniqueness of this interior problem (Natterer, 2001). It is well known that g(ρ,θ) is smooth and bounded inside the ROI if R(s,[var phi]) is continuous and bounded. Previously, we proved that both g(ρ,θ)ρ and g(ρ,θ)ρθ are smooth and bounded inside the ROI (Yu and Wang, 2009). In the this note we will prove that g(ρ,θ)[equivalent]C for ρ[set membership][0,a) if and only if C=0.

The organization of this report is as follows. In the next section, we first prove that the constant C must be zero in the symmetric case using Taylor's expansion technique and Weierstrass's approximation theorem, and then extend it to the general non-symmetric case. In the third section, we numerically verify that a non-zero constant cannot satisfy the so-called Picard condition in the framework of singular value decomposition. In the last section, we discuss the relationship between this report and our previous work as well as future directions

II. Theoretical Analysis

Assuming that a circularly symmetric artifact image g(ρ,θ)=g(ρ) is reconstructed from a projection dataset R(s,[var phi]) = R\(s) (a <|s|≤1,[var phi][set membership][0,π)) and R(s)[equivalent]0 for s [set membership][–a,a] (Figure 1) and g(ρ)[equivalent]C (ρ[set membership][0,a) ), where C is a non-zero constant. Then, by the definition of the line integral we must have the following relationship for g(ρ) (ρ[set membership][a,1]):

R~(s)=a1g(ρ)dρ2s2+Ca2s2=a1ρg(ρ)dρρ2s2+Ca2s2=0,s[0,a].
(1)
Figure 1
Circularly symmetric artifact image g.

Our main theoretical results are as follows.

Theorem 2.1. Assuming that a circularly symmetric artifact image g(ρ,θ) = g(ρ) is reconstructed from a projection dataset R(s,[var phi])=R(s) (|s|≤1,[var phi][set membership][0,π)) and R(s) [equivalent] 0 for s [set membership][–a,a]. If g(ρ) is a square integrable function on [0,1] with g(ρ)[equivalent]C (ρ[set membership][0,a) ), then C=0.

Proof: Eq. (1) can be changed to

Ca1(sa)2+a1g(ρ)dρ1(sρ)2=0,0sa.
(2)

Denote (2k–1)!!=1×3×5---× (2k–1) and k!=1×2×3---×k, we have the following Taylor's expansions

11(sρ)2=1+k=1(2k1)!!2kk!(sρ)2k,
(3)
1(sa)2=112k1k=1(2k1)!!2kk!(sa)2k.
(4)

Inserting Eqs. (3) and (4) into Eq. (2), we have

a1g(ρ)(aρ)2kdρ=Ca2k1,k0.
(5)

Making the variable transformation t=a/ρ, Eq. (5) becomes

a1g(at)t2k2dt=C2k1,k0.
(6)

Let

g~(t)={g(at)Ct[a,1]Ct[0,a)},
(7)

Eq. (6) implies

01g~(t)t2kdt=0,k0.
(8)

Considering the even function ĝ(t)=(g(|t|) on the interval [–1,1], we have

11g^(t)tkdt=0,k0.
(9)

By Weierstrass's approximation theorem and the fact that the continuous function space is dense in the square integrable function space on [–1,1], there exist a polynomial series {pn(t)}n=1 on interval [–1,1] such that

11g^(t)Pn(t)2dt1n.
(10)

This leads to

11g^2(t)dt=limn11g^(t)pn(t)dt=0.
(11)

Hence, ĝ(t)=0 almost everywhere on [–1,1], which implies C= 0.

The above theorem can be extended to the general non-symmetric case.

Theorem 2.2. Assuming that a circularly non-symmetric artifact image g(ρ,θ) is reconstructed from a projection dataset R(s,[var phi]) (|s|≤1,[var phi][set membership][0,π) ) and R(s,[var phi])[equivalent]0 for s[set membership][–a,a]. If g(ρ,θ) is a square integrable function on ρ[set membership][0,1] and θ[set membership][0,2π] with g(ρ,θ)[equivalent]C (ρ[set membership][0,a) ), then C=0.

Proof: Let us denote g(ρ,θ) in the Cartesian coordinates as gd(x, y) with x = ρcosθ and y = ρsinθ. Define g(ρ,θ)=12π02πg(ρ,θ+ϕ)dϕ, correspondingly we have

gd(x,y)=12π02πgd(ρcos(θ+ϕ),ρsin(θ+ϕ))dϕ=12π02πgd(xcosϕysinϕ,xsinϕ+ycosϕ)dϕ.
(12)

Clearly, g(ρ,θ) is a circularly symmetric square integrable function and satisfies g(ρ,θ)[equivalent]C (ρ[set membership][0,a) ). For s [set membership][–a,a] the Radon transform R(s,[var phi]) of g(ρ,θ) can be expressed as

R(s,φ)=1s21s2gd(s(cosφ,sinφ)+t(sinφ,cosφ))dt=1s21s2gd(scosφtsinφ,ssinφ+tcosφ)dt=1s21s2(12π02πgd(scos(φ+ϕ)tsin(φ+ϕ),ssin(φ+ϕ)+tcos(φ+ϕ))dϕ)dt=12π02π(1s21s2gd(scos(φ+ϕ)tsin(φ+ϕ),ssin(φ+ϕ)+tcos(φ+ϕ))dt)dϕ=12π02π(1s21s2gd(s(cos(φ+ϕ),sin(φ+ϕ))+t(sin(φ+ϕ),cos(φ+ϕ)))dt)dϕ=12π02πR~(s,φ+ϕ)dϕ=0,
(13)

where Eq.(12) has been used. By Theorem 2.1., we have C=0.

Remark: For real-world applications, all the reconstructed CT artifact images are bounded function and hence square integrable. Therefore, our assumption of square integrable functions can cover all the practical cases.

III. Numerical Analysis

Assuming that an integral kernel K(s,ρ) defined on Is×Iρ is continuous, square integrable and non-degenerate, by means of the singular value expansion (see Chapter 1.2.1 in (Hansen, 1998) ), we have

K(s,ρ)=i=1σiui(s)vi(ρ),
(14)

where σ1 ≥σ2 ≥---> 0 , and ui(s) and vi(ρ) are orthonormal basis functions with respect to the conventional inner product

(ui,uj)=Isui(t)uj(t)d={1ifi=j0ifij},
(15)

and

(vi,vj)=Iρvi(t)vj(t)dt={1ifi=j0ifij}.
(16)

For any continuous and square integral function q(s), s[set membership]Is, the Fredholm integral equation of the first kind

IρK(s,ρ)g(ρ)dρ=q(s)
(17)

can be expressed as

i=1σi(vi,g)ui(s)=i=1(ui,q)ui(s),sIs.
(18)

If there exists solutions for Eq. (17), one of them would be in the following form:

g(ρ)=i=1(ui,q)σivi(ρ),ρIρ.
(19)

Eqs. (18) and (19) hold for uniform convergence (see Chapter 2.4 and Theorem 8.3.2 in (Smithies, 1958)). Although K(s,ρ) and q(s) are continuous and square integrable, g(ρ) expressed by Eq. (19) is not necessarily continuous or square integrable. In order that there exists a square integrable solution g(ρ) to Eq. (17), the right side of Eq. (19) must satisfy the Picard condition (see Chapter 1.2.3 in (Hansen, 1998)):

i=1((ui,q)σi)2<.
(20)

This means that the (ui, q) must decay faster than the singular value σi. This requirement is identical to that q(s) must belong to the range of the kernel K(s,ρ).

Let K(s,ρ)=ρρ2s2, Iρ = [a, 1], Is = [0, a] and q(s)=Ca2s2, Eq. (1) can be considered as a a Fredholm integral equation of the first kind as defined in Eq. (17). If there exists a square integrable solution g(ρ) to Eq. (1), it should satisfy the Picard condition. However, it is not easy to perform a thorough theoretical analysis to verify the Picard condition for our problem. Nevertheless, the discrete Picard condition is often used in practical applications to study the existence of a solution to Eq. (17) using the singular value decomposition (SVD) approach. As pointed out by Hansen (Hansen, 1988), although a discrete Picard condition does not exist in a strict mathematical sense, it is certainly informative when numerically solving a Fredholm integral equation of the first kind Eq.(17).

In our case, let us discrete s as sn=a(n0.5)N, n = 1,2,---,N and ρ as ρm=a+(m0.5)(1a)M, m = 1,2,---,M. Accordingly, K(s,ρ), g(ρ) and q(s) are discretized as K(snm), g(ρm), q(sn). Hence, we have

KG=Q,
(21)

where

K=1aM(K(s1,ρ1)K(s1,ρm)K(s1,ρM)K(sn,ρ1)K(sn,ρm)K(sn,ρM)K(sN,ρ1)K(sN,ρm)K(sN,ρM)),
(22)
G=(g(ρ1)g(ρm)g(ρM))T,
(23)

and

Q=(q(s1)q(sn)q(sN))T.
(24)

Let M=N, the matrix K has the following SVD form (Hansen, 1988):

K=UΛVT=m=1MσmumvmT,
(25)

where U=(u1 --- um --- uM) and V = (v1 --- vm --- vM) are orthogonal matrices, and Λ is a diagonal matrix whose diagonal elements satisfying σ1≥σ2≥---≥ 0. The solution to Eq. (20) can be written as

G=i=mMumTQσmvm.
(26)

According to the Picard condition, we can compute umTQσm to test the existence of a solution to Eq. (1).

Without loss of generality, we set a=0.5 and q(s)=a2s2. For different combinations of sampling point numbers M and N, we computed and plotted lgumTQσm in Figure 2, which clearly show that the singular value σm decays faster than umTQ, suggesting that there is no square integrable solution for C ≠ 0 because any non-zero constant does not lead to a sufficiently fast decay of umTQ. These numerical results are consistent to our theoretical analysis in the preceding section that if g(ρ) is square integrable constrained by an interior scan dataset then we must have C = 0 in Eq. (1).

Figure 2
Curves of lgumTQσm for different sampling point numbers (a) M =N=50, (b) 100, (c) 200 and (d) 400, respectively.

IV. Discussions and Conclusion

Although we assumed that an object to be reconstructed was piecewise constant both inside and outside an ROI, the piecewise constancy is only needed inside the ROI in the proofs in (Yu and Wang, 2009; Han et al., 2009). In other words, the piecewise constancy outside the ROI is not necessary. Note that a generalized total variation minimization theorem was proved for compressed sensing based interior tomography in (Han et al., 2009), without involving the Dirac delta function. Because most objects in CT applications can be approximately modeled as piecewise constant, the piecewise constancy assumption is quite reasonable (Wang et al., 2004), and our compressed sensing based interior tomography approach is practically useful. Most importantly, our finding suggests a new direction of interior tomography, and may lead to powerful methods for better interior reconstruction. Currently, we are working to relax the piecewise constancy assumption so that real-world images can be modeled with higher accuracy for new interior reconstruction results.

As discussed in (Yu and Wang, 2009), the key idea of compressed sensing based interior tomography is to define an appropriate sparsifying transform and an associated objective function. Then, the minimization of the objective function will lead to the true ROI image without any ambiguity. For that purpose, the most important step is to prove that the artifact image inside the ROI must be zero when its total variation is zero. Actually, it might appear in (Yu and Wang, 2009) that the artifact image inside the ROI can be any constant when the total variation is zero. Hence, it is the major contribution of this report that has proved the constant being necessarily zero, providing the uniqueness of our interior reconstruction. Using this methodology, we are working to extend compressed sensing based interior reconstruction under different conditions and for other modalities.

In conclusion, we have studied the characteristic of the artifact image inside an interior ROI subject to the local projection data constraint. Our theoretical and numerical analysis has demonstrated that the artifact image inside the ROI can be expressed as a constant if and only if the constant is zero. This analysis has consolidated our proofs of Theorems 2.2 and 2.3 in (Yu and Wang, 2009).

Acknowledgment

This work is partially supported by NIH/NIBIB grants (EB002667, EB004287, EB007288), a grant from Toshiba Medical Research Institute USA, Inc., NSFC (60325101, 60532080, 60628102), Key Laboratory of Machine Perception (Ministry of Education) of Peking University, and Microsoft Research of Asia. The first author thanks Dr. Per Christian Hansen for offering a hardcopy of his paper (Hansen, 1988).

References

  • 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]
  • Han W, Yu H, Wang G. A total variation minimization theorem for compressed sensing based tomography, under review. 2009. [PMC free article] [PubMed]
  • Hansen PC. Computation of the singular value expansion. Computing. 1988;40:185–99.
  • Hansen PC. Rank-Deficient and Discrete Ill-Posed Problems: Numerical Aspects of Linear Inversion. Society for Industrial and Applied Mathematics; Philadelphia: 1998.
  • 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. Society for Industrial and Applied Mathematics; Philadelphia: 2001.
  • Smithies F. Integral Equations. Cambridge University Press; 1958.
  • Wang G, Li Y, Jiang M. Uniqueness theorems in bioluminescence tomography. Medical Physics. 2004;31:2289–99. [PubMed]
  • Ye Y, Yu H, Wei Y, Wang G. A general local reconstruction approach based on a truncated Hilbert transform. International Journal of Biomedical Imaging. 2007;2007 Article ID: 63634, 8 pages. [PMC free article] [PubMed]
  • Yu H, Wang G. Compressed sensing based Interior tomography. Phys Med Biol. 2009;54(9):2791–2805. [PMC free article] [PubMed]
  • Yu H, Ye Y, 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]