PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Inverse Probl. Author manuscript; available in PMC 2010 April 20.
Published in final edited form as:
Inverse Probl. 2010 January 1; 26(3): 350131–3501329.
doi:  10.1088/0266-5611/26/3/035013
PMCID: PMC2857378
NIHMSID: NIHMS192276

High Order Total Variation Minimization for Interior Tomography

Abstract

Recently, an accurate solution to the interior problem was proposed based on the total variation (TV) minimization, assuming that a region of interest (ROI) is piecewise constant. In this paper, we generalize that assumption to allow a piecewise polynomial ROI, introduce the high order TV (HOT), and prove that an ROI can be accurately reconstructed from projection data associated with x-rays through the ROI through the HOT minimization if the ROI is piecewise polynomial. Then, we verify our theoretical results in numerical simulation.

Keywords: Interior tomography, compressive sensing, total variation (TV), high order TV

I. Introduction

While classic CT theory targets exact reconstruction of a whole cross-section or an entire object from un-truncated projections, practical applications often focus on much smaller regions of interest (ROIs). Classic CT theory cannot exactly reconstruct an internal ROI only from truncated projections associated with x-rays through the ROI because this interior problem does not have a unique solution [1-2]. When applying traditional CT algorithms for interior reconstruction from truncated projection data, features outside the ROI may seriously disturb inside features and often hide diagnostic information.

Over past decades, lambda tomography has been extensively studied [3-5][6-10] but it only targets gradient-like features and has not been well received in the biomedical community. In recent work, Ye et al. found that the interior problem can be exactly and stably solved if a sub-region in the ROI is known [11]. Similar results were also independently reported by other researchers [12-13]. However, it can be difficult to obtain precise prior knowledge of a sub-region in many cases such as cardiac CT perfusion studies.

The classic Nyquist sampling theorem states that one must sample a signal at least twice faster than its bandwidth to keep all the information. Surprisingly, an emerging theory – compressive sensing (CS) – promises to capture compressible signals at a much lower rate than the Nyquist rate and yet allows accurate reconstruction of such a signal from a limited number of samples [14-15]. The main idea of CS is that most signals are sparse in appropriate orthonormal systems, that is, a majority of their coefficients are close or equal to zero. Typically, CS starts with taking a limited amount of representative measures, and then recovers a signal with an overwhelming probability via the 1 norm minimization. Since the number of samples is limited, recovery of the signal or image would involve solving an underdetermined matrix equation. That is, there are many candidate solutions that can all well fit the incomplete measurement. Thus, some additional constraint is needed to select the “best” candidate. While the classic solution is to minimize the 2 norm, Donoho, Candes, Romberg and Tao showed that finding the candidate with the minimum 1 norm, which leads to the TV minimization in some cases, is most likely the correct choice [14-15].

Under the guidance of CS theory, we recently demonstrated that exact interior reconstruction could be achieved if an ROI is piecewise constant [16-19]. Inspired by these results, here we extend their work to utilize more general prior knowledge for interior tomography. Specifically, we will prove that if an ROI is piecewise polynomial, then the interior problem can be accurately solved through the high order TV (HOT) minimization.

The organization of this paper is as follows. In the next section, the interior problem will be defined. In the third section, our theoretical findings will be presented. In the fourth section, numerical results will be described to support our methodology. In the last section, related issues will be discussed.

II. Interior Problem

Without loss of generality, we assume the following conditions throughout this paper:

  • Condition (1) An object image f0(x) is compactly supported on a disc ΩA={x=(x1,x2)R2:x<A}, where A is a positive constant. Furthermore, f0(x) is a piecewise smooth function; that is, ΩA can be partitioned into finitely many sub-domains {Dj}j=1N0, such that f0(x) is smooth with bounded derivatives in each Dj;
  • Condition (2) An internal ROI is a smaller disc Ωa={x=(x1,x2)R2:x<a}, as shown in Fig. 1, where a is a positive constant and a < A;
    Figure 1
    Region of interest (ROI) inside a compact support.
  • Condition (3) Projections through the ROI
    Rf0(s,θ)=f0(sθ+tθ)dt,a<s<a,θS1,
    (1)
    are available, where θ = (cos [var phi], sin [var phi]), θ[perpendicular] = (−sin [var phi], cos [var phi]), 0 ≤ [var phi] < 2π.

The interior problem is to find an image f(x) such that

  • Condition (4) f(x) is a piecewise smooth function and compactly supported on the disc ΩA;
  • Condition (5) Rf(s,θ) = Rf0(s,θ), −a < s < a, θ [set membership] S1.

It is well known that under the conditions (4) and (5) the interior problem does not have a unique solution [1-2]. The following theorem characterizes the structure of solutions to the interior problem.

Theorem 1. Any image f(x) satisfying Conditions (4) and (5) can be written as f(x) = f0(x) + u(x) for xR2, where u(x) is an analytic function in the disc Ωa, and

Ru(s,θ)=0,a<s<a,θS1.

We call such an image f(x) a candidate image, and correspondingly u(x) an ambiguity image.

Proof: Let u(x) = f(x)–f0(x). Clearly, u(x) is a piecewise smooth function and compactly supported on the disc ΩA, and

Ru(s,θ)=0,a<s<a,θS1.
(2)

By the well-known Radon inversion formula, we have

u(x)=14π202πsasRu(θ,s)xθsdsdφ,xΩa.
(3)

By the Taylor expansion

1xθs=1s(1xθs)=k=0(xθ)ksk+1,
(4)

we obtain

u(x)=k=0k1+k2=kck1,k2x1k1x2k2.
(5)

Hence, the function u(x) is an analytic function in Ωa. Examples of non-zero ambiguity images can be found in [1].

From now on, let u(x) always represent an ambiguity image unless otherwise stated. Theorem 1 suggests that the interior problem can be solved by utilizing appropriate a priori knowledge. In the following, we will introduce the concept of the high order TV (HOT) and minimize it to solve the interior problem exactly under the assumption that an image f0(x) is a piecewise polynomial function in an ROI Ωa.

III. Theoretical Analysis

If an object image f0(x) is a piecewise polynomial function in an ROI Ωa, we can prove that f0(x) is the only candidate image that minimizes the HOT to be defined by Eqs. (44-45) and (81-82). First, let us prove that if an ambiguity image is a polynomial function in Ωa, then it must be zero. This result will be formally stated as Theorem 2. In order to prove Theorem 2, we will need Lemmas 1 and 2.

Lemma 1. Suppose that a is a positive constant. If

  • (a) g(z) is an analytic function in C\(,a][a,+);
  • (b) p(x) is a polynomial function;
  • (c) g(x)=1πp.vt<ap(t)xtdt, for x [set membership] (−a,a),

then limy0+Im(g(x+iy))=p(x), for x(,a)(a,+).

Proof: Since g(z) is uniquely determined by its value on interval (−a,a) and

g(x)=1πp.vt<ap(t)xtdt,forx(a,a),
(6)

g(z) is uniquely determined by p(x) for x [set membership] (−a,a). Let g(z) be denoted by gn(z) with respect to p(x) = xn. Because of the linear relationship between g(z) and p(x), it suffices to prove the lemma in the case p(x) = xn (n = 0, 1, 2,···). That is, we need to prove

limy0+Im(gn(x+iy))=xn,forx(,a)(a,+),n=0,1,2,.
(7)

Let us proceed by induction with respect to n . When n = 0 , that is, p(x) = 1, we have

