- Email: [email protected]
Fast and exact 2d image reconstruction based on Hakopian interpolation
Applied Numerical Mathematics 121 (2017) 185–197
Contents lists available at ScienceDirect
Applied Numerical Mathematics www.elsevier.com/locate/apnum
Fast and exact 2d image reconstruction based on Hakopian interpolation ✩ Xuenan Sun a,∗ , Xuezhang Liang b a b
School of Mathematics and Statistics, Northeast Normal University, Changchun, 130024, China Institute of Mathematics, Jilin University, Changchun, 130022, China
a r t i c l e
i n f o
Article history: Received 8 September 2016 Received in revised form 20 May 2017 Accepted 13 July 2017 Available online 21 July 2017
a b s t r a c t A new algorithm for the reconstruction of two-dimensional (2D) images from projections is given. The algorithm is based on the expansions of Hakopian interpolation polynomial into Chebyshev–Fourier series. The computer simulation experiments show that the new algorithm is effective. © 2017 IMACS. Published by Elsevier B.V. All rights reserved.
Keywords: Hakopian interpolation Chebyshev–Fourier series Image reconstruction Radon transformation
1. Introduction The fundamental problem of the computed tomography (CT) is to reconstruct an image from its Radon projections. More precisely, we consider the case of a two dimensional image, which is described by a density function f (x, y ) defined on the unit disk D = {(x, y )|x2 + y 2 ≤ 1} of the R2 plane. Let f be any absolutely integrable function on D. Then a Radon projection of f is a line integral as follows,
Rθ ( f ; t ) :=
f (x, y )dxdy , ∀ 0 ≤ θ ≤ 2π
and
−1 ≤ t ≤ 1,
I (θ,t )
where I (θ, t ) = {(x, y ) ∈ R2 | x cos θ + y sin θ = t } ∩ D is a line segment inside D. By the Radon transform f → {Rθ ( f ; t )}, f is associated to a family of univariate functions of t, parameterized by θ . The essential problem of image reconstruction is to recover f from its Radon projections. It is well known that the set of Radon projections {Rθ ( f ; t ) | |t | ≤ 1, 0 ≤ θ ≤ 2π } determines f completely. This was proved by Radon [14] and John [7] for differentiable functions. Furthermore, if f ∈ L 1 (R2 ) has compact support in D, then f is uniquely determined by any infinite set of Radon projections [18]. However, only a finite number of projections can be measured in practice. Thus, the main problem of image reconstruction is to obtain an effective algorithm that produces a good approximation to f from a large collection of its Radon projections. For the practical reconstruction of two-dimensional images from their measured projections, one most frequently applies the method of filtered backprojection (FBP) and its modifications [15]. These algorithms can be viewed as computer imple-
✩ This work was partially supported by the National Science Foundation of Jilin Province (No. 20170101044JC), the Education Department of Jilin Province (No. JJKH20170904KJ) and by the National Natural Science Foundation of China (No. 11171056). Corresponding author. E-mail address: [email protected] (X. Sun).
*
http://dx.doi.org/10.1016/j.apnum.2017.07.002 0168-9274/© 2017 IMACS. Published by Elsevier B.V. All rights reserved.
186
X. Sun, X. Liang / Applied Numerical Mathematics 121 (2017) 185–197
mentation of an inversion of the Radon transform through the central slice theorem [1]. They use a discrete convolution of the projections with a filter function either in the spatial domain or in the frequency domain, plus backprojection. Another class of solutions utilizes a decomposition of the Radon transform into circular harmonics and certain polynomials [4,3]. The Chebyshev polynomials of the second kind play a key role in this context [2,20,12]. In [3], an algorithm was described based on the decomposition of the projections into Chebyshev polynomials of the second kind, which uses the backprojection technique. Consequently, the method can be considered as a Chebyshev-domain filtered backprojection (CDFBP). The method combines the exactness of the polynomial approach with the speed of FBP. In [19], an algorithm called OPED was provided, as it was based on orthogonal polynomial expansion on the unit disk. The numerical tests showed that the algorithm was fast and stable, and it produced high quality images [21]. In comparison with FBP, these algorithms are the kind which have slightly higher mathematical complexity and less intuitive expression. In tomography, the data often consists of values of linear integrals over segments. In many situations, a table of mean values of a function of d variables on (d − 1)-dimensional hyperplane is considered to be the most natural type of data for multivariate functions. Hakopian’s famous interpolation formula [5] is an important reason to take this approach. Hakopian proved that for any given n distinct points X 1 , · · · X n on the boundary of a disk D, the set of integrals of f over all the linear segments [ X i , X j ] determines uniquely every polynomial f of degree n − 2. A Lagrange representation of the bivariate Hakopian interpolation was given in [8]. In fact, the first method of recovering a polynomial of degree n from its integrals over chords defined by equally spaced points { X i } on the boundary of the unit disk was described by Marr in [12]. Interpolation methods based on integrals over chords can be used for approximate reconstruction of functions from their Radon transforms. Because of the importance of such recovery methods for applications in tomography, they have been intensively studied (see [12,11]). In this paper, we develop an efficient and simple reconstruction algorithm based on the Hakopian interpolation method, which uses the Chebyshev–Fourier analysis on the unit disk. Therefore, the new algorithm can be implemented by fast discrete sine transform and an interpolation step. Indeed, it combines the exactness of the polynomial approach with the speed of FBP. The rest of this paper is organized as follows: In Section 2, we discuss the background of Hakopian interpolation and orthogonal polynomials. Section 3 is devoted to the Chebyshev–Fourier orthogonal expansions of Hakopian interpolation are studied. Readers who are not interested in the mathematical details may skip most of this section until its end, where a step-by-step outline of the algorithm is given. Reconstructions are made of the Shepp and Logan phantom in Section 4. The numerical example is given to show that the algorithm is effective. 2. Preliminaries In this section, we present some known results, which are useful later. To begin with, we denote by n (R2 ) the following set of all real algebraic polynomials in two variables of total degree n,
n (R2 ) :=
C j ,k x j yk
x, y , C j ,k ∈ R .
j +k≤n
We shall consider real functions f on the plane R2 . Also, a point x of R2 is denoted by x = (x, y ) . Throughout this paper, the functions f which are to be approximated are assumed to have their support in the unit disk D := {(x, y ) ∈ R2 |x2 + y 2 ≤ 1}. Denote by C (D) the space of all continuous functions defined on D, and define
E m ( f ) :=
inf
max | f (x) − p (x)|, ∀ f ∈ C (D).
p ∈m (R2 ) x∈D
For a given natural number n ≥ 2, we choose the nodes of Hakopian interpolation on the disk which are defined as follows:
xk = (cos
2kπ n
, sin
2kπ n
) ,
k = 1, 2, · · · , n .
Introduce the following integral functional:
1 f (xk + t (xl − xk ))dt , ∀ f ∈ C (D),
f = [xk ,xl ]
0
and denote by λ(k, l) the Radon projections of two univariate function f (x), that is
λ(k, l) = [xk ,xl ]
f,
1 ≤ k < l ≤ n.
X. Sun, X. Liang / Applied Numerical Mathematics 121 (2017) 185–197
187
Moreover, let
U m−1 (t ) = T m (t ),
T m (t ) = cos m arccos t , R m−1 (t ) =
T m (t ) + T m−1 (t ) t+1
S m−1 (t ) =
,
T m (t ) − T m−1 (t )
x = (x, y ) = (r cos θ, r sin θ) , and t = r cos (
t−1 k+l n
,
π − θ).
Next, we recall the following known results in [8,9]. Lemma 2.1. For any f (x) ∈ C (R2 ), there exists a unique polynomial P (x) = H n ( f ; x) ∈ n−2 (R2 ), such that
f − P = 0,
0 ≤ k < l ≤ n.
[xk ,xl ]
Here P (x) = H n ( f ; x) denotes the Hakopian interpolation polynomial of f (x, y ) ∈ C (D) with respect to the n nodes. Furthermore, we have the following representation:
H n ( f ; x) =
L k,l (x) ·
1≤k
f,
(1)
[xk ,xl ]
where (1) If n is even and let n = 2m, then
L k,l (x) =
sin2 l−nk π T m (t )
m2 (t − cos l−nk π )
L k,l (x) = −
2U m−1 (t ) −
sin2 l−nk π U m−1 (t ) m2 (t − cos l−nk π )
2T m (t ) +
T m (t ) t − cos l−nk π
when l − k is odd,
,
(1 − t 2 )U m−1 (t ) m2 (t − cos l−nk π )
(2)
,
when l − k is even.
(3)
(2) If n is odd and let n = 2m − 1, then
L k,l (x) =
l−k 2 2(sin 2m −1 π ) R m−1 (t )
l−k (2m − 1)2 (t − cos 2m −1 π )
2(1 + t ) R m −1 (t )
− L k,l (x) =
l−k 2 2(sin 2m −1 π ) S m−1 (t )
l−k (2m − 1)2 (t − cos 2m −1 π )
l−k (1 + cos 2m −1 π ) R m−1 (t )
,
when l − k is odd,
,
when l − k is even.
l−k t − cos 2m −1 π
2(1 − t ) S m −1 (t )
−
l−k (1 − cos 2m −1 π ) S m−1 (t ) l−k t − cos 2m −1 π
Lemma 2.2. There exists an absolute constant A, such that for arbitrary x ∈ D, f (x) ∈ C (D) and n ≥ 2, we have
| f (x) − H n ( f ; x)| ≤ A E n−2 ( f ) · n log n. Note that, the Radon projection is easily computed in the special case when f has a “plane wave” structure. More precisely, let un (t ) denote the Chebyshev polynomial of the second kind, that is
1 un (t ) = √ D n (arccos t ),
π
where
D n (θ) =
sin(n + 1)θ sin θ
,
n = 0, 1 , · · · .
188
X. Sun, X. Liang / Applied Numerical Mathematics 121 (2017) 185–197
These classical polynomials constitute a complete orthonormal system U = {un−1 (t )}n∞=1 in l2ω (−1, 1) with the weight
√
ω(t ) = 2 1 − t 2 .
After substitution t = x · ξ , for x ∈ D and ξ ∈ S 1 (S 1 is the set of unit vector). Then un (x · ξ ) generate a family of complete orthonormal systems of ridge polynomials in L 2 (D) as follows:
U S 1 = {{un−1 (x · ξ )}n∞=1 }ξ ∈ S 1 . The fundamental properties of this system are expressed by the following [13,16]: (1) Orthogonality relation
un (x · ξ ) p (x)dx = 0,
∀ p (x) ∈ n−1 (R2 ),
n = 1, 2, · · · ,
∀ξ ∈ S 1 ,
D
or un (x · ξ )⊥n−1 (R2 ) in L 2 (D). (2) Inner products of Chebyshev ridge polynomials
un (x · ξ )un (x · η)dx = D
√
2π
un (1)
√ =
π
n+1
un (ξ · η),
∀ξ, η ∈ S 1 .
(4)
τn± denote the subspace of trigonometric polynomials a(·) on S 1 , of degree n and satisfying a(−ξ ) ≡
Furthermore, let (−1)n a(ξ ). Then
1
un (ξ · η)
π un (ξ · η)a(η)dη = a(ξ ), a(ξ ) ∈ τn± , ξ ∈ S 1 ,
S1
and in particular
1
2
√
un (ξ1 · η)un (η · ξ2 )dη = un (ξ1 · ξ2 ),
π
ξ 1 , ξ2 ∈ S 1 .
S1
(3) Integral representation Given a function f (x) ∈ L 2 (D), consider the following Chebyshev–Fourier coefficients, depending on ξ ∈ S 1 as a parameter:
an ( f , ξ ) =
f (y)un (y · ξ )dy,
n = 0, 1 , · · · , ξ ∈ S 1 .
D
The following relation represents the integral form of Chebyshev ridge polynomial Fourier series, which is in fact Radon inversion formula for f (x) expressed via Chebyshev ridge polynomials: ∞ 1 f (x) = n an−1 ( f , ξ )un−1 (x · ξ )dξ. 2π L 2 (D)
n =1
S1
3. Chebyshev–Fourier expansions of Hakopian interpolation Our reconstruction algorithm will be based on this Chebyshev–Fourier orthogonal expansion. Radon–Fourier analysis via Chebyshev ridge polynomials is crucial in the derivation for the algorithm. In order to get the expansions of Hakopian interpolation, we need the following lemma [10]. Lemma 3.1. For every natural number n ≥ 2 and 1 ≤ j < k ≤ n, we have
π 0
π 0
cos2 n2 θ cos θ − cos
k− j n
π
sin2 n2 θ cos θ − cos
k− j n
π
d θ = 0,
if k − j is odd,
(5)
d θ = 0,
if k − j is even.
(6)
X. Sun, X. Liang / Applied Numerical Mathematics 121 (2017) 185–197
189
By making use of the orthogonality of Jacobi polynomials with respect to the weight function (1 − x)±1/2 (1 + x)±1/2 , the proof of Lemma 3.1 is not difficult. Here we omit the proof of it [10]. For simplicity, we set
k+l
η = cos
n
t = r cos
π , sin
l+k n
k+l n
= (cos β, sin β) ∈ S 1 ,
π
c = ωk ,l (cos
π −θ
l−k n
π ).
= r cos θ cos β + r sin θ sin β
= x cos β + y sin β = x · η, and let
ϕ (t ) :=
ωk ,l (t ) ωk ,l (cos l−nk π )
ωk ,l (t )
=
c
= Lk,l (x).
Now, from the orthogonality of the Chebyshev polynomials un (t ), it follows that
ϕ (t ) =
n −2
c j u j (t ),
(7)
j =0
with the coefficients given by
cj =
1
1 − t 2 ϕ (t )u j (t )dt
−1
=
1 c
1
1 − t 2 ωk ,l (t )u j (t )dt
−1
π
1
= √ c
π
sin θ ωk ,l (cos θ) sin( j + 1)θ dθ
0
j+1
= √ c
π
ωk,l (cos θ) cos( j + 1)θ dθ.
π 0
In the following, we compute the coefficient for the Hakopian interpolation expansion in two cases, respectively. Case 1: If n is even, let n = 2m. It is easy to show that
ωk,l (t ) =
⎧ 2 t) ⎨c 0 cos (m arccos l−k t −cos n π 2 t) ⎩c 0 sin (m arccos t −cos l−nk π
c = ωk ,l (cos
l−k n
π ) = c0
when l − k is odd, when l − k is even. m2
sin2 l−nk
π
,
where c 0 is a constant. Set α = l−nk π . From Lemma 3.1, we obtain that
( j + 1)c 0 cj = √ c
π
π
0
=
( j + 1)c 0 √ 2c
2π
π
0
=
( j + 1)c 0 √ 4c
2π
π
0
cos2 mθ cos θ − cos α cos2 mθ cos θ − cos α
cos( j + 1)θ dθ
(cos( j + 1)θ − cos( j + 1)α )dθ
sin( j + 1) θ −2 α sin θ −2 α
·
sin( j + 1) θ +2 α sin θ +2 α
dθ
190
X. Sun, X. Liang / Applied Numerical Mathematics 121 (2017) 185–197
( j + 1)c 0 = √
π
sin( j + 1)(θ − α ) sin(θ − α )
π
2c
0
=
√ ( j + 1)c 0 π
sin( j + 1)α
·
c
sin α
·
sin( j + 1)θ sin θ
dθ
,
and consequently,
ϕ (t ) =
n −2
√ ( j + 1)c 0 π
·
c
j =0
sin( j + 1)α sin α
· u j (t ) .
If let t s = cos ns π , then
ϕ (t ) =
n −2 ( j + 1)π c 0
c
j =0
u j (tl−k )u j (t ) .
Therefore,
L k,l (x)u i (x · ξ )dx =
n −2 ( j + 1)π c 0
c
j =0
D
= = Then, since by
c0 c
=
4 sin2 α , n2
(i + 1)π c 0 c
c0 π
√
c
π
u j (x · η)u i (x · ξ )dx
u j (tl−k )
√ u i (tl−k )
D
π
i+1
u i (ξ · η)
u i (tl−k )u i (ξ · η).
we can further have
L k,l (x)u i (x · ξ )dx = 4π
sin l−nk π sin(i + 1) l−nk π n2
D
u i (ξ · η).
With the aid of (6), similarly, we can obtain the same result when l − k is even. Case 2: If n is odd, let n = 2m − 1, then
ωk,l (t ) =
⎧ 2 ⎨c 1 (cos(m arccos t )+cos((ml−−k 1) arccos t )) (t +1)(t −cos n π ) ⎩c 1 (cos(m arccos t )−cos((ml−−k 1) arccos t ))2 (t −1)(t −cos n π )
when l − k is odd, when l − k is even,
and
c = ωk ,l (cos
l−k n
π ) = c1
(2m − 1)2 2 sin2 l−nk π
,
where c 1 is a constant. Analogously to the Case 1, we will only derive the case that l − k is odd. From Lemma 3.1 we obtain that
( j + 1)c 1 cj = √ c
π
π
0
=
( j + 1)c 1 √ c
2π
π
0
=
( j + 1)c 1 √ c
π
π
0
=
2( j + 1)π c 1 c
(cos mθ + cos((m − 1)θ))2 cos( j + 1)θ dθ (cos θ + 1)(cos θ − cos α ) (cos 2m2−1 )2 cos( j + 1)θ dθ cos θ − cos α sin( j + 1)(θ − α ) sin(θ − α ) u j (tl−k ),
·
sin( j + 1)θ sin θ
dθ
X. Sun, X. Liang / Applied Numerical Mathematics 121 (2017) 185–197
and therefore, we can get the orthogonal expansions of
ϕ (t ) =
n −2 2( j + 1)π c 1
c
j =0
191
ϕ (t ):
u j (tl−k )u j (t ).
(8)
It follows that
L k,l (x)u i (x · ξ )dx =
n −2
2( j + 1)π c 1 c
j =0
D
3 2
2(π ) c 1
= and because of
c1 c
=
2 sin2 α , n2
c
u j (x · η)u i (x · ξ )dx
u j (tl−k ) D
u i (tl−k )u i (ξ · η),
we have
L k,l (x)u i (x · ξ )dx = 4π
sin l−nk π sin(i + 1) l−nk π n2
D
u i (ξ · η).
Combining Case 1 and Case 2, for arbitrary natural number n, we obtain that
L k,l (x)u i (x · ξ )dx = 4π
sin l−nk π sin(i + 1) l−nk π n2
D
u i (ξ · η).
(9)
Therefore, with the aim of (9), we have the Chebyshev–Fourier coefficient of H n ( f ; x):
⎛
⎜ L k,l (x) ⎝
ai ( H n ( f ; x), ξ ) =
1≤k
D
⎛
=
⎝
1≤k
=
⎞
⎟
f ⎠ u i (x · ξ )dx
[xk ,xl ]
⎞
L k,l (x)u i (x · ξ )dx⎠ ·
D
4π
f
[xk ,xl ]
sin l−nk π sin(i + 1) l−nk π n2
1≤k
u i (ξ · η)
f
[xk ,xl ]
Especially when i ≥ n − 1, ai = 0. Then,
H n ( f ; x) =
=
=
n −2 1
2π 2 n2 4
ω =1 n −2
π
n2
S1
u ω−1 (ξ · η)u ω−1 (x · ξ )dξ
ω
ω =1
√
aω−1 ( H n ( f ; x), ξ )u ω−1 (x · ξ )dξ
ω
1≤k
S1
n −2
ω =1
ω
sin
l−k
1≤k
n
π sin ω
l−k n
sin
l−k n
π sin ω
l−k n
π
f [xk ,xl ]
π u ω−1 (t )
f,
[xk ,xl ]
where t = x cos( l+nk π ) + y sin( l+nk π ). Without loss of generality, we suppose that n is even and let n = 2m (the similar result can be obtained if n is odd). Then,
H n ( f ; x) =
n −2 4
n2
ω =1
ω
1
+
2
(10)
192
X. Sun, X. Liang / Applied Numerical Mathematics 121 (2017) 185–197
where
1
=
m m
t = x cos
2
2i − 1 2m
j =1 i =1
sin
=
2j −1 2m
m m −1
sin
j =1 i =1
t = x cos
j m
i
π sin ω
2i − 1 2m
π + y sin
π u ω−1 (t )
2j −1 2m
j = 1, 2, · · · , m ,
π sin ω π u ω−1 (t ) m
π + y sin
(11)
[x j−i ,xi + j−1 ]
π
i
m
f
j m
f
(12)
[x j−i ,x j+i ]
j = 1, 2, · · · , m .
π
k In addition, the index k of x in (11) and (12) will probably be negative or exceed n. In these two cases, we define k±n x =x . That is, in 1 , if k = j − i < 0, then we set k = j − i + n; if l = i + j − 1 > n, then we set l = i + j − 1 − n. In 2, if k = j − i < 0, then we set k = j − i + n; if l = i + j > n, then we set l = i + j − n. For simplicity, we introduce the following notations k
(1 )
c ω, j =
m
2i − 1
sin
n
i =1
(2 )
c ω, j =
m −1 i =1
sin
2i n
π sin ω
π sin ω
2i n
2i − 1 n
π
f, [x j−i ,xi + j−1 ]
π
f. [x j−i ,xi + j ]
Then, these findings can be translated into a simple reconstruction algorithm that consists of the following steps: (1)
(2)
Step 1. Computing c ω, j and c ω, j , respectively, by projections. Step 2. Computing the summation of λ1 (t , j ) = Step 3. Computing the summation H n ( f ; x) =
n −2
ωcω(1,)j u ω−1 (t ) and λ2 (t , j ) =
ω=1 √ m 4 π n2
j =1
n −2
ω=1
ωcω(2,)j u ω−1 (t ).
{λ1 (t , j ) + λ2 (t , j )}.
In reconstruction of step 1, the c ω, j and c ω, j can be directly calculated by a discrete sine transform of sin 2in−1 π × 2i [x j −i ,xi + j −1 ] f and sin n π [x j −i ,xi + j ] f . Very efficient algorithms requiring only on the order of (m + 1) · log(m + 1) operations can be applied for this purpose [3]. Since the sine transform is its own inverse, we can also use it to speed up the computation of λ1 (tk , j ) and λ2 (tk , j ) (step 2), where tk = cos knπ , k = 1, 2, · · · , n − 1. Values at other points t, which are needed for the reconstruction, are obtained by linear interpolation. (1)
(2)
4. Examples of the algorithm To demonstrate the practical usefulness of the algorithm, we apply it to the reconstruction of the Shepp–Logan head phantom [17] using line integrals. This phantom head is well known because of its use in testing the accuracy of reconstruction algorithms with regards to their ability to reconstruct cross sections of the human head with X-ray tomography. A particular advantage of this phantom is that it is constructed from the superposition of 10 ellipses of homogeneous “density” and thus, based on the linearity of the Radon Transform, the projections of this construction can be specified analytically. This phantom has highly singularities as the image contains jumps at the boundary of every ellipse in the image, include the one on the boundary. The function that represents the image is a step function, which is not continuous at the boundary of each of the ellipses. We do the reconstruction on a 256 × 256 grid using 256 projection angles. The original image and the reconstruction are given in Figs. 1 and 2. In addition to the Shepp–Logan phantom, we reconstruct another image on a 280 × 280 grid using 256 projection angles. The original function named “intersected ellipse” is defined by ourselves. The original image and the reconstruction are shown in Figs. 3 and 4. Furthermore, a more quantitative way of evaluating pictures is the following. Select a column of pixels, which is such that it goes through a number of interesting features in the original. For example, in Shepp–Logan head phantom the 80th
X. Sun, X. Liang / Applied Numerical Mathematics 121 (2017) 185–197
193
Fig. 1. Shepp–Logan head phantom.
Fig. 2. Reconstruction image by this algorithm.
Fig. 3. The original intersected ellipse.
of the 256 columns goes through the above ellipse, the left rotated ellipse, and the left horizontal ellipse. One way to evaluate the quality of reconstruction is compare the graphs of the 256 pixel densities for this column in the original and the reconstruction. Fig. 5 shows the values of the original and the reconstruction along the 80th column. At the same time, we introduce two picture distance measures to measure the closeness of the reconstruction to the original [6]. We use t u , v and r u , v to denote the densities of the vth pixel of the uth row of the test phantom and reconstruction respectively, and t to denote the average of the densities in the test phantom. We assume that both pictures are l × l. The measures are defined in the following:
194
X. Sun, X. Liang / Applied Numerical Mathematics 121 (2017) 185–197
Fig. 4. Reconstruction image by this algorithm.
Fig. 5. Line plots of the 80th column of phantom and reconstruction.
⎡
l l
(t u , v − ru , v ) ⎢ ⎢ u =1 v =1 l l ⎣ (t u , v − t )2
2
d=⎢
⎤ 12 ⎥ ⎥ ⎥ , ⎦
u =1 v =1
l
r=
l
u =1 v =1
|t u , v − ru , v |
l l u =1 v =1
. |t u , v |
These two measures emphasize different aspects of picture quality. The first one, d, is a normalized root mean square distance measure. A large difference in a few places causes the value of d to be large. The second one, r, is a normalized mean absolute distance measure. As opposed to d, it emphasizes the importance of a lot of small errors rather than of a few large errors [6]. Picture distance measures for reconstructions from different projection data are shown in Tables 1, 2, 3 and 4. Table 1 and Table 2 show the computational results of the Shepp–Logan phantom. The computational results of the “intersected ellipse” are given in Table 3 and Table 4. From these tables we can see that the reconstruction obtained from Hakopian interpolation algorithm compares well with reconstructions from FBP method, and the reconstructions from the Hakopian algorithm and CD-FBP method are almost the same. Moreover, the calculation time for the reconstructions (256 × 256 Shepp–Logan phantom) is given in Table 5. We can use it as a reference to the computational complexity of the three methods. Table 5 shows that the computational time of FBP method is shorter than the other ones. Meanwhile, the computational time of the Hakopian algorithm can be acceptable.
X. Sun, X. Liang / Applied Numerical Mathematics 121 (2017) 185–197
195
Table 1 Picture distance measures for Shepp–Logan phantom. Number of projection angles
120 128 256
CD-FBP algorithm
Hakopian algorithm
d
r
d
r
0.197571 0.191253 0.144232
0.100863 0.095916 0.051576
0.203349 0.197040 0.144907
0.106673 0.104356 0.068250
Table 2 Picture distance measures for Shepp–Logan phantom. Number of projection angles
120 128 256
FBP algorithm
Hakopian algorithm
d
r
d
r
0.406432 0.395613 0.360563
0.370060 0.359921 0.352121
0.203349 0.197040 0.144907
0.106673 0.104356 0.068250
Table 3 Picture distance measures for “intersected ellipse”. Number of projection angles
64 128 256
CD-FBP algorithm
Hakopian algorithm
d
r
d
r
0.130408 0.093564 0.071143
0.042846 0.027616 0.018604
0.135765 0.097128 0.073810
0.047598 0.030880 0.021450
d
r
d
r
0.144536 0.100776 0.097522
0.048256 0.030940 0.023176
0.135765 0.097128 0.073810
0.047598 0.030880 0.021450
Table 4 Picture distance measures for “intersected ellipse”. Number of projection angles
64 128 256
FBP algorithm
Hakopian algorithm
Table 5 The computational time of three methods. Algorithm
128 projection angles
256 projection angles
512 projection angles
FBP CD-FBP Hakopian
about 2 s about 9 s about 11 s
about 2 s about 12 s about 39 s
about 3 s about 13 s about 93 s
Sometimes, people are only interested in a local region of the reconstructed image. So we reconstruct the left and right rotated ellipse of the Shepp–Logan from the local projection data. The projections are line integrals within the region of interest. The others are zero outside of the local region. The left local region is 45 pixels in radius in a 200 × 200 image. The right local region is 35 pixels in radius. The local reconstructions are shown in Figs. 6 and 7. The figures show that the algorithm based on Hakopian interpolation is able to reconstruct the region of interest. All reconstructions were done on a 1.80 GHz Intel(R) Core(TM) i5-3337U CPU personal computer. The programming language is Visual C++6.0, and the operating system is Windows 7 Professional. We used the FFT for discrete sine transform in the package FFTW (http://www.fftw.ort/). 5. Conclusions A fast and easy to implement algorithm for 2D image reconstruction has been developed. Mathematically, basis is a decomposition of Hakopian interpolation into Chebyshev polynomials. To the best of our knowledge the described implementation in form of this algorithm is new, although there is a loose relationship to other recent developments ([3,2,20]). The main advantage of the algorithm as compared with the FBP method is its exactness. We introduce a fast implementation of Hakopian interpolation algorithm by using FFT and a linear interpolation step. We hope that our algorithm will be useful for both practical and theoretical problems.
196
X. Sun, X. Liang / Applied Numerical Mathematics 121 (2017) 185–197
Fig. 6. The local reconstruction of left rotated ellipse.
Fig. 7. The local reconstruction of right rotated ellipse.
Acknowledgements The authors thank the anonymous reviewers for their thorough reading of the paper and their constructive comments. Appendix A In this part, we need to explain the relation between line integrals and projection angles. As well as we know, every straight line in the plane will be characterized as the locus of points satisfying
x cos θ + y sin θ = t ,
(A.1)
for some real distance parameter t and angle θ . For convenience, we allow t to assume negative values in (A.1), so that every straight line has exactly two distinct parametric representations, (t , θ) and (t , θ ), with t = −t, θ = θ + π (mod 2π ). We mean the set of parameters by C,
C = {(t , θ)| − 1 ≤ t ≤, 0 ≤ θ ≤ 2π }. Then every chord, or straight line intersecting D, is represented by exactly two points on C. Now, we return to the algorithm described in this paper. We consider an interpolation nodes set, X n = {xk |xk = (cos 2(k−n1)π , sin 2(k−n1)π ) , k = 1, 2, · · · , n}, of n points spaced uniformly around the circumference of D. We can suppose that these points are numbered from 1 to n, in counterclockwise order, and for simplicity assume that point 1 lies along the positive x-axis. Then, xk are the coordinates of the k-th point. Therefore, the n(n − 1)/2 possible chords joining pairs of points in X n correspond to a mesh of n(n − 1) points on C by (A.1). To describe these, we define
tl = cos(π l/n), l = 1, 2, · · · , n − 1,
(A.2)
and
θk,l =
π n
(2k + l − 2), k = 1, 2, · · · , n, l = 1, 2, · · · , n − 1.
(A.3)
X. Sun, X. Liang / Applied Numerical Mathematics 121 (2017) 185–197
197
Then, the point (tl , θk,l ) on C corresponds to the chord joining the two points in X n numbered k1 and k2 , where
k1 = k < k2 = k + l,
if k + l = n,
k1 = k + l − n < k2 = k,
if k + l > n.
Thus, we connect the number of line integrals with the number of projection angels. This situation is the same as the discrete case in [12]. References [1] H.H. Barrett, The Radon transform and its applications, in: E. Wolf (Ed.), Progress in Optics XXI, Elsevier, Amsterdam, 1984. [2] B. Bojanov, I.K. Georgieva, Interpolation by bivariate polynomials based on Radon projections, Stud. Math. 162 (2004) 141–160. [3] Thomas Bortfeld, Uwe Oelfke, Fast and exact 2D image reconstruction by means of Chebyshev decomposition and backprojection, Phys. Med. Biol. 44 (1999) 1105–1120. [4] A.M. Cormack, Representation of a function by its line integrals with some radiological applications, J. Appl. Phys. 34 (1963) 2722–2727. [5] H.A. Hakopian, Multivariate divided difference and multivariate interpolation of Lagrange and Hermite type, J. Approx. Theory 34 (1982) 286–305. [6] G.T. Herman, Image Reconstruction from Projections: The Fundamentals of Computerized Tomography, Academic Press, 1980, pp. 65–67. [7] F. John, Abhängigkeiten zwischen den Flächenintegralen einer stetigen Funktion, Ber. Verh. Sächs. Akad. 69 (1917) 262–277. [8] X.Z. Liang, On Hakopian interpolation in the disk, Approx. Theory Appl. 2 (1986) 37–45. [9] X.Z. Liang, C.M. Lü, On the convergence of Hakopian interpolation and cubature, J. Approx. Theory 88 (1997) 28–46. [10] X.Z. Liang, C.M. Lü, On the integral convergence of Kergin interpolation on the disk, J. Comput. Appl. Math. 95 (1998) 45–63. [11] B. Logan, L. Shepp, Optimal reconstruction of a function from its projections, Duke Math. J. 42 (1975) 645–659. [12] R. Marr, On the reconstruction of a function on a circular domain from a sampling of its line integrals, J. Math. Anal. Appl. 45 (1974) 357–374. [13] K.I. Oskolkov, Ridge approximation, Chebyshev–Fourier analysis and optimal quadrature formulas, Tr. Mat. Inst. Steklova 219 (1997) 269–285 (Mi tm620). [14] J. Radon, Über die Bestimmung von Funktionen durch ihre Integralwerte längs gewisser Mannigfaltigkeriten, Math. Ann. 111 (1935) 541–559. [15] G.N. Ramachandran, A.V. Lakshminarayanan, Three dimensional reconstruction from radiographs and electron micrographs: application of convolutions instead of Fourier transforms, Proc. Natl. Acad. Sci. 68 (1971) 2236–2240. [16] T.J. Rivlin, Chebyshev Polynomial: From Approximation Theory to Algebra and Number Theory, 2nd edition, Wiley, New York, 1990. [17] L.A. Shepp, B.F. Logan, The Fourier reconstruction of a head section, IEEE Trans. Nucl. Sci. NS-21 (1974) 21–43. [18] D.C. Solmon, The X-ray transform, J. Math. Anal. Appl. 56 (1976) 61–83. [19] Yuan Xu, A new approach to the reconstruction of images from Radon projections, Adv. Appl. Math. 36 (2006) 388–420. [20] Yuan Xu, Olge Tischenko, Fast oped algorithm for reconstruction of images from Radon data, East J. Approx. 12 (2007) 427–444. [21] Yuan Xu, O. Tischenko, C. Heoschen, Fast implementation of the image reconstruction algorithm OPED, in: Medical Imaging 2009: Physics of Medical Imaging, Proc. SPIE 7258 (2009) 72585F.
Contents lists available at ScienceDirect
Applied Numerical Mathematics www.elsevier.com/locate/apnum
Fast and exact 2d image reconstruction based on Hakopian interpolation ✩ Xuenan Sun a,∗ , Xuezhang Liang b a b
School of Mathematics and Statistics, Northeast Normal University, Changchun, 130024, China Institute of Mathematics, Jilin University, Changchun, 130022, China
a r t i c l e
i n f o
Article history: Received 8 September 2016 Received in revised form 20 May 2017 Accepted 13 July 2017 Available online 21 July 2017
a b s t r a c t A new algorithm for the reconstruction of two-dimensional (2D) images from projections is given. The algorithm is based on the expansions of Hakopian interpolation polynomial into Chebyshev–Fourier series. The computer simulation experiments show that the new algorithm is effective. © 2017 IMACS. Published by Elsevier B.V. All rights reserved.
Keywords: Hakopian interpolation Chebyshev–Fourier series Image reconstruction Radon transformation
1. Introduction The fundamental problem of the computed tomography (CT) is to reconstruct an image from its Radon projections. More precisely, we consider the case of a two dimensional image, which is described by a density function f (x, y ) defined on the unit disk D = {(x, y )|x2 + y 2 ≤ 1} of the R2 plane. Let f be any absolutely integrable function on D. Then a Radon projection of f is a line integral as follows,
Rθ ( f ; t ) :=
f (x, y )dxdy , ∀ 0 ≤ θ ≤ 2π
and
−1 ≤ t ≤ 1,
I (θ,t )
where I (θ, t ) = {(x, y ) ∈ R2 | x cos θ + y sin θ = t } ∩ D is a line segment inside D. By the Radon transform f → {Rθ ( f ; t )}, f is associated to a family of univariate functions of t, parameterized by θ . The essential problem of image reconstruction is to recover f from its Radon projections. It is well known that the set of Radon projections {Rθ ( f ; t ) | |t | ≤ 1, 0 ≤ θ ≤ 2π } determines f completely. This was proved by Radon [14] and John [7] for differentiable functions. Furthermore, if f ∈ L 1 (R2 ) has compact support in D, then f is uniquely determined by any infinite set of Radon projections [18]. However, only a finite number of projections can be measured in practice. Thus, the main problem of image reconstruction is to obtain an effective algorithm that produces a good approximation to f from a large collection of its Radon projections. For the practical reconstruction of two-dimensional images from their measured projections, one most frequently applies the method of filtered backprojection (FBP) and its modifications [15]. These algorithms can be viewed as computer imple-
✩ This work was partially supported by the National Science Foundation of Jilin Province (No. 20170101044JC), the Education Department of Jilin Province (No. JJKH20170904KJ) and by the National Natural Science Foundation of China (No. 11171056). Corresponding author. E-mail address: [email protected] (X. Sun).
*
http://dx.doi.org/10.1016/j.apnum.2017.07.002 0168-9274/© 2017 IMACS. Published by Elsevier B.V. All rights reserved.
186
X. Sun, X. Liang / Applied Numerical Mathematics 121 (2017) 185–197
mentation of an inversion of the Radon transform through the central slice theorem [1]. They use a discrete convolution of the projections with a filter function either in the spatial domain or in the frequency domain, plus backprojection. Another class of solutions utilizes a decomposition of the Radon transform into circular harmonics and certain polynomials [4,3]. The Chebyshev polynomials of the second kind play a key role in this context [2,20,12]. In [3], an algorithm was described based on the decomposition of the projections into Chebyshev polynomials of the second kind, which uses the backprojection technique. Consequently, the method can be considered as a Chebyshev-domain filtered backprojection (CDFBP). The method combines the exactness of the polynomial approach with the speed of FBP. In [19], an algorithm called OPED was provided, as it was based on orthogonal polynomial expansion on the unit disk. The numerical tests showed that the algorithm was fast and stable, and it produced high quality images [21]. In comparison with FBP, these algorithms are the kind which have slightly higher mathematical complexity and less intuitive expression. In tomography, the data often consists of values of linear integrals over segments. In many situations, a table of mean values of a function of d variables on (d − 1)-dimensional hyperplane is considered to be the most natural type of data for multivariate functions. Hakopian’s famous interpolation formula [5] is an important reason to take this approach. Hakopian proved that for any given n distinct points X 1 , · · · X n on the boundary of a disk D, the set of integrals of f over all the linear segments [ X i , X j ] determines uniquely every polynomial f of degree n − 2. A Lagrange representation of the bivariate Hakopian interpolation was given in [8]. In fact, the first method of recovering a polynomial of degree n from its integrals over chords defined by equally spaced points { X i } on the boundary of the unit disk was described by Marr in [12]. Interpolation methods based on integrals over chords can be used for approximate reconstruction of functions from their Radon transforms. Because of the importance of such recovery methods for applications in tomography, they have been intensively studied (see [12,11]). In this paper, we develop an efficient and simple reconstruction algorithm based on the Hakopian interpolation method, which uses the Chebyshev–Fourier analysis on the unit disk. Therefore, the new algorithm can be implemented by fast discrete sine transform and an interpolation step. Indeed, it combines the exactness of the polynomial approach with the speed of FBP. The rest of this paper is organized as follows: In Section 2, we discuss the background of Hakopian interpolation and orthogonal polynomials. Section 3 is devoted to the Chebyshev–Fourier orthogonal expansions of Hakopian interpolation are studied. Readers who are not interested in the mathematical details may skip most of this section until its end, where a step-by-step outline of the algorithm is given. Reconstructions are made of the Shepp and Logan phantom in Section 4. The numerical example is given to show that the algorithm is effective. 2. Preliminaries In this section, we present some known results, which are useful later. To begin with, we denote by n (R2 ) the following set of all real algebraic polynomials in two variables of total degree n,
n (R2 ) :=
C j ,k x j yk
x, y , C j ,k ∈ R .
j +k≤n
We shall consider real functions f on the plane R2 . Also, a point x of R2 is denoted by x = (x, y ) . Throughout this paper, the functions f which are to be approximated are assumed to have their support in the unit disk D := {(x, y ) ∈ R2 |x2 + y 2 ≤ 1}. Denote by C (D) the space of all continuous functions defined on D, and define
E m ( f ) :=
inf
max | f (x) − p (x)|, ∀ f ∈ C (D).
p ∈m (R2 ) x∈D
For a given natural number n ≥ 2, we choose the nodes of Hakopian interpolation on the disk which are defined as follows:
xk = (cos
2kπ n
, sin
2kπ n
) ,
k = 1, 2, · · · , n .
Introduce the following integral functional:
1 f (xk + t (xl − xk ))dt , ∀ f ∈ C (D),
f = [xk ,xl ]
0
and denote by λ(k, l) the Radon projections of two univariate function f (x), that is
λ(k, l) = [xk ,xl ]
f,
1 ≤ k < l ≤ n.
X. Sun, X. Liang / Applied Numerical Mathematics 121 (2017) 185–197
187
Moreover, let
U m−1 (t ) = T m (t ),
T m (t ) = cos m arccos t , R m−1 (t ) =
T m (t ) + T m−1 (t ) t+1
S m−1 (t ) =
,
T m (t ) − T m−1 (t )
x = (x, y ) = (r cos θ, r sin θ) , and t = r cos (
t−1 k+l n
,
π − θ).
Next, we recall the following known results in [8,9]. Lemma 2.1. For any f (x) ∈ C (R2 ), there exists a unique polynomial P (x) = H n ( f ; x) ∈ n−2 (R2 ), such that
f − P = 0,
0 ≤ k < l ≤ n.
[xk ,xl ]
Here P (x) = H n ( f ; x) denotes the Hakopian interpolation polynomial of f (x, y ) ∈ C (D) with respect to the n nodes. Furthermore, we have the following representation:
H n ( f ; x) =
L k,l (x) ·
1≤k
f,
(1)
[xk ,xl ]
where (1) If n is even and let n = 2m, then
L k,l (x) =
sin2 l−nk π T m (t )
m2 (t − cos l−nk π )
L k,l (x) = −
2U m−1 (t ) −
sin2 l−nk π U m−1 (t ) m2 (t − cos l−nk π )
2T m (t ) +
T m (t ) t − cos l−nk π
when l − k is odd,
,
(1 − t 2 )U m−1 (t ) m2 (t − cos l−nk π )
(2)
,
when l − k is even.
(3)
(2) If n is odd and let n = 2m − 1, then
L k,l (x) =
l−k 2 2(sin 2m −1 π ) R m−1 (t )
l−k (2m − 1)2 (t − cos 2m −1 π )
2(1 + t ) R m −1 (t )
− L k,l (x) =
l−k 2 2(sin 2m −1 π ) S m−1 (t )
l−k (2m − 1)2 (t − cos 2m −1 π )
l−k (1 + cos 2m −1 π ) R m−1 (t )
,
when l − k is odd,
,
when l − k is even.
l−k t − cos 2m −1 π
2(1 − t ) S m −1 (t )
−
l−k (1 − cos 2m −1 π ) S m−1 (t ) l−k t − cos 2m −1 π
Lemma 2.2. There exists an absolute constant A, such that for arbitrary x ∈ D, f (x) ∈ C (D) and n ≥ 2, we have
| f (x) − H n ( f ; x)| ≤ A E n−2 ( f ) · n log n. Note that, the Radon projection is easily computed in the special case when f has a “plane wave” structure. More precisely, let un (t ) denote the Chebyshev polynomial of the second kind, that is
1 un (t ) = √ D n (arccos t ),
π
where
D n (θ) =
sin(n + 1)θ sin θ
,
n = 0, 1 , · · · .
188
X. Sun, X. Liang / Applied Numerical Mathematics 121 (2017) 185–197
These classical polynomials constitute a complete orthonormal system U = {un−1 (t )}n∞=1 in l2ω (−1, 1) with the weight
√
ω(t ) = 2 1 − t 2 .
After substitution t = x · ξ , for x ∈ D and ξ ∈ S 1 (S 1 is the set of unit vector). Then un (x · ξ ) generate a family of complete orthonormal systems of ridge polynomials in L 2 (D) as follows:
U S 1 = {{un−1 (x · ξ )}n∞=1 }ξ ∈ S 1 . The fundamental properties of this system are expressed by the following [13,16]: (1) Orthogonality relation
un (x · ξ ) p (x)dx = 0,
∀ p (x) ∈ n−1 (R2 ),
n = 1, 2, · · · ,
∀ξ ∈ S 1 ,
D
or un (x · ξ )⊥n−1 (R2 ) in L 2 (D). (2) Inner products of Chebyshev ridge polynomials
un (x · ξ )un (x · η)dx = D
√
2π
un (1)
√ =
π
n+1
un (ξ · η),
∀ξ, η ∈ S 1 .
(4)
τn± denote the subspace of trigonometric polynomials a(·) on S 1 , of degree n and satisfying a(−ξ ) ≡
Furthermore, let (−1)n a(ξ ). Then
1
un (ξ · η)
π un (ξ · η)a(η)dη = a(ξ ), a(ξ ) ∈ τn± , ξ ∈ S 1 ,
S1
and in particular
1
2
√
un (ξ1 · η)un (η · ξ2 )dη = un (ξ1 · ξ2 ),
π
ξ 1 , ξ2 ∈ S 1 .
S1
(3) Integral representation Given a function f (x) ∈ L 2 (D), consider the following Chebyshev–Fourier coefficients, depending on ξ ∈ S 1 as a parameter:
an ( f , ξ ) =
f (y)un (y · ξ )dy,
n = 0, 1 , · · · , ξ ∈ S 1 .
D
The following relation represents the integral form of Chebyshev ridge polynomial Fourier series, which is in fact Radon inversion formula for f (x) expressed via Chebyshev ridge polynomials: ∞ 1 f (x) = n an−1 ( f , ξ )un−1 (x · ξ )dξ. 2π L 2 (D)
n =1
S1
3. Chebyshev–Fourier expansions of Hakopian interpolation Our reconstruction algorithm will be based on this Chebyshev–Fourier orthogonal expansion. Radon–Fourier analysis via Chebyshev ridge polynomials is crucial in the derivation for the algorithm. In order to get the expansions of Hakopian interpolation, we need the following lemma [10]. Lemma 3.1. For every natural number n ≥ 2 and 1 ≤ j < k ≤ n, we have
π 0
π 0
cos2 n2 θ cos θ − cos
k− j n
π
sin2 n2 θ cos θ − cos
k− j n
π
d θ = 0,
if k − j is odd,
(5)
d θ = 0,
if k − j is even.
(6)
X. Sun, X. Liang / Applied Numerical Mathematics 121 (2017) 185–197
189
By making use of the orthogonality of Jacobi polynomials with respect to the weight function (1 − x)±1/2 (1 + x)±1/2 , the proof of Lemma 3.1 is not difficult. Here we omit the proof of it [10]. For simplicity, we set
k+l
η = cos
n
t = r cos
π , sin
l+k n
k+l n
= (cos β, sin β) ∈ S 1 ,
π
c = ωk ,l (cos
π −θ
l−k n
π ).
= r cos θ cos β + r sin θ sin β
= x cos β + y sin β = x · η, and let
ϕ (t ) :=
ωk ,l (t ) ωk ,l (cos l−nk π )
ωk ,l (t )
=
c
= Lk,l (x).
Now, from the orthogonality of the Chebyshev polynomials un (t ), it follows that
ϕ (t ) =
n −2
c j u j (t ),
(7)
j =0
with the coefficients given by
cj =
1
1 − t 2 ϕ (t )u j (t )dt
−1
=
1 c
1
1 − t 2 ωk ,l (t )u j (t )dt
−1
π
1
= √ c
π
sin θ ωk ,l (cos θ) sin( j + 1)θ dθ
0
j+1
= √ c
π
ωk,l (cos θ) cos( j + 1)θ dθ.
π 0
In the following, we compute the coefficient for the Hakopian interpolation expansion in two cases, respectively. Case 1: If n is even, let n = 2m. It is easy to show that
ωk,l (t ) =
⎧ 2 t) ⎨c 0 cos (m arccos l−k t −cos n π 2 t) ⎩c 0 sin (m arccos t −cos l−nk π
c = ωk ,l (cos
l−k n
π ) = c0
when l − k is odd, when l − k is even. m2
sin2 l−nk
π
,
where c 0 is a constant. Set α = l−nk π . From Lemma 3.1, we obtain that
( j + 1)c 0 cj = √ c
π
π
0
=
( j + 1)c 0 √ 2c
2π
π
0
=
( j + 1)c 0 √ 4c
2π
π
0
cos2 mθ cos θ − cos α cos2 mθ cos θ − cos α
cos( j + 1)θ dθ
(cos( j + 1)θ − cos( j + 1)α )dθ
sin( j + 1) θ −2 α sin θ −2 α
·
sin( j + 1) θ +2 α sin θ +2 α
dθ
190
X. Sun, X. Liang / Applied Numerical Mathematics 121 (2017) 185–197
( j + 1)c 0 = √
π
sin( j + 1)(θ − α ) sin(θ − α )
π
2c
0
=
√ ( j + 1)c 0 π
sin( j + 1)α
·
c
sin α
·
sin( j + 1)θ sin θ
dθ
,
and consequently,
ϕ (t ) =
n −2
√ ( j + 1)c 0 π
·
c
j =0
sin( j + 1)α sin α
· u j (t ) .
If let t s = cos ns π , then
ϕ (t ) =
n −2 ( j + 1)π c 0
c
j =0
u j (tl−k )u j (t ) .
Therefore,
L k,l (x)u i (x · ξ )dx =
n −2 ( j + 1)π c 0
c
j =0
D
= = Then, since by
c0 c
=
4 sin2 α , n2
(i + 1)π c 0 c
c0 π
√
c
π
u j (x · η)u i (x · ξ )dx
u j (tl−k )
√ u i (tl−k )
D
π
i+1
u i (ξ · η)
u i (tl−k )u i (ξ · η).
we can further have
L k,l (x)u i (x · ξ )dx = 4π
sin l−nk π sin(i + 1) l−nk π n2
D
u i (ξ · η).
With the aid of (6), similarly, we can obtain the same result when l − k is even. Case 2: If n is odd, let n = 2m − 1, then
ωk,l (t ) =
⎧ 2 ⎨c 1 (cos(m arccos t )+cos((ml−−k 1) arccos t )) (t +1)(t −cos n π ) ⎩c 1 (cos(m arccos t )−cos((ml−−k 1) arccos t ))2 (t −1)(t −cos n π )
when l − k is odd, when l − k is even,
and
c = ωk ,l (cos
l−k n
π ) = c1
(2m − 1)2 2 sin2 l−nk π
,
where c 1 is a constant. Analogously to the Case 1, we will only derive the case that l − k is odd. From Lemma 3.1 we obtain that
( j + 1)c 1 cj = √ c
π
π
0
=
( j + 1)c 1 √ c
2π
π
0
=
( j + 1)c 1 √ c
π
π
0
=
2( j + 1)π c 1 c
(cos mθ + cos((m − 1)θ))2 cos( j + 1)θ dθ (cos θ + 1)(cos θ − cos α ) (cos 2m2−1 )2 cos( j + 1)θ dθ cos θ − cos α sin( j + 1)(θ − α ) sin(θ − α ) u j (tl−k ),
·
sin( j + 1)θ sin θ
dθ
X. Sun, X. Liang / Applied Numerical Mathematics 121 (2017) 185–197
and therefore, we can get the orthogonal expansions of
ϕ (t ) =
n −2 2( j + 1)π c 1
c
j =0
191
ϕ (t ):
u j (tl−k )u j (t ).
(8)
It follows that
L k,l (x)u i (x · ξ )dx =
n −2
2( j + 1)π c 1 c
j =0
D
3 2
2(π ) c 1
= and because of
c1 c
=
2 sin2 α , n2
c
u j (x · η)u i (x · ξ )dx
u j (tl−k ) D
u i (tl−k )u i (ξ · η),
we have
L k,l (x)u i (x · ξ )dx = 4π
sin l−nk π sin(i + 1) l−nk π n2
D
u i (ξ · η).
Combining Case 1 and Case 2, for arbitrary natural number n, we obtain that
L k,l (x)u i (x · ξ )dx = 4π
sin l−nk π sin(i + 1) l−nk π n2
D
u i (ξ · η).
(9)
Therefore, with the aim of (9), we have the Chebyshev–Fourier coefficient of H n ( f ; x):
⎛
⎜ L k,l (x) ⎝
ai ( H n ( f ; x), ξ ) =
1≤k
D
⎛
=
⎝
1≤k
=
⎞
⎟
f ⎠ u i (x · ξ )dx
[xk ,xl ]
⎞
L k,l (x)u i (x · ξ )dx⎠ ·
D
4π
f
[xk ,xl ]
sin l−nk π sin(i + 1) l−nk π n2
1≤k
u i (ξ · η)
f
[xk ,xl ]
Especially when i ≥ n − 1, ai = 0. Then,
H n ( f ; x) =
=
=
n −2 1
2π 2 n2 4
ω =1 n −2
π
n2
S1
u ω−1 (ξ · η)u ω−1 (x · ξ )dξ
ω
ω =1
√
aω−1 ( H n ( f ; x), ξ )u ω−1 (x · ξ )dξ
ω
1≤k
S1
n −2
ω =1
ω
sin
l−k
1≤k
n
π sin ω
l−k n
sin
l−k n
π sin ω
l−k n
π
f [xk ,xl ]
π u ω−1 (t )
f,
[xk ,xl ]
where t = x cos( l+nk π ) + y sin( l+nk π ). Without loss of generality, we suppose that n is even and let n = 2m (the similar result can be obtained if n is odd). Then,
H n ( f ; x) =
n −2 4
n2
ω =1
ω
1
+
2
(10)
192
X. Sun, X. Liang / Applied Numerical Mathematics 121 (2017) 185–197
where
1
=
m m
t = x cos
2
2i − 1 2m
j =1 i =1
sin
=
2j −1 2m
m m −1
sin
j =1 i =1
t = x cos
j m
i
π sin ω
2i − 1 2m
π + y sin
π u ω−1 (t )
2j −1 2m
j = 1, 2, · · · , m ,
π sin ω π u ω−1 (t ) m
π + y sin
(11)
[x j−i ,xi + j−1 ]
π
i
m
f
j m
f
(12)
[x j−i ,x j+i ]
j = 1, 2, · · · , m .
π
k In addition, the index k of x in (11) and (12) will probably be negative or exceed n. In these two cases, we define k±n x =x . That is, in 1 , if k = j − i < 0, then we set k = j − i + n; if l = i + j − 1 > n, then we set l = i + j − 1 − n. In 2, if k = j − i < 0, then we set k = j − i + n; if l = i + j > n, then we set l = i + j − n. For simplicity, we introduce the following notations k
(1 )
c ω, j =
m
2i − 1
sin
n
i =1
(2 )
c ω, j =
m −1 i =1
sin
2i n
π sin ω
π sin ω
2i n
2i − 1 n
π
f, [x j−i ,xi + j−1 ]
π
f. [x j−i ,xi + j ]
Then, these findings can be translated into a simple reconstruction algorithm that consists of the following steps: (1)
(2)
Step 1. Computing c ω, j and c ω, j , respectively, by projections. Step 2. Computing the summation of λ1 (t , j ) = Step 3. Computing the summation H n ( f ; x) =
n −2
ωcω(1,)j u ω−1 (t ) and λ2 (t , j ) =
ω=1 √ m 4 π n2
j =1
n −2
ω=1
ωcω(2,)j u ω−1 (t ).
{λ1 (t , j ) + λ2 (t , j )}.
In reconstruction of step 1, the c ω, j and c ω, j can be directly calculated by a discrete sine transform of sin 2in−1 π × 2i [x j −i ,xi + j −1 ] f and sin n π [x j −i ,xi + j ] f . Very efficient algorithms requiring only on the order of (m + 1) · log(m + 1) operations can be applied for this purpose [3]. Since the sine transform is its own inverse, we can also use it to speed up the computation of λ1 (tk , j ) and λ2 (tk , j ) (step 2), where tk = cos knπ , k = 1, 2, · · · , n − 1. Values at other points t, which are needed for the reconstruction, are obtained by linear interpolation. (1)
(2)
4. Examples of the algorithm To demonstrate the practical usefulness of the algorithm, we apply it to the reconstruction of the Shepp–Logan head phantom [17] using line integrals. This phantom head is well known because of its use in testing the accuracy of reconstruction algorithms with regards to their ability to reconstruct cross sections of the human head with X-ray tomography. A particular advantage of this phantom is that it is constructed from the superposition of 10 ellipses of homogeneous “density” and thus, based on the linearity of the Radon Transform, the projections of this construction can be specified analytically. This phantom has highly singularities as the image contains jumps at the boundary of every ellipse in the image, include the one on the boundary. The function that represents the image is a step function, which is not continuous at the boundary of each of the ellipses. We do the reconstruction on a 256 × 256 grid using 256 projection angles. The original image and the reconstruction are given in Figs. 1 and 2. In addition to the Shepp–Logan phantom, we reconstruct another image on a 280 × 280 grid using 256 projection angles. The original function named “intersected ellipse” is defined by ourselves. The original image and the reconstruction are shown in Figs. 3 and 4. Furthermore, a more quantitative way of evaluating pictures is the following. Select a column of pixels, which is such that it goes through a number of interesting features in the original. For example, in Shepp–Logan head phantom the 80th
X. Sun, X. Liang / Applied Numerical Mathematics 121 (2017) 185–197
193
Fig. 1. Shepp–Logan head phantom.
Fig. 2. Reconstruction image by this algorithm.
Fig. 3. The original intersected ellipse.
of the 256 columns goes through the above ellipse, the left rotated ellipse, and the left horizontal ellipse. One way to evaluate the quality of reconstruction is compare the graphs of the 256 pixel densities for this column in the original and the reconstruction. Fig. 5 shows the values of the original and the reconstruction along the 80th column. At the same time, we introduce two picture distance measures to measure the closeness of the reconstruction to the original [6]. We use t u , v and r u , v to denote the densities of the vth pixel of the uth row of the test phantom and reconstruction respectively, and t to denote the average of the densities in the test phantom. We assume that both pictures are l × l. The measures are defined in the following:
194
X. Sun, X. Liang / Applied Numerical Mathematics 121 (2017) 185–197
Fig. 4. Reconstruction image by this algorithm.
Fig. 5. Line plots of the 80th column of phantom and reconstruction.
⎡
l l
(t u , v − ru , v ) ⎢ ⎢ u =1 v =1 l l ⎣ (t u , v − t )2
2
d=⎢
⎤ 12 ⎥ ⎥ ⎥ , ⎦
u =1 v =1
l
r=
l
u =1 v =1
|t u , v − ru , v |
l l u =1 v =1
. |t u , v |
These two measures emphasize different aspects of picture quality. The first one, d, is a normalized root mean square distance measure. A large difference in a few places causes the value of d to be large. The second one, r, is a normalized mean absolute distance measure. As opposed to d, it emphasizes the importance of a lot of small errors rather than of a few large errors [6]. Picture distance measures for reconstructions from different projection data are shown in Tables 1, 2, 3 and 4. Table 1 and Table 2 show the computational results of the Shepp–Logan phantom. The computational results of the “intersected ellipse” are given in Table 3 and Table 4. From these tables we can see that the reconstruction obtained from Hakopian interpolation algorithm compares well with reconstructions from FBP method, and the reconstructions from the Hakopian algorithm and CD-FBP method are almost the same. Moreover, the calculation time for the reconstructions (256 × 256 Shepp–Logan phantom) is given in Table 5. We can use it as a reference to the computational complexity of the three methods. Table 5 shows that the computational time of FBP method is shorter than the other ones. Meanwhile, the computational time of the Hakopian algorithm can be acceptable.
X. Sun, X. Liang / Applied Numerical Mathematics 121 (2017) 185–197
195
Table 1 Picture distance measures for Shepp–Logan phantom. Number of projection angles
120 128 256
CD-FBP algorithm
Hakopian algorithm
d
r
d
r
0.197571 0.191253 0.144232
0.100863 0.095916 0.051576
0.203349 0.197040 0.144907
0.106673 0.104356 0.068250
Table 2 Picture distance measures for Shepp–Logan phantom. Number of projection angles
120 128 256
FBP algorithm
Hakopian algorithm
d
r
d
r
0.406432 0.395613 0.360563
0.370060 0.359921 0.352121
0.203349 0.197040 0.144907
0.106673 0.104356 0.068250
Table 3 Picture distance measures for “intersected ellipse”. Number of projection angles
64 128 256
CD-FBP algorithm
Hakopian algorithm
d
r
d
r
0.130408 0.093564 0.071143
0.042846 0.027616 0.018604
0.135765 0.097128 0.073810
0.047598 0.030880 0.021450
d
r
d
r
0.144536 0.100776 0.097522
0.048256 0.030940 0.023176
0.135765 0.097128 0.073810
0.047598 0.030880 0.021450
Table 4 Picture distance measures for “intersected ellipse”. Number of projection angles
64 128 256
FBP algorithm
Hakopian algorithm
Table 5 The computational time of three methods. Algorithm
128 projection angles
256 projection angles
512 projection angles
FBP CD-FBP Hakopian
about 2 s about 9 s about 11 s
about 2 s about 12 s about 39 s
about 3 s about 13 s about 93 s
Sometimes, people are only interested in a local region of the reconstructed image. So we reconstruct the left and right rotated ellipse of the Shepp–Logan from the local projection data. The projections are line integrals within the region of interest. The others are zero outside of the local region. The left local region is 45 pixels in radius in a 200 × 200 image. The right local region is 35 pixels in radius. The local reconstructions are shown in Figs. 6 and 7. The figures show that the algorithm based on Hakopian interpolation is able to reconstruct the region of interest. All reconstructions were done on a 1.80 GHz Intel(R) Core(TM) i5-3337U CPU personal computer. The programming language is Visual C++6.0, and the operating system is Windows 7 Professional. We used the FFT for discrete sine transform in the package FFTW (http://www.fftw.ort/). 5. Conclusions A fast and easy to implement algorithm for 2D image reconstruction has been developed. Mathematically, basis is a decomposition of Hakopian interpolation into Chebyshev polynomials. To the best of our knowledge the described implementation in form of this algorithm is new, although there is a loose relationship to other recent developments ([3,2,20]). The main advantage of the algorithm as compared with the FBP method is its exactness. We introduce a fast implementation of Hakopian interpolation algorithm by using FFT and a linear interpolation step. We hope that our algorithm will be useful for both practical and theoretical problems.
196
X. Sun, X. Liang / Applied Numerical Mathematics 121 (2017) 185–197
Fig. 6. The local reconstruction of left rotated ellipse.
Fig. 7. The local reconstruction of right rotated ellipse.
Acknowledgements The authors thank the anonymous reviewers for their thorough reading of the paper and their constructive comments. Appendix A In this part, we need to explain the relation between line integrals and projection angles. As well as we know, every straight line in the plane will be characterized as the locus of points satisfying
x cos θ + y sin θ = t ,
(A.1)
for some real distance parameter t and angle θ . For convenience, we allow t to assume negative values in (A.1), so that every straight line has exactly two distinct parametric representations, (t , θ) and (t , θ ), with t = −t, θ = θ + π (mod 2π ). We mean the set of parameters by C,
C = {(t , θ)| − 1 ≤ t ≤, 0 ≤ θ ≤ 2π }. Then every chord, or straight line intersecting D, is represented by exactly two points on C. Now, we return to the algorithm described in this paper. We consider an interpolation nodes set, X n = {xk |xk = (cos 2(k−n1)π , sin 2(k−n1)π ) , k = 1, 2, · · · , n}, of n points spaced uniformly around the circumference of D. We can suppose that these points are numbered from 1 to n, in counterclockwise order, and for simplicity assume that point 1 lies along the positive x-axis. Then, xk are the coordinates of the k-th point. Therefore, the n(n − 1)/2 possible chords joining pairs of points in X n correspond to a mesh of n(n − 1) points on C by (A.1). To describe these, we define
tl = cos(π l/n), l = 1, 2, · · · , n − 1,
(A.2)
and
θk,l =
π n
(2k + l − 2), k = 1, 2, · · · , n, l = 1, 2, · · · , n − 1.
(A.3)
X. Sun, X. Liang / Applied Numerical Mathematics 121 (2017) 185–197
197
Then, the point (tl , θk,l ) on C corresponds to the chord joining the two points in X n numbered k1 and k2 , where
k1 = k < k2 = k + l,
if k + l = n,
k1 = k + l − n < k2 = k,
if k + l > n.
Thus, we connect the number of line integrals with the number of projection angels. This situation is the same as the discrete case in [12]. References [1] H.H. Barrett, The Radon transform and its applications, in: E. Wolf (Ed.), Progress in Optics XXI, Elsevier, Amsterdam, 1984. [2] B. Bojanov, I.K. Georgieva, Interpolation by bivariate polynomials based on Radon projections, Stud. Math. 162 (2004) 141–160. [3] Thomas Bortfeld, Uwe Oelfke, Fast and exact 2D image reconstruction by means of Chebyshev decomposition and backprojection, Phys. Med. Biol. 44 (1999) 1105–1120. [4] A.M. Cormack, Representation of a function by its line integrals with some radiological applications, J. Appl. Phys. 34 (1963) 2722–2727. [5] H.A. Hakopian, Multivariate divided difference and multivariate interpolation of Lagrange and Hermite type, J. Approx. Theory 34 (1982) 286–305. [6] G.T. Herman, Image Reconstruction from Projections: The Fundamentals of Computerized Tomography, Academic Press, 1980, pp. 65–67. [7] F. John, Abhängigkeiten zwischen den Flächenintegralen einer stetigen Funktion, Ber. Verh. Sächs. Akad. 69 (1917) 262–277. [8] X.Z. Liang, On Hakopian interpolation in the disk, Approx. Theory Appl. 2 (1986) 37–45. [9] X.Z. Liang, C.M. Lü, On the convergence of Hakopian interpolation and cubature, J. Approx. Theory 88 (1997) 28–46. [10] X.Z. Liang, C.M. Lü, On the integral convergence of Kergin interpolation on the disk, J. Comput. Appl. Math. 95 (1998) 45–63. [11] B. Logan, L. Shepp, Optimal reconstruction of a function from its projections, Duke Math. J. 42 (1975) 645–659. [12] R. Marr, On the reconstruction of a function on a circular domain from a sampling of its line integrals, J. Math. Anal. Appl. 45 (1974) 357–374. [13] K.I. Oskolkov, Ridge approximation, Chebyshev–Fourier analysis and optimal quadrature formulas, Tr. Mat. Inst. Steklova 219 (1997) 269–285 (Mi tm620). [14] J. Radon, Über die Bestimmung von Funktionen durch ihre Integralwerte längs gewisser Mannigfaltigkeriten, Math. Ann. 111 (1935) 541–559. [15] G.N. Ramachandran, A.V. Lakshminarayanan, Three dimensional reconstruction from radiographs and electron micrographs: application of convolutions instead of Fourier transforms, Proc. Natl. Acad. Sci. 68 (1971) 2236–2240. [16] T.J. Rivlin, Chebyshev Polynomial: From Approximation Theory to Algebra and Number Theory, 2nd edition, Wiley, New York, 1990. [17] L.A. Shepp, B.F. Logan, The Fourier reconstruction of a head section, IEEE Trans. Nucl. Sci. NS-21 (1974) 21–43. [18] D.C. Solmon, The X-ray transform, J. Math. Anal. Appl. 56 (1976) 61–83. [19] Yuan Xu, A new approach to the reconstruction of images from Radon projections, Adv. Appl. Math. 36 (2006) 388–420. [20] Yuan Xu, Olge Tischenko, Fast oped algorithm for reconstruction of images from Radon data, East J. Approx. 12 (2007) 427–444. [21] Yuan Xu, O. Tischenko, C. Heoschen, Fast implementation of the image reconstruction algorithm OPED, in: Medical Imaging 2009: Physics of Medical Imaging, Proc. SPIE 7258 (2009) 72585F.