A Proposal for Combining FBP and ART in CT-reconstruction

2003 Meeting on Fully 3D-reconstruction, St Malo, June 30 –July 4, 2003 1(4)

A Proposal for Combining FBP and ART in CT-reconstruction

Per-Erik Danielsson and Maria Magnusson Seger

Computer Vision Laboratory, Dept. of EE, Linkoping University, SE-581 83 Linkoping, Sweden ped@isy.liu.se; maria@isy.liu.se

Abstract. This paper presents a proposal (not entirely new) for combining analytical and algebraic reconstruction techniques. Such a combination bears the promise to improve the image quality of fast but non-exact, reconstruction of the filtered backprojection type. The difference between the present proposal and traditional ART is that we compute a full error image with FBP, applied to projection differences, to update the solution in each iteration step. The main road-block seems to be the same that has been an obstacle for many ART-algorithms in CT applications, namely that the forward projections are subjected to ailiasing, which tend to override the intended benefits of the updating loop. We present an analysis of this problem and indicate some possible solutions.

The iteration loop

Iterative methods such as Algebraic Reconstruction Techniques (ART), in combination with analytical Fourier methods, such as Filtered Backprojection (FBP), have been considered in the past by several workers in the field. It comes rather naturally to suggest ART as a post-processing step that could possibly improve a reconstruction result, primarily obtained by FBP, or equivalently, to propose an initial FBP reconstruction step to give the ART process a head start. Surprisingly few articles on this matter can be found in the CT- literature, however. In this paper we like to propose the reconstruction scheme shown in

Fig. 1. This is not very different from the ILIN180-method, investigated by Nuyts et al in [1]. Similarities and differences will be discussed shortly.