g0(x)=1πp.vt<a1xtdt=1π(ln(a+x)ln(ax)),forx(a,a).
(8)

Hence,

g0(z)=1π(ln(a+z)ln(az)),forzC\(,a][a,+),
(9)

where the principal value interval for the logarithm function is (−π,π). That is,

ln(a+z)=lna+z+iarg(a+z),π<arg(a+z)<π,forzC\(,a],
(10)

and

ln(az)=lnaz+iarg(az),π<arg(az)<π,forzC\[a,+).
(11)

As shown in Fig. 2, arg(az1) approaches −π when arg(a + z1) is close to −π, and arg(az2) approaches π when arg(a + z2) is close to π. Clearly,

Im(g0(z))=1π(arg(a+z)arg(az)),forzC\(,a][a,+).
(12)

Since

limy0+arg(a+x+iy)=π,limy0+arg(axiy)=0,forx(,a)
(13)

and

limy0+arg(a+x+iy)=0,limy0+arg(axiy)=π,forx(a,),
(14)

we have

limy0+Im(g0(x+iy))=1,forx(,a)(a,+).
(15)

For n > 0 and x [set membership](−a,a), we have

gn(x)=1πp.vt<atnxtdt=1πp.vt<a(tnxtn1)+xtn1xtdt=1πt<atn1dt+xπp.vt<atn1xtdt=1nπ(an(a)n)+xgn1(x).
(16)

Therefore,

gn(z)=1nπ(an(a)n)+zgn1(z),forzC\(,a][a,+),
(17)

and

Im(gn(z))=Im(zgn1(z)),forzC\(,a][a,+).
(18)

For x(,a)(a,+) and y > 0 , we have

Im(gn(x+iy))=xIm(gn1(x+iy))+yRe(gn1(x+iy)),
(19)

and

limy0+Im(gn(x+iy))=xlimy0+Im(gn1(x+iy))=x2limy0+Im(gn2(x+iy))==xnlimy0+Im(g0(x+iy))=xn,
(20)

which completes the proof of Lemma 1.

Figure 2
Illustration of (a) arg(a+z) and (b) arg(a-z).

Lemma 2. Suppose that a and A are positive constants with a < A. If a function v(x) satisfies

  • (d) v(x) is bounded with supp v [subset, dbl equals] [−A,A];
  • (e) v(x) = p(x) for x [set membership] (−a,a), where p(x) is a polynomial function;
  • (f) Hv(x) = 0, for x [set membership] (−a,a), where Hv(x) is the Hilbert transform of v(x), that is,
    Hv(x)=1πp.v.Rv(t)xtdt,
    then v(x) = 0, for x [set membership] (−∞,∞).

Proof: The function (see Corollary 4.1.2 on page 40 in [20])

g(z)=1πtav(t)tzdt,
(21)

is analytic on C\(,a][a,+). Since

1t(x+iy)=tx+iy(tx)2+y2,
(22)

we have

Im(g(x+iy))=1πtay(tx)2+y2v(t)dt=1πRy(tx)2+y2v~(t)dt=1πR1(txy)2+1v~(t)d(txy)=1πR1τ2+1v~(x+τy)dτ,fory>0,
(23)

where

v~(x)={v(x),x(,a)(a,)0,x[a,a]}.
(24)

Because v(x) can be written as

v~(x)=1πR1τ2+1v~(x)dτ,
(25)

we have

limy0+RIm(g(x+iy))v~(x)dxlimy0+{1πRv~(+τy)v~()1τ2+1dτ}=0,
(26)

where

v~(+τy)v~()1=Rv~(x+τy)v~(x)dx.
(27)

Therefore,

limy0+Im(g(x+iy))=v~(x)=v(x),
(28)

for x(,a)(a,) almost everywhere (a.e.).

By Condition (f), we have

g(x)=1πp.vRv(t)txdt1πp.vt<av(t)txdt=Hv(x)+1πp.vt<ap(t)xtdt=1πp.vt<ap(t)xtdt,forx(a,a).
(29)

By Lemma 1, we obtain

limy0+Im(g(x+iy))=p(x),forx(,a)(a,+).
(30)

Combining Eq. (28), (30) and Condition (e), we have

v(x)=p(x),fora.e.x(,).
(31)

Condition (d) and Eq. (31) imply p(x) [equivalent] 0. Hence, v(x) = 0, for a.e. x [set membership] (−∞,∞), which completes the proof.

Theorem 2. If an ambiguity image u(x) satisfies

  • (g) u(x) = p(x) for x [set membership]Ωa, where p(x) is a 2-D polynomial function;
  • (h) Ru(s,θ) = 0, for s [set membership] (−a,a), θ [set membership] S1,

then u(x) = 0.

Proof: As illustrated in Fig. 3, for an arbitrary [var phi]0 [set membership][0,π), let Lθ0 denote the line through the origin and tilted by θ0 = (cos [var phi]0, sin [var phi]0). When u(x) is restricted to the line Lθ0, it can be expressed as

vθ0(t)=u(t(cosφ0,sinφ0)),t(,).
(32)

By the relationship between the differentiated backprojection of projection data and the Hilbert transform of an image [21-22], we have

Hvθ0(t)=12πφ0π2φ0+π2Ru(s,θ)ss=xθ0tθdφ,
(33)

where xθ0t=t(cosφ0,sinφ0), θ = (cos [var phi], sin [var phi]).

Figure 3
Radial line through the origin.

By Condition (g) and Eq. (32), we have

vθ0(t)=p(t(cosφ0,sinφ0)),fort(a,a),
(34)

where p(t(cos [var phi]0, sin [var phi]0)) is a polynomial function with respect to t. By Condition (h) and Eq. (33), we have

Hvθ0(t)=0,fort(a,a).
(35)

By Lemma 2, we obtain that vθ0t = 0, which implies that u(x) = 0.

Now, let us analyze a candidate image under the assumption that f0(x) is a piecewise polynomial function in Ωa.

Theorem 3. Suppose that an object image f0(x) is a piecewise polynomial function in Ωa. If a candidate image f(x) is also a piecewise polynomial function in Ωa, then f(x) = f0(x).

Proof: By Theorem 1, we have that f(x) = f0(x) + u(x), where u(x) is an analytic function in Ωa. On the other hand, by the assumption of Theorem 3, u(x) = f(x) − f0(x) is also a piecewise polynomial function in Ωa. Combining these two facts, u(x) must be a polynomial function in Ωa. Then, by Theorem 2 we have u(x) = 0.

Theorem 3 suggests that f0(x) can be singled out among all candidate images by minimizing some kind of high order TV. For simplicity, let us first assume that an object image f0(x) is a piecewise linear function in Ωa; that is, Ωa can be decomposed into finitely many sub-domains: {Ωi}i=1m (Fig. 4) such that

f0(x)=fi(x)=aix1+bix2+ci,forxΩi,1im,
(36)

and each sub-domain Ωi is adjacent to its neighboring sub-domains Ωj with piecewise smooth boundaries Γi,j,, j [set membership] Ni.

Figure 4
ROI consisting of 7 sub-domains.

If m = 1, that is, f0(x) is a linear function in Ωa, then any candidate image f(x) = f0(x) + u(x) is a smooth function in Ωa. We define the second order TV for any candidate image f(x) of this type as

HOT2(f)=Ωamin(maxθS1dfdθ(x),maxθS1d2fdθ2(x))dx1dx2.
(37)

where θ = (cos [var phi], sin [var phi]), [var phi] [set membership] [0, 2 π), dfdθ(x)=dfx,θ(t)dtt=0, d2fdθ2(x)=d2fx,θ(t)dt2t=0 and fx,θ(t) = f(x + ), t [set membership] (−(ε,ε).

Note that if fx1, fx2 are continuous at x, then

dfdθ(x)=fx1(x)cosφ+fx2(x)sinφ.forθ=(cosφ,sinφ)S1.
(38)

Hence,

maxθS1dfdθ(x)=(fx1(x))2+(fx2(x))2.
(39)

There also exists a similar formula for maxθS1d2fdθ2(x), as represented in the following Lemma.

Lemma 3. Suppose that f(x) is a candidate image. If 2fx12, 2fx1x2 and 2fx22 are continuous at x, then

maxθS1d2fdθ2(x)=14(2fx12(x)2fx22(x))2+(2fx1x2(x))2+122fx12(x)+2fx22(x).
(40)

Proof: If 2fx12, 2fx1x2 and 2fx22 are continuous at x, then

d2fdθ2(x)=(x1cosφ+x2sinφ)2f(x)=2fx12(x)cos2φ+22fx1x2(x)cosφsinφ+2fx22(x)sin2φ=2fx12(x)(1+cos2φ2)+2fx1x2(x)sin2φ+2fx22(x)(1cos2φ2)=12(2fx12(x)+2fx22(x))+12(2fx12(x)2fx22(x))cos2φ+2fx1x2(x)sin2φ.
(41)

Thus, d2fdθ2(x) reaches its maximum when

(cos2φ,sin2φ)=sgn(12(2fx12(x)+2fx22(x)))(12(2fx12(x)2fx22(x)),2fx1x2(x))14(2fx12(x)2fx22(x))2+(2fx1x2(x))2.
(42)

This leads to Eq. (40).

Combining Eqs. (37), (39) and (40), we have

HOT2(f)=Ωamin{14(2fx122fx22)2+(2fx1x2)2+122fx12+2fx22,(fx1)2+(fx2)2}dx1dx2.
(43)

where f(x) = f0(x) + u(x), and f0(x) is a linear function in Ωa.

If m > 1, dfdθ(x) or d2fdθ2(x) may not exist on the boundaries i=1mj>i,jNiΓi,j Hence, we need modify the definition Eq. (37) of the second order TV to incorporate jumps of the image f(x) at Γi,j into the second order TV. To address the discontinuity of the first and second order derivatives across internal boundaries, we define the second order TV as the limit of the following sum:

HOT2(f)=lim sup(max1kMdiam{Qk})0k=1MIk(f),
(44)

where {Qk}k=1M is an arbitrary partition of Ωa, diam(Qk) the diameter of Qk, and

Ik(f)=min(supθS1supgC0(Qk),g(x)1Qkf(x)dgdθ(x)dx,supθS1supgC0(Qk),g(x)1Qkf(x)d2gdθ2(x)dx).
(45)

In the following Lemma, we will derive an explicit formula for HOT2(f ).

Lemma 4. Suppose that an object image f0(x) is a piecewise linear function in Ωa as defined by Eq. (36). For any candidate image f(x) = f0(x) + u(x), we have

HOT2(f)=i=1mj>i,jNiΓi,j(aiaj)x1+(bibj)x2+(cicj)ds+Ωamin{14(2fx122fx22)2+(2fx1x2)2+122fx12+2fx22,(fx1)2+(fx2)2}dx1dx2,
(46)

where the second term is a Lebesgue integral.

Proof: Note that f(x) = f0(x) + u(x), where u(x) is an analytic function. Let {Qk}k=1M be an arbitrary partition of Ωa. First, let us consider Qk which covers a common boundary Γi,j of a pair of neighboring sub-domains Ωi and Ωj, and is contained in ΩiΩj (Fig. 5 (a)). Let θN be the normal vector of the curve Γi,j pointing from Ωj towards Ωi, and <θ,θN > the angle between the vectors θ and θN. For gC0(Qk), we have

Qkf(x)dgdθ(x)dx=QkΩif(x)dgdθ(x)dx+QkΩjf(x)dgdθ(x)dx.
(47)

Performing 2-D integration by parts on the two terms of the right-hand side of Eq. (47) respectively, and utilizing the fact that g(x) = 0 near the boundary of Qk, we have

QkΩif(x)dgdθ(x)dx>Γi,jQk(fi(x)+u(x))g(x)cos<θ,θN>dsQkΩidfdθ(x)g(x)dx,
(48)
QkΩjf(x)dgdθ(x)dx=Γi,jQk(fj(x)+u(x))g(x)cos<θ,θN>dsQkΩjdfdθ(x)g(x)dx.
(49)

Inserting Eqs. (48) and (49) into Eq. (47), we obtain

Qkf(x)dgdθ(x)dx=Γi,jQk(fj(x)fi(x))g(x)cos<θ,θN>dsQkΩidfdθ(x)g(x)dxQkΩjdfdθ(x)g(x)dx.
(50)

We have

supgC0(Qk),g(x)1Γi,jQk(fj(x)fi(x))g(x)cos<θ,θN>ds=Γi,jQk(fj(x)fi(x))cos<θ,θN>ds,
(51)
supθS1supgC0(Qk),g(x)1Γi,jQk(fj(x)fi(x))g(x)cos<θ,θN>ds=Γi,jQkfj(x)fi(x)ds+o(1)Γi,jQk.
(52)

where o(1) represents an infinitesimal, satisfying lim(max1kMdiam{Qk})0o(1)=0. Besides, we have

QkΩidfdθ(x)g(x)dx=O(1)QkΩi,
(53)
QkΩidfdθ(x)g(x)dx=O(1)QkΩi,
(54)

where O(1) represents a quantity bounded by a constant which does not depend on Qk. Summing Eqs. (53-54) up, we obtain

QkΩidfdθ(x)g(x)dx+QkΩidfdθ(x)g(x)dx=O(1)Qk.
(55)

Combining Eqs. (50), (52) and (55), we obtain

supθS1supgC0(Qk),g(x)1Qkf(x)dgdθ(x)dx=Γi,jQkfj(x)fi(x)ds+O(1)Qk+o(1)Γi,jQk.
(56)

Repeatedly performing 2-D integration by parts, for gC0(Qk), we have

Qkf(x)d2gdθ2(x)dx=QkΩif(x)d2gdθ2(x)dx+QkΩjf(x)d2gdθ2(x)dx=Γi,jQk(fj(x)fi(x))dgdθ(x)cos<θ,θN>dsQkΩidfdθ(x)dgdθ(x)dxQkΩjdfdθ(x)dgdθ(x)dx=Γi,jQk(fj(x)fi(x))dgdθ(x)cos<θ,θN>dsΓi,jQk(dfjdθ(x)dfidθ(x))g(x)cos<θ,θN>ds+QkΩid2fd2θ(x)g(x)dx+QkΩjd2fdθ2(x)g(x)dx=Γi,jQk(fj(x)fi(x))dgdθ(x)cos<θ,θN>dsΓi,jQk(dfjdθ(x)dfidθ(x))g(x)cos<θ,θN>ds+Qkd2fdθ2(x)g(x)dx.
(57)

For the second term on the right-hand side of Eq. (57), we have

supgC0(Qk),g(x)1Γi,jQk(dfjdθ(x)dfidθ(x))g(x)cos<θ,θN>ds=Γi,jQkdfjdθ(x)dfidθ(x)cos<θ,θN>ds,
(58)
supθS1supgC0(Qk),g(x)1Γi,jQk(dfjdθ(x)dfidθ(x))g(x)cos<θ,θN>ds=Γi,jQkmaxθS1(dfjdθ(x)dfidθ(x)cos<θ,θN>)ds+o(1)Γi,jQk.
(59)

For the third term on the right-hand side of Eq. (57), we have

Qkd2fdθ2(x)g(x)dx=O(1)Qk.
(60)

Next, let us evaluate supθS1supgC0(Qk),g(x)1Qkf(x)d2gdθ2(x)dx in the following two cases.

Figure 5
Sub-domain Qk of (a) the first and (b) second types.

Case 1: If there exists some xΓi,jQk such that fj(x) − fi(x) ≠ 0, then for θ [set membership] S1, an arbitrary real constant C and an arbitrary compact set KQk, there exists a gC0(Qk) satisfying |g(x|≤1 and

dgdθ(x)=C(fj(x)fi(x)),forxΓi,jK,and
(61)

we have

supθS1supgC0(Qk),g(x)1Γi,jQk(fj(x)fi(x))dgdθ(x)cos<θ,θN>ds=+.
(62)

Combining Eqs. (59), (60) and (62), we obtain

supθS1supgC0(Qk),g(x)1Qkf(x)d2gdθ2(x)dxsupθS1supgC0(Qk),g(x)1Γi,jQk(fj(x)fi(x))dgdθ(x)cos<θ,θN>dssupθS1supgC0(Qk),g(x)1Γi,jQk(dfjdθ(x)dfidθ(x))g(x)cos<θ,θN>dssupθS1supgC0(Qk),g(x)1Qkd2fdθ2(x)g(x)dx=+.
(63)

Case 2: Otherwise, fj(x) − fi(x) = 0 for any xΓi,jQk, and the first term on the right-hand side of Eq. (57) equals 0. Then, we have

supθS1supgC0(Qk),g(x)1Qkf(x)d2gdθ2(x)dx=Γi,jQkmaxθS1(dfjdθ(x)dfidθ(x)cos<θ,θN>)ds+O(1)Qk+o(1)Γi,jQk.
(64)

In this case, we also have

Γi,jQkfj(x)fi(x)ds=0,
(65)
supθS1supgC0(Qk),g(x)1Qkf(x)dgdθ(x)dx=Γi,jQkfj(x)fi(x)ds+O(1)Qk+o(1)Γi,jQk=O(1)Qk+o(1)Γi,jQk.
(66)

Combining Eqs. (63), (64) and (66), we arrive at the following conclusion: If Qk covers the common boundary Γi,j of a pair of neighboring sub-domains Ωi and Ωj and is contained in ΩiΩj, then

Ik(f)=min(supθS1supgC0(Qk),g(x)1Qkf(x)dgdθ(x)dx,supθS1supθC0(Qk),g(x)1Qkf(x)d2gdθ2(x)dx)=supθS1supgC0(Qk),g(x)1Qkf(x)dgdθ(x)dx=Γi,jQkfj(x)fi(x)ds+O(1)Qk+o(1)Γi,jQk=Γi,jQk(aiaj)x1+(bibj)x2+(cicj)ds+O(1)Qk+o(1)Γi,jQk.
(67)

Next, let us consider Qk that is completely contained in some sub-domain Ωi (Fig. 5(b)). For gC0(Qk), performing 2-D integration by parts again we have

Qkf(x)dgdθ(x)dx=Qkdfdθ(x)g(x)dx,
(68)
Qkf(x)d2gdθ2(x)dx=Qkd2fdθ2(x)g(x)dx.
(69)

Hence,

supgC0(Qk),g(x)1Qkf(x)dgdθ(x)dx=Qkdfdθ(x)dx,
(70)
supgC0(Qk),g(x)1Qkf(x)d2gdθ2(x)dx=Qkd2fdθ2(x)dx.
(71)

Furthermore,

supθS1supgC0(Qk),g(x)1Qkf(x)dgdθ(x)dx=QkmaxθS1dfdθ(x)dx+o(1)Qk,
(72)
supθS1supgC0(Qk),g(x)1Qkf(x)d2gdθ2(x)dx=QkmaxθS1d2fdθ2(x)dx+o(1)Qk.
(73)

Combining Eqs. (72-73), we have

Ik(f)=min(supθS1supgC0(Qk),g(x)1Qkf(x)dgdθ(x)dx,supθS1supgC0(Qk),g(x)1Qkf(x)d2gdθ2(x)dx)=Qkmin(maxθS1dfdθ(x)dx,maxθS1d2fdθ2(x))dx+o(1)Qk=Qkmin{14(2fx122fx22)2+(2fx1x2)2+122fx12+2fx22,(fx1)2+(fx2)2}dx1dx2+o(1)Qk,
(74)

where we have used the fact that f0(x) is a piecewise linear function, 2fx12, 2fx1x2 and 2fx22 are continuous at xΩa\i=1mj>1,jNiΓi,j. Combining (67) and (74), we obtain

k=1MIk(f)=i=1mj>i,jNiΓi,j(aiaj)x1+(bibj)x2+(cicj)ds+Ωamin{14(2fx122fx22)2+(2fx1x2)2+122fx12+2fx22,(fx1)2+(fx2)2}dx1dx2+O(1)(i=1mj>i,jNiΓi,j)(max1kMdiam{Qk})+o(1)Ωa.
(75)

Therefore, we have

HOT2(f)=lim sup(max1kMdiam{Qk})0k=1MIk(f)=i=1mj>i,jNiΓi,j(aiaj)x1+(bibj)x2+(cicj)ds+Ωamin{14(2fx122fx22)2+(2fx1x2)2+122fx12+2fx22,(fx1)2+(fx2)2}dx1dx2.
(76)

which completes the proof of Lemma 4.

For the second term of Eq. (46), because of the zero Lebesgue measure of i=1mj>i,jNiΓi,j, even if

min{14(2fx12(x)2fx22(x))2+(2fx1x2(x))2+122fx12(x)+2fx22(x),(fx1(x))2+(fx2(x))2}
(77)

does not exist for xi=1mj>i,jNiΓi,j, the value of

Ωamin{14(2fx122fx22)2+(2fx1x2)2+122fx12+2fx22,(fx1)2+(fx2)2}dx1dx2
(78)

will remain well defined in the sense of Lebesgue integral.

Now, our main contribution can be presented in the following Theorem.

Theorem 4. Suppose that an object image f0(x) is a piecewise linear function in Ωa as defined by Eq. (36). If h(x) is a candidate image and HOT2(h)=minf=f0+u1HOT2(f) where u1(x) is an arbitrary ambiguity image, then h(x) = f0(x) for x [set membership]Ωa.

Proof: Let h(x) = f0(x) + u(x) for some ambiguity image u . By Lemma 4, we have

HOT2(f)i=1mj>i,jNiΓi,j(aiaj)x1+(bibj)x2+(cicj)ds,forf(x)=f0(x)+u1(x),
(79)

where u1(x) is an arbitrary ambiguity image, and

HOT2(f0)=i=1mj>i,jNiΓi,j(aiaj)x1+(bibj)x2+(cicj)ds.
(80)

Hence,

minf=f0+u1HOT2(f)=i=1mj>i,jNiΓi,j(aiaj)x1+(bibj)x2+(cicj)ds.
(81)

Because HOT2(h)=minf=f0+uHOT2(f), we have

Ωamin{14(2hx122hx22)2+(2hx1x2)2+122hx12+2hx22,(hx1)2+(hx2)2}dx1dx2=0.
(82)

Therefore,

min{14(2hx12(x)2hx22(x))2+(2hx1x2(x))2+122hx12(x)+2hx22(x),(hx1(x))2+(hx2(x))2}=0,forxΩa\i=1mj>i,jNiΓi,j;
(83)

that is,

14(2hx12(x)2hx22(x))2+(2hx1x2(x))2+122hx12(x)+2hx22(x)=0or(hx1(x))2+(hx2(x))2=0,forxΩa\i=1mj>i,jNiΓi,j.
(84)

We assert that

14(2hx12(x)2hx22(x))2+(2hx1x2(x))2+122hx12(x)+2hx22(x)=0,forx=(x1,x2)Ωa\i=1mj>i,jNiΓi,j.
(85)

Otherwise, there exists some x0Ωa\i=1mj>i,jNiΓi,j such that

(14(2hx122hx22)2+(2hx1x2)2+122hx12+2hx22)x=x0>0,
(86)

By continuity, there exists a neighborhood of x0 denoted by Ωx0 such that

Ωx0Ωa\i=1mj>i,jNiΓi,j,
(87)

and

14(2hx12(x)2hx22(x))2+(2hx1x2(x))2+122hx12(x)+2hx22(x)>0,forxΩx0.
(88)

By Eq. (84), we have

(hx1(x))2+(hx2(x))2=0,forxΩx0,
(89)

i.e.,

hx1(x)=hx2(x)=0,forxΩx0.
(90)

Therefore,

2hx12(x)=2hx1x2(x)=2hx22(x)=0,forxΩx0.
(91)

This is in contradiction to Eq. (88). Eq. (85) implies that

2hx12(x)=2hx1x2(x)=2hx22(x)=0,forxΩa\i=1mj>i,jNiΓi,j.
(92)

Because f0(x) is a piecewise linear function in Ωa, it follows from Eq. (92) that

2ux12(x)=2ux1x2(x)=2ux22(x)=0,forxΩa\i=1mj>i,jNiΓi,j.
(93)

Due to the analyticity of u(x) by Theorem 1, we have

2ux12(x)=2ux1x2(x)=2ux22(x)=0,forxΩa.
(94)

Hence,

u(x)=ax1+bx2+c,forxΩa.
(95)

By Theorem 2, we obtain u(x) = 0 and h(x = f0(x).

Generally speaking, we assume that an object image f0(x) is a piecewise n-th (some n ≥ 1) order polynomial function in Ωa; that is, Ωa can be decomposed into finitely many sub-domains {Ωi}i=1m (Fig. 4) such that

f0(x)=fi(x),forxΩi,1im,
(96)

where fi(x) is a n-th order polynomial function, and each sub-domain Ωi is adjacent to its neighboring sub-domains Ωj with piecewise smooth boundaries Γi,j, j [set membership] Ni. We can similarly define the n+1-th order TV for any candidate image f(x) as follows:

HOTn+1(f)=lim sup(max1kMdiam{Qk})0k=1MIk(f),
(97)

where {Qk}k=1M is an arbitrary partition of Ωa, and

Ik(f)=min(supθS1supgC0(Qk),g(x)1Qkf(x)dgdθ(x)dx,supθS1supgC0(Qk),g(x)1Qkf(x)dn+1gdθn+1(x)dx),
(98)

where

dn+1gdθn(x)=dn+1gx,θ(t)dtnt=0.
(99)

We have the following generalized version of Lemma 4.

Lemma 5. Suppose that an object image f0(x) is a piecewise n-th (n ≥ 1) order polynomial function in Ωa, as defined by Eq. (96). For any candidate image f(x) = f0(x) + u(x), we have

HOTn+1(f)=i=1mj>i,jNiΓi,jfi(x)fj(x)ds+Ωamin(maxθS1dfdθ(x),maxθS1dn+1fdθn+1(x))dx1dx2,
(100)

where the second term is a Lebesgue integral.

Proof: We proceed by arguments similar to that for Lemma 4. Again, f(x) = f0(x) + u(x), and u(x) is an analytic function. Let {Qk}k=1M be an arbitrary partition of Ωa. First, let us consider Qk covering the common boundary Γi,j of a pair of neighboring sub-domains Ωi and Ωj and is contained in ΩiΩj (Fig. 5 (a)). As in the proof of Lemma 4, for gC0(Qk) we have

supθS1supgC0(Qk),g(x)1Qkf(x)dgdθ(x)dx=Γi,jQkfj(x)fi(x)ds+O(1)Qk+o(1)Γi,jQk.
(101)

Repeatedly performing 2-D integrations by parts, for gC0(Qk) we have

Qkf(x)dn+1gdθn+1(x)dx=QkΩif(x)dn+1gdθn+1(x)dx+QkΩjf(x)dn+1gdθn+1(x)dx=Γi,jQk(fj(x)fi(x))dngdθn(x)cos<θ,θN>dsQkΩidfdθ(x)dngdθn(x)dxQkΩjdfdθ(x)dngdθn(x)dx=Γi,jQk(fj(x)fi(x))dngdθn(x)cos<θ,θN>dsΓi,jQk(dfjdθ(x)dfidθ(x))dn1gdθn1(x)cos<θ,θN>ds+QkΩid2fdθ2(x)dn1gdθn1(x)dx+QkΩjd2fdθ2(x)dn1gdθn1(x)dx=l=1n(1)lΓi,jQk(dlfjdθl(x)dlfidθl(x))dnlgdθnl(x)cos<θ,θN>ds+(1)n+1(QkΩidn+1fdθn+1(x)g(x)dx+QkΩjdn+1fdθn+1(x)g(x)dx)=l=1n1(1)lΓi,jQk(dlfjdθl(x)dlfidθl(x))dnlgdθnl(x)cos<θ,θN>ds+(1)nΓi,jQk(dnfjdθn(x)dnfidθn(x))g(x)cos<θ,θN>ds+(1)n+1Qkdn+1fdθn+1(x)g(x)dx,
(102)

where θN is the normal vector of the boundary curve Γi,j pointing from Ωj towards Ωi, and <θ,θN > the angle between the vectors θ and θN. For the second term on the right-hand side of Eq. (102), we have

supgC0(Qk),g(x)1(1)nΓi,jQk(dnfjdθn(x)dnfidθn(x))g(x)cos<θ,θN>ds=Γi,jQkdnfjdθn(x)dnfidθn(x)cos<θ,θN>ds,
(103)
supθS1supgC0(Qk),g(x)1(1)nΓi,jQk(dnfjdθn(x)dnfidθn(x))g(x)cos<θ,θN>ds=Γi,jQkmaxθS1(dnfjdθn(x)dnfidθn(x)cos<θ,θN>)ds+o(1)Γi,jQk.
(104)

For the third term on the right-hand side of Eq. (102), we have

(1)n+1Qkdn+1fdθn+1(x)g(x)dx=O(1)Qk.
(105)

We proceed to evaluate supθS1supgC0(Qk),g(x)1Qkf(x)dn+1gdθn+1(x)dx in the following two cases.

Case 1: If there exist some θ [set membership] S1 and some xΓi,jQk such that max0ln1dlfjdθl(x)dlfjdθl(x)0, then for an arbitrary real constant series {Cl}l=1n1 and an arbitrary compact set KQk, there exists a gC0(Qk) satisfying |g(x) and

dnlgdθn1(x)=Cl(1)l(dlfjdθl(x)dlfidθl(x)),for1ln1,xΓi,jK,and
(106)

we have

supθS1supgC0(Qk),g(x)1l=1n1(1)lΓi,jQk(dlfjdθl(x)dlfidθl(x))dn1gdθnl(x)cos<θ,θN>ds=+.
(107)

Combining Eqs. (104), (105) and (107), we obtain

supθS1supgC0(Qk),g(x)1Qkf(x)dn+1gdθn+1(x)dxsupθS1supgC0(Qk),g(x)1l=1n1(1)lΓi,jQk(dlfjdθl(x)dlfidθl(x))dnlgdθnl(x)cos<θ,θN>dssupθS1supgC0(Qk),g(x)1(1)nΓi,jQk(dnfjdθn(x)dnfidθn(x))cos<θ,θN>dssupθS1supgC0(Qk),g(x)1Qkdn+1fdθn+1(x)g(x)dx=+.
(108)

Case 2: Otherwise, dlfjdθl(x)dlfjdθl(x)=0 for xΓi,jQk, θ [set membership] S1, 0 ≤ ln −1. The first term on the right-hand side of Eq. (102) equals 0. Then, we have

supθS1supgC0(Qk),g(x)1Qkf(x)dn+1gdθn+1(x)dx=Γi,jQkmaxθS1(dnfjdθn(x)dnfidθn(x)cos<θ,θN>)ds+O(1)Qk+o(1)Γi,jQk.
(109)

In this case we also have

Γi,jQkfj(x)fi(x)ds=0,
(110)
supθS1supgC0(Qk),g(x)1Qkf(x)dgdθ(x)dx=Γi,jQkfj(x)fi(x)ds+O(1)Qk+o(1)Γi,jQk=O(1)Qk+o(1)Γi,jQk.
(111)

Combining Eqs. (108), (109) and (111), we have the following conclusion: If Qk covers the common boundary Γi,j of a pair of neighboring sub-domains Ωi and Ωj and is contained in ΩiΩj, then

Ik(f)=min(supθS1supgC0(Qk),g(x)1Qkf(x)dgdθ(x)dx,supθS1supgC0(Qk),g(x)1Qkf(x)dn+1gdθn+1(x)dx)=supθS1supgC0(Qk),g(x)1Qkf(x)dgdθ(x)dx=Γi,jQkfj(x)fi(x)ds+O(1)Qk+o(1)Γi,jQk.
(112)

Next, Let us consider Qk that is completely contained in some sub-domain Ωi (Fig. 5 (b)). For gC0(Qk), performing 2-D integration by parts we have

Qkf(x)dgdθ(x)dx=Qkdfdθ(x)g(x)dx,
(113)
Qkf(x)dn+1gdθn+1(x)dx=(1)n+1Qkdn+1fdθn+1(x)g(x)dx.
(114)

Hence,

supgC0(Qk),g(x)1Qkf(x)dgdθ(x)dx=Qkdfdθ(x)dx,
(115)
supgC0(Qk),g(x)1Qkf(x)dn+1gdθn+1(x)dx=Qkdn+1fdθn+1(x)dx.
(116)

Furthermore,

supθS1supgC0(Qk),g(x)1Qkf(x)dgdθ(x)dx=QkmaxθS1dfdθ(x)dx+o(1)Qk,
(117)
supθS1supgC0(Qk),g(x)1Qkf(x)dn+1gdθn+1(x)dx=QkmaxθS1dn+1fdθn+1(x)dx+o(1)Qk.
(118)

Combining Eqs. (117-118), we have

Ik(f)=min(supθS1supgC0(Qk),g(x)1Qkf(x)dgdθ(x)dx,supθS1supgC0(Qk),g(x)1Qkf(x)dn+1gdθn+1(x)dx)=min(QkmaxθS1dfdθ(x)dx,QkmaxθS1dn+1fdθn+1(x)dx)+o(1)Qk=Qkmin(maxθS1dfdθ(x)dx,maxθS1dn+1fdθn+1(x))dx+o(1)Qk.
(119)

In Eqs. (113-119), we have utilized the fact that maxθS1dn+1fdθn+1(x) is continuous in each Ωi(1 ≤ im ). This fact can be derived from the following equation

dn+1fdθn+1(x)=(x1cosφ+x2sinφ)n+1f(x)=l=0n+1Cn+1ln+1fx1lx2n+1l(x)coslφsinn+1lφ,forxΩi,1im,θ=(cosφ,sinφ)S1.
(120)

Combining Eqs. (112) and 119), we obtain

k=1MIk(f)=i=1mj>i,jNiΓi,jfi(x)fj(x)ds+Ωamin(maxθS1dfdθ(x),maxθS1dn+1fdθn+1(x))dx1dx2+O(1)(i=1mj>i,jNiΓi,j)(max1kMdiam{Qk})+o(1)Ωa.
(121)

Therefore, we have

HOT2(f)=lim sup(max1kMdiam{Qk})0k=1MIk(f)=i=1mj>i,jNiΓi,jfi(x)fj(x)ds+Ωamin(maxθS1dfdθ(x),maxθS1dn+1fdθn+1(x))dx1dx2.
(122)

which completes the proof of Lemma 5.

Although it is difficult to have an explicit formula for maxθS1dn+1fdθn+1(x) with xΩa\i=1mj>1,jNiΓi,j in the cases of n ≥ 2, we can still have the following Theorem as a generalized version of Theorem 4.

Theorem 5. Suppose that an object function f0(x) is a piecewise n-th (n ≥1) order polynomial function in Ωa, as defined by Eq. (80). If h(x) is a candidate image and HOTn+1(h)=minf=f0+u1HOTn+1(f) where u1(x) is an arbitrary ambiguity image, then h(x) = f0(x) for x [set membership] Ωa.

Proof: We proceed by arguments similar to that for Theorem 4. Let h(x) = f0(x) + u(x) for some ambiguity image u(x). By Lemma 5, we have

HOTn+1(f)i=1mj>i,jNiΓi,jfi(x)fj(x)ds,forf=f0(x)+u1(x),
(123)

where u1(x) is an arbitrary ambiguity image, and

HOTn+1(f0)=i=1mj>i,jNiΓi,jfi(x)fj(x)ds.
(124)

Hence,

minf=f0+u1HOTn+1(f)=i=1mj>i,jNiΓi,jfi(x)fj(x)ds.
(125)

Because HOTn+1(h)=minf=f0+uHOTn+1(f), we have

ΩaΩamin(maxθS1dhdθ(x),maxθS1dn+1hdθn+1(x))dx1dx2=0.
(126)

Therefore,

min(maxθS1dhdθ(x),maxθS1dn+1hdθn+1(x))=0,forxΩa\i=1mj>i,jNiΓi,j.
(127)

That is,

maxθS1dhdθ(x)=0ormaxθS1dn+1hdθn+1(x)=0,forxΩa\i=1mj>i,jNiΓi,j.
(128)

We assert that

maxθS1dn+1hdθn+1(x)=0,forx=(x1,x2)Ωa\i=1mj>i,jNiΓi,j.
(129)

Otherwise, there exists some x0Ωa\i=1mj>1,jNiΓi,j such that

maxθS1dn+1hdθn+1(x0)>0,
(130)

By continuity, there exists a neighborhoodof x0 denoted by Ωx0 such that

Ωx0Ωa\i=1mj>i,jNiΓi,j,
(131)

and

maxθS1dn+1hdθn+1(x)>0,forxΩx0.
(132)

By Eq. (128), we have

maxθS1dhdθ(x)=0,forxΩx0,
(133)

e.g.,

dhdθ(x)=hx1(x)cosφ+hx2(x)sinφ=0,forxΩx0,θ=(cosφ,sinφ)S1.
(134)

From Eq.(134) we have

hx1(x)=hx2(x)=0,forxΩx0.
(135)

Therefore,

n+1hx1lx2n+1l(x)=0,forxΩx0,0ln+1,
(136)

which leads to

dn+1hdθn+1(x)=l=0n+1Cn+1ln+1hx1lx2n+1l(x)coslφsinn+1lφ=0,forxΩx0,θ=(cosφ,sinφ)S1.
(137)

Eq. (137) means that

maxθS1dn+1hdθn+1(x)=0,forxΩx0.
(138)

This is in contradiction to Eq. (132). Eq. (129) implies that

n+1hx1lx2n+1l(x)=0,forxΩa\i=1mj>i,jNiΓi,j,0ln+1.
(139)

Because f0(x) is a piecewise n-th (n ≥ 1) order polynomial function in Ωa, it follows from Eq. (139) that

n+1hx1lx2n+1l(x)=0,forxΩa\i=1mj>i,jNiΓi,j,0ln+1.
(140)

Due to the analyticity of u(x) by Theorem 1, we have

n+1hx1lx2n+1l(x)=0,forxΩa,0ln+1.
(141)

Hence, u(x) is a n-th (n ≥ 1) order polynomial function in Ωa. By Theorem 2, we obtain u(x) = 0 and h(x) = f0(x).

IV. Numerical Simulation

To verify our theoretical findings described in the preceding section , we developed a HOT minimization based interior tomography algorithm. This iterative algorithm is a modification of our previously reported TV minimization based interior tomography algorithm [17]. Our HOT minimization algorithm consists of two major steps. In the first step, the ordered-subset simultaneous algebraic reconstruction technique (OS-SART) [23] is used to reconstruct a digital image fm,n = f(mΔ,nΔ) based on all the available local projections, where Δ represents the sampling interval, m and n are integers. In the second step, the discrete HOT of fm,n is minimized using the standard steepest descent method. These two steps are iteratively performed in an alternating manner. The major difference between this algorithm and that in [17] lies in the formulas for computing steepest descent directions. Here we assume that 2fx12, 2fx1x2, 2fx22 are continuous at x, and we have

HOT2=Ωa(14(2fx122fx22)2+(2fx1x2)2+122fx12+2fx22)dx1dx2Δ2m,n(14(fm+1,n+fm1,nfm,n+1fm,n1Δ2)2+(fm+1,n+1+fm1,n1fm+1,n1fm1,n+14Δ2)2+12fm+1,n+fm1,n2fm,nΔ2+fm,n+1+fm,n12fm,nΔ2)=m,n((fm+1,n+fm1,nfm,n+1fm,n12)2+(fm+1,n+1+fm1,n1fm+1,n1fm1,n+14)2+fm+1,n+fm1,n+fm,n+1+fm,n14fm,n2)=m,n(Vm,na+Vm,nb),
(142)

where

Vm,na=(fm+1,n+fm1,nfm,n+1fm,n12)2+(fm+1,n+1+fm1,n1fm+1,n1fm1,n+14)2,
(143)
Vm,nb=fm+1,n+fm1,n+fm,n+1+fm,n14fm,n2.
(144)

Therefore,

HOT2fm,n=Dm+1,na,1+Dm1,na,1Dm,n+1a,1Dm,n1a,1+Dm+1,n+1a,2+Dm1,n1a,2Dm+1,n1a,2Dm1,n+1a,2,4Dm,nb+Dm+1,nb+Dm1,nb+Dm,n+1b+Dm,n1b,
(145)

where

Dm,na,1=(fm+1,n+fm1,nfm,n+1fm,n1)4Vm,na,
(146)
Dm,na,2=(fm+1,n+1+fm1,n1fm+1,n1fm1,n+1)16Vm,na,
(147)
Dm,nb={0.5if(fm+1,n+fm1,n+fm,n+1+fm,n14fm,n)00.5if(fm+1,n+fm1,n+fm,n+1+fm,n14fm,n)<0}.
(148)

For the other essential details of our algorithm, the reader may check[17].

A 2D Shepp-Logan phantom was modified to test our algorithm. The modified phantom included a set of piecewise linear ellipses with the parameters listed in Table 1. For each ellipse, the function was first defined inside a compact support Ωe={(x1,x2)R2:x12a12+x22a22<1} with the value (x2ra2+1)μ, then its center was translated to (x10,x20) and rotated by an angle ω. The units for a1, a2 and (x10,x20) were in mm and the unit for ω in degree.

Table 1
Parameters of the 2D modified Shepp-Logan phantom.

In our numerical simulation, we assumed a circular scanning locus of radius 570 mm and a fan-beam imaging geometry. We used an equi-spatial detector array of length 100 mm. 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 consisted of 300 elements, each of which had a 0.33 mm aperture. This configuration covered a circular field of view (FOV) of radius 49.8 mm. For a full scan, we equi-angularly collected 720 projections.

In the aforementioned FOV, we first created a high-resolution phantom image in a 2048×2048 matrix, and numerically generated truncated projections through an ROI using a well-known ray-tracing technique. Then, ROI images were reconstructed in a 256×256 matrix o using our proposed HOT minimization algorithm and a local FBP method from truncated projections after smooth extrapolation into missing data. Representative reconstructed images are in Fig. 6. The typical profiles are in Fig. 7. The convergence curves in terms of truncated projection discrepancy and relative reconstruction error are in Fig. 8. It can be seen in Figs. 6--88 that the reconstructed image using our proposed HOT minimization algorithm is in an excellent agreement with the truth inside the ROI.

Figure 6
Reconstruction of the modified Shepp-Logan phantom with linearly varying shadings over subdomains. (a) The original phantom, (b) the reconstruction with the local FBP (after smooth data extrapolation), and (c) the reconstruction with the proposed HOT ...
Figure 7
Representative profiles along (a) horizontal and (b) vertical lines in Figure 6. The horizontal and vertical axes represent the pixel coordinate and functional value respectively, with the thick line on the horizontal axis indicating the ROI.
Figure 8
Convergence curves of the proposed HOT minimization algorithm. While the horizontal axis represents the iteration index, the vertical axis shows (a) the projection error and (b) reconstruction error (%) respectively.

To evaluate the noise characteristics and demonstrate the stability of our proposed HOT minimization algorithm, Poisson noise [24] was added to the simulated projection data. In this process, 105 photons were emitted from the x-ray source towards each detector aperture. Then, we repeated the above reconstruction procedures and produced the typical images and profiles as shown in Figs. 9 and and10,10, along withthe convergence curves in Fig. 7. These figures demonstrate that the proposed algorithm not only produces satisfactory image accuracy but also suppresses image noise. This is not surprising, since the TV minimization technique was initially proposed to reduce noise [25].

Figure 9
Counterpart of Figure 6, with data noise reflected in the local FBP and HOT minimization reconstructions. The display window is [0.1, 0.4].
Figure 10
Counterpart of Figure 7, with data noise reflected in Figure 9 (b) and (c).

V. Discussions and Conclusion

In Theorem 5, we have claimed that if f0(x) is a piecewise n-th (n ≥ 1) order polynomial function in Ωa, then it can be uniquely determined by minimizing its (n+1)-th order HOT. To apply our theoretical result in practical applications, it is critical to estimate the maximum polynomials order n inside an ROI. In principle, there would be no problem if n is overestimated. Actually, as long as n captures the main features of image variations, even if it is underestimated, we can decompose f0(x) into a n-th polynomial plus a residual error item much smaller than the primary term. Our hypothesis is that the reconstructed results can be written as the n-th polynomial plus a small artifact function. The higher the order n is, the smaller the artifact function becomes. When n is equal to or larger than the real polynomial order, the artifact function would be zero.

In practical applications, high order derivatives can be implemented as high order differences. Given the sensitivity of the discrete high order differencing operation, we suggest that the polynomial order be initially set to 1 ≤ n ≤ 3, and refined as needed. Note that in many cases we do have prior knowledge about images in a particular context. Hence, we can often estimate the polynomial order reliably. For example, in non-destructive testing machinery parts are fairly homogeneous, and the piecewise constant or piecewise linear assumption should work well.

As an important molecular imaging modality, single-photon emission computed tomography (SPECT) is to reconstruct a radioactive source distribution within a patient or animal. Different from the line integral model for x-ray imaging, SPECT projections are modeled as exponentially attenuated Radon transform data [26-28]. Extending x-ray interior tomography results [11-12], we previously proved that accurate and stable interior SPECT reconstruction is feasible from uniformly attenuated local projection data aided by prior knowledge of a sub-region [29]. Naturally, it would be useful to extend the theoretical results reported in this paper to SPECT. That is, it is possible to reconstruct a SPECT ROI accurately only from the uniformly attenuated local projections by minimizing the HOT if the distribution function is piecewise polynomial in the ROI. We are actively working along this direction.

To verify the correctness of our theoretical analysis in Section III, we have developed an iterative algorithm in Section IV. Nevertheless, our numerical implementation has been relatively preliminary. Neither the codes nor the control parameters have been optimized. Currently, we are developing efficient yet robust algorithms for real-world applications. In particular, we may merge the HOT minimization and iterative reconstruction into a more integrated procedure. Hopefully, we will be able to report more in follow-up papers. A theoretically interesting topic is about the convergence of the proposed algorithm. Previously, we studied the convergence of the well-known SART technique [30], and a general Landweber scheme [23] [30]. Although these results are relevant to the proposed algorithm, its convergence analysis may be much more challenging, and will be investigated in the future.

In conclusion, we have extended our TV minimization based interior tomography into a HOT minimization framework. Our work has indicated that an interior ROI can be accurately reconstructed from truncated projection data through the ROI if it can be partitioned into finitely many sub-domains over each of which a polynomial function well represents the image to be reconstructed. This HOT approach may find important applications in biomedical CT such as cardiac perfusion studies, and could inspire development of other imaging modalities as well.

Acknowledgement

We would like to thank anonymous reviewers for their advice and suggestions. This work was partially supported by NSFC (60325101, 60532080, 60628102, 60872078), Key Laboratory of Machine Perception (Ministry of Education) of Peking University, Microsoft Research of Asia, and NIH/NIBIB grants (EB002667, EB004287, EB007288).

References

1. Natterer F. Classics in Applied Mathematics. Society for Industrial and Applied Mathematics; Philadelphia: 2001. The mathematics of computerized tomography.
2. Hamaker C, Smith KT, Solmon DC, Wagner SC. The Divergent beam X-ray transform. Rocky Mountain Journal of Mathematics. 1980;10(1):253–283.
3. Faridani A, Finch DV, Ritman EL, Smith KT. Local tomography II. SIAM J. Appl. Math. 1997;57(4):1095–1127.
4. Faridani A, Ritman EL, Smith KT. Local tomography. SIAM J. Appl. Math. 1992;52:459–484.
5. Ramm AG, Katsevich AI. The Randon transform and local tomography. CRC Press; Boca Raton: 1996.
6. Yu HY, Ye YB, Wang G. Practical cone-beam lambda tomograpy. Medical Physics. 2006;33(10):3640–3646. [PubMed]
7. Katsevich A. Improved cone beam local tomography. Inverse Problems. 2006;22(2):627–643.
8. Katsevich A. Motion compensated local tomography. Inverse Problems. 2008;24(4) Paper ID: 045012 (21pp)
9. Ye YB, Yu HY, Wang G. Cone-beam pseudo-lambda tomography. Inverse Problems. 2007;23(1):203–215.
10. Yu HY, Wang G. A general formula for fan-beam lambda tomography. International Journal of Biomedical Imaging. 2006 2006, Article ID: 10427, 9 pages. [PMC free article] [PubMed]
11. Ye YB, Yu HY, Wei YC, 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]
12. Kudo H, Courdurier M, Noo F, Defrise M. Tiny a priori knowledge solves the interior problem in computed tomography. Phys. Med. Biol. 2008;53(9):2207–2231. [PubMed]
13. 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]
14. Candes EJ, Romberg J, Tao T. Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information. IEEETransactions on Information Theory. 2006;52(2):489–509.
15. Donoho DL. Compressed sensing. IEEE Transactions on Information Theory. 2006;52(4):1289–1306.
16. 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, Editor. 2008.
17. Yu HY, Wang G. Compressed sensing based Interior tomography. Phys Med Biol. 2009;54(9):2791–2805. [PMC free article] [PubMed]
18. Yu HY, Jiang JS, Jiang M, Wang G. Further analysis on compressed sensing based interior tomography. Phys Med Biol. 2009;54(18):N425–N432. [PMC free article] [PubMed]
19. Han WM, Yu HY, Wang G. A total variation minimization theorem for compressed sensing based tomography. International Journal of Biomedical Imaging. 2009 2009: Article ID: 125871, 3 pages. [PMC free article] [PubMed]
20. Courdurier M. PhD. Dissertation. University of Washington-Seattle; 2007. Restricted Measurements for the X-ray Transform.
21. Gel'fand IM, Graev MI. Crofton function and inversion formulas in real integral geometry. Functional Analysis and its Applications. 1991;25(1):1–5.
22. 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]
23. Wang G, Jiang M. Ordered-subset simultaneous algebraic reconstruction techniques (OS-SART). Journal of X-ray Science and Technology. 2004;12(3):169–177.
24. Yu HY, Ye YB, Zhao SY, Wang G. Local ROI reconstruction via generalized FBP and BPF algorithms along more flexible curves. International Journal of Biomedical Imaging. 2006;2006(2) Article ID: 14989,7 pages. [PMC free article] [PubMed]
25. Rudin LI, Osher S, Fatemi E. Nonlinear total variation based noise removal algorithms. Physica D. 1992;60(1-4):259–268.
26. Rullgard H. An explicit inversion formula for the exponential Radon transform using data from 180~. Ark. Mat. 2004;42:353–362.
27. Rullgard H. Stability of the inverse problem for the attenuated Radon transform with 180 degrees data. Inverse Problems. 2004;20(3):781–797.
28. Noo F, Defrise M, Pack JD, Clackdoyle R. Image reconstruction from truncated data in single-photon emission computed tomography with uniform attenuation. Inverse Problems. 2007;23(2):645–667.
29. Yu HY, Yang JS, Jiang M, Wang G. Interior SPECT-exact and stable ROI reconstruction from uniformly attenuated local projections. International Journal for Numerical Methods in Engineering. 2009;25(6):693–710. [PMC free article] [PubMed]
30. Jiang M, Wang G. Convergence studies on iterative algorithms for image reconstruction. IEEE Trans. Med. Imaging. 2003;22:569–579. [PubMed]