Definitions and notations for our proposal are inserted at various places in Fig.1, which should be interpreted as follows. f(x,y,z) is the original volume and

) , , ( t s

pθ are the measured projections. Assume that the latter are utilized in a non-exact reconstruction attempt, e.g. for space-invariant one-dimensional filtering in 3D-reconstruction from helical cone-beam projections. In any case, such filtering followed by backprojection produces the volume . The iterative part of the algorithm then produces forward projected data 1 f ) , s t , (

pi θ , using the same geometry as was employed for the measured datap(θ,t,s). The projection differences ∆pi(θ,t,s)are utilized for

FBP-reconstruction, and the resulting volume(∆f)1 is accumulated to yielding . Note that the iterative loop is not the conventional Kaczmarz method [4] but a full FBP-reconstruction. Also, no updating of takes place until all projections have been generated from its present result, a principle, which is sometimes called simultaneous ART.

1

f f2

i

f

Fig. 1. Proposal for iterative filtered back-projection. Ideally, in an exact reconstruction situation the filtered

back-projection delivers the wanted object function up to a precision that is determined by the detector resolution and number of projections. In operator notation we write the exact reconstruction asBHP f(x,y)= f(x,y). In case the reconstruction is

non-exact we have B so that we do not obtain exactly but, say,

-1 P ≠ H f f f δ + = B f f1 HP ≡ , (1) Filtering

H

) , , ( t s qθ ) , , (x y z f in ) , , ( t s pθ

H

B

i

f

) , , ( t s pi θ i f ) (∆ ) , , ( t s p_i θ ∆

α

(

t s

)

qi θ, , ∆ + + Backproj

.

Proj.

P

_Proj.

_P

Filtering +

-2003 Meeting on Fully 3D-reconstruction, St Malo, June 30 –July 4, -2003 2(4) where δ is the artifact portion of the FBP-f

reconstructed image. In the following iterations we get

(2)

[

]

,... , , i , f δ ( f ) f f ( f ) f ( f f f δ f ) f f ( f ) f ( f f f δ f f δ f δ f f f f f ) f f ( f ) f ( f f i i ) i i i i i 2 1 0 2 For 1 1 1 3 2 2 2 2 3 2 2 1 1 1 1 1 1 1 1 2 = + = − + ≡ ∆ + = + = − + ≡ ∆ + = − = − − − = − + = ≡ = − + ≡ ∆ + = + − + BH P P P P H B P H B P H B P P P P H B in in

Here, should be spelled artifacts of the artifacts, i.e. the artifacts produced by the non-exact reconstruction given the artifact image

f 2 δ

f

δ as input. The artifact component changes sign for each iteration, and in practice it might be wise to employ only a fraction

α

of the full values in the volume(∆f )_i. Let us interpret (1) as the linear system (3). Here, is an -dimensional vector and the linear operator

f N

P H

B is the sum of one identity and one error matrix of size NxN as in (4),

f

f

f

= I

+

δ

P

H

B

, (3) δ + = I P H B . (4)

Intuitively, it seems that the iteration (2) should converge under very liberal and often fulfilled conditions. The convergence rate should be exponential in nature, but will certainly depend strongly on the degree of correctness in the filterH , which is just to say that the convergence must be proportional to the exactness of the inversion processBHP . The iteration (2) belongs to the class of generalized Landweber iterations, the convergence rate of which is discussed at considerable length by Pan et al. in [2]. Another treatment of this convergence problem is found in Press et al. [3]. However, at the present stage we do not have experiments to compare with either of these theoretical results.

Nuyts et al. [1] investigated an algorithm called LIN180, using the updating loop in Fig.1 with the purpose to improve the quality of single detector-row helical CT-reconstruction algorithms. Therefore, the first reconstruction step in LIN180 employs rebinning of measured projection data, from the original projection ray geometry, where the source path and its fan beam rays are “skewing” themselves through a cross-sectional slice, to a parallel set of projection rays all of which are parallel and in-slice over 180 degrees. The negative effect of this rebinning step seems to overshadow the potential improvements. While the

investigation [1] attempted to use ART updating to compensate for a mismatch in projection geometry, we are instead aiming to compensate a mathematically non-exact analytical reconstruction, typically embodied in an oversimplified filtering step. To underline our assumptions, the recursion expression in the bottom line of (2) is clearly dependent on the identity P_in≡P in the top line. Hence, the possibility to identify [2] with the Landweber iteration would be severely tarnished byP_in≠P .

Aliasing in forward projections of sampled

images

As can be understood from the previous discussion, we will always make sure that the basic ray geometry of measured projection data and the ones used for subsequent projections and backprojections will be identical. Even so, can we be certain, that in Fig.1 and the iteration (2), the projections p(θ,t,s)from the object volume have been obtained with an operator that is mathematically identical to the following ones generated from the sampled reconstructions? Fig.2, top row, middle, shows the image , which is a sub-sampled and smoothed version of the mathematically defined Shepp-Logan phantom. Second row shows the result of a standard filtered back-projection reconstruction applied to projection data generated from

(

x fs ,y

)

( )

x y

fs , . The image quality for this image is

appallingly low, much lower than what is obtained with the same reconstruction applied to projections from the original non-sampled phantom. We may rightfully suspect that the technique to compute forward projections from the available samples was inadequate. The culprit is a form of ailiasing that we will discuss in qualitative terms with the help of Fig. 3 and 4. Observations of this kind are not overly common. However, Matej and Lewitt motivate their choice of blobs with similar reasoning in [5].

The projections to be used in the second row of Fig.2 are computed from the sampled image fs

( )

x,y in a common pixel-driven fashion. Each pixel value is transported from the pixel center along a ray to the detector axis, and distributed to the nearest two detector positions with the window (interpolation) function w(t)=Λ(t/∆). This is

i) to assume a continuous function

, wherea is

rotationally symmetric and the inverse Radon transform of the triangle function and then

) , ( ) , ( ) , ( ~ y x f y x a y x f = ∗ _s ( yx, ) ) / ( ∆ Λ t ~

ii) to compute the Radon transform of f( yx, )

2003 Meeting on Fully 3D-reconstruction, St Malo, June 30 –July 4, 2003 3(4)

L Sinogram from the continuous phantom ∆ = ∆ ₃2 t M Sampled and smoothed phantom R Sampled, smoothed, and zoomed phantom

( )

∆ ∆ Λ = / ) (t t w , ∆t =∆ L Sinogram differences between sampled and continuous phantoms M Reconstruction result R Zoomed result

(

)

θ θ cos cos / ) ( ∆ ∆ Λ = t t w ∆ = ∆ ₃2 t L Sinogram differences M Reconstruction result R Zoomed result ) ( ) (t w t w = _sc eq. (5) ∆ = ∆ ₃2 t L Sinogram differences M Reconstruction result R Zoomed result

Figure 2. Reconstructions of the sampled Shepp-Logan phantom.

1 − ∆ u v s f 1 − ∆ ) (τ ) , ( vu F ) ( sinc ) (τ =∆ 2τ∆ H 1 8∆− 1 5∆− 1 2∆− 1 − ∆

τ

Nyquist boundary for _{Copies of}

Zero-crossings for W

Figure 3. The Fourier domain F_s( vu, ) of the sampled image

) , ( yx

fs overlayed with three radial lines on which we find the

contributions from three different projections.

Figure 4. The sinc2_{-function and the repeated}

2003 Meeting on Fully 3D-reconstruction, St Malo, June 30 –July 4, 2003 4(4)

Due to the Fourier slice theorem, we should be able to retrieve the Fourier transforms P

( )

τ,θ

Fs

of these projections in the 2D-Fourier transform of the image in Fig.3. But since the projections have been convolved with the same triangle function for all angles in the above step ii), the corresponding Fourier components have been modulated with the same sinc-function

for all angles in Fig.4. The period of this function (see the two circles in Fig.3) matches the double Nyquist limits of the Fourier transform in the x- and y-directions, which means that aliasing from the DC-component and the lower frequencies in these projections are much subdued. But this is not so for other directions. The strongest aliasing appears in the 45 degree direction, where the first repeated instance of

, including its large DC-component, is intercepted by ) , ( vu

(

x y f_s , ) ( 2 _τ∆ ) , ( vu

)

sinc ) (τ =∆ W F_s ) , ( vu F

( )

t,π₄

P close to the maximum of the second lobe of the sinc2_{-function. The effect is}

directly visible in Fig. 2, second row, left, as a checkered band in the zoomed difference sinogram. We have investigated two alternative methods to combat the above aliasing phenomena. The first one is to interpolate with a triangle function of width ∆cosθ which is to employ the method due to Joseph [7]. The second is to use a more sophisticated interpolation function with very small contributions outside the main lobe. The generalized Kaiser-Bessel function was proposed by Lewitt [6] for latter purpose. Here, we will use what we call the four-point sincot-function [8], which suppresses the side lobes outside

sc

w 1 −

∆ very good. With α=0.605 in the following formula the largest side-lobe in the power spectrum peaks at less than -65 dB.

(

)

(5) 0 Else , 2 , 4 cot sin 2 cos ) 1 ( 4 1 ) ( = ∆ < ∆ ∆ ∆ − + = sc sc w t t t t t w α α π π π

To suppress aliasing from the part of the main lobe that lies outside the Nyquist falls limit, we also have to decrease the sampling density along the detector axis to ∆ = ₃2∆

t

/ (

. The results of these two attempts are shown in the two bottom rows of Fig.2. Note that in both cases, the back-projection in the reconstruction is still done with the linear interpolation kernel Λ t ∆). Although the linear interpolation à la Joseph in the third row produces much better image quality than the result in the second row, it is still inferior to the result in the fourth row, where the computation cost is

approximately twice as high. We expect that forward projections of this type will be necessary to include in the combination of FBP and ART proposed in the first part of this abstract.

Acknowledgement

Some valuable comments on the subject of this paper was given by Gunter Lauritsch. Financial support by the Swedish Foundation for Strategic Research under the COMEX project and by Siemens Medical Solutions, Forchheim, Germany is gratefully acknowledged.

References

[1] J.Nuyts, B. De Man, P.Dupont, M.Defrise, P.Suetens, and L.Mortelmans. Iterative reconstruction for helical CT: a simulation study. Phys. Med. Biol. 43 (1998) 729-737

[2] T-S. Pan, A.E Yagle, N.H. Clinthorne, and W.L. Rogers, Acceleration and filtering in the generalized Landweber iteration using a variable shaping matrix. IEEE Trans. Med. Im. 12 (1993) 278-286

[3] W. H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical Recipies, Second Edition, pp 55-58, Cambridge University Press, 1992

[4] S. Kaczmarz, Angenaherte Auflösung von Systemen linear Gleichungen, Bull. Acad. Pol. Sci. Lett. A., 6-8A (1937) 355-357

[5] S. Matej and R. Lewitt. Efficient 3D grids for image reconstruction using spherically-symmetric volume elements. IEEE Trans. Nuuc. Sci., 42 (1995) 1361-137.

[6] R. Lewitt. Multidimensional digital image representations using generalized Kaiser-Bessel window functions. J.Opt. Soc. Am. A., 7 (1990) 1834 -1846

[7] P.M. Joseph, An improved algorithm for reprojecting rays through pixel images. IEEE Trans. Med. Im. 1 (1982) 192-196

[8] M. Magnusson Seger. Three-dimensional reconstruction from cone-beam data using an efficient Fourier technique combined with a special interpolation filter. Phys. Med. Biol. 43 (1998) 951-959