Image Reconstruction for Fan Beam X-Ray Tomography

By Laurie Patterson,2015-03-15 13:31
14 views 0
Image Reconstruction for Fan Beam X-Ray Tomography

    Image Reconstruction for Fan Beam X-Ray Tomography

    Gerneral Using Hankel Transform Pair (1)

    Shuang Ren Zhao

    Doubletask, Toronto, Canada


    A new integral transform pair has been introduced for the reconstruction of slices from fan beam projection data. This integral transform plays the same role as the Hankel transform (or Fourier Bessel transform) in the parallel beam case and includes the Hankel transform as a special case. The relations between the new integral transform and the Hankel transform pair are derived and the implementation of the algorithm is described. Memory space requirements and speed of the transform are discussed and the advantage of precalculations are mentioned. The transform pair is especially useful in case the two dimensional object slice to be reconstructed has only few low frequency components in the azimutal direction. Furthermore the method is capable of treating the case where non analytic a priori information is to be taken into account. In this case projections into the measurement directions are performed and the differences between these projection data and the measured data are reconstructed in order to improve the object data. Alternating transformations between projection and object space are efficient and fast. Finally some results of simulations including the influence of noise are shown and possible applications of the method are discussed.

    Keywords;Image reconstruction, fan-beam, tomography, Hankel transform

1. Introduction

    The problem of reconstructing an object from a set of its projections arises in among other

    fields, computer-aided tomography CAT, radio astronomy, electron microscopy, and spotlight-mode

    synthetic aperture radar. Traditional analytic algorithms of reconstruction for parallel or fan beam

    projections are the convolution backprojection method [1-6], the Fourier method[7] and the Hankel

    transform method [8]. The convolution-backprojection method is fast and widely used. The Fourier

    method is faster than the convolution-backprojection method, but the reconstruction quality is often not

    satisfactory because interpolation in the Fourier domain easily introduces artifacts[9]. Recently a

    Fourier method has been introduced for the fan beam case [10] by our group. It has no interpolation

in the Fourier domain, so that good reconstruction quality can be obtained. For a special case where

the object to be reconstructed has no high frequencies along the angle;? direction, the Hankel transform

    method can be used and good quality of reconstruction can be achieved by few projections. However,

up to now the Hankel method is only used for parallel beam projections. This paper is extending the

    Hankel method for the fan beam case. We derived a new integral transform pair which plays exactly

the same role in the fan-beam case as the Hankel transform pair in the parallel-beam case.

2. Two integral formulas for fan beam tomography


    The spirit of the Hankel transform method for parallel beam tomography is first to obtain the

    integral transform from the object function f(r,?) to the projections P?(u) and the inverse transform

    from the projections P?(u) to the object function f(r,?), then expand these into Fourier series. We

     follow this spirit to obtain the general Hankel Transform pair which can be applied to fan beam


    From Ref. [11] the method of back projection for the case of equally spaced collinear

    detectors is available 2;;;;; (1) f(r,?,?;; 1;;R'??sh(s'-s)dsd?; ;;; U2 0 -; where

     ;; 1 (2) h(t)= |?|exp(j?t)d? 4;-;;

    U=1+? sin(?;?,) s'=r cos(?;?,;U (3)

    D (4) R'??s)= R??s) D2+s2

    and f(r,?) is the image function to be reconstructed in polar coordinates (r,?,( R?(s) denotes a fan

    projection as show in Fig.1 where s is the distance along the straight line corresponding to the detector


     bank and;? is the rotation angle ,;?=D , where D is the distance from the source point;? to the origin



     ?;objection fuction.

     s D;?;

    (r,? ) u



    Fig 1. Fan beam geometry.


    Define the one dimensional Fourier transform of R'??s) as

     (5) G'???,??Fs(R'??s))



    (.)=;;(.)exp(-j?s)ds; sF-1?(.)? 1;;(.)exp(j?s)ds ?Fs


Combining (1) , (2) and (5) and changing the order of the integral ,we have ;; 2;; 1;;; f(r,?,?;??,;4;,; U2

    0 ?;


    Equation (4), (5) and (6) gives an integral transform from R??s) to f(r,?,, now we derive the other

     integral transform from f(r,?, to R??s). The parallel beam projection P??u) is[12] ;; R ??,; P?(u)=;;;;; f(r,?,??u-r cos(?;?,,rdrd?; ??, r=0

     the fan beam projection R??s)[11] is

     R??s)=P?(u)|u=D sin(??s));???!??s) (8) where

    ??tan-1(s/D) (9) considering u=D sin?;


from [13] we have

    1 ??u-r cos(?;?,,?;(11) Ucos?;

    where U, s' are defined in equation (3) and

    D cos??;?:;,;


    Using equations (8) (9) and (11) ,equation (5) becomes




     - ;;;;;; R ;;1 U cos?;?;-;;;??,







    G'?? 1 exp(j?s')d??d?


    ?s-s') G'??,? D2+s2R?(s)exp(-j?s)ds

    ?s-s')rdrd?exp;j?s]ds ?cos? f(r,?


    R ;; 1 f(r,?exp[-j?s']rdrd??;?:~,;U



     where f(r,?,?, if r>R. Equation (6) and (13) are two basic integral formulas. (6) can be used for

image reconstruction. (13) can be used to calculate the projection data from the image function as

required for simulation. These two integral formulas can be expanded using Fourier series, ;; ?:?,;f(r,?,?;;fm(r)exp(j m;?,; m=-;; ;; (15a) R'??s)=;;R'm(s)exp(j m;?,; m=-;; 2;; 1 (15b) R'm(s)= R';?s)exp(-j m?d?2;? 0


    ?:?a) G'?(?,?;;G'm(?,exp(j m;?,; m=-;;

     combining definition (5) and equations (15) (16) gives,


    (16b) G'm(?,?;;R'm(s)exp(-j s;?,ds -;;

     Now we obtain the new integral transform pair

    ;; ;1 (17a) fm(r)=2 exp(j m2)G'm(?Hm(r,?,:?d??;



     ;;; r sin(?,;1;; 1 exp[j(m;?;;?:?b) ?,?d?Hm(r,?,?;U 2;U?;


and a ; (18a) fm(r)Im(r,?,rdr G'm(?)=2; exp(-j m;,



     ;;; r sin(?,;1;; 1 exp[-j(m;?;;?:(b) ?,?d?Im(r,?,?U 2;;?U



In these two formulas the substitution;?;???!;;; has been used and U=1+?cos(?,(;


    (4) (15b) (16b) (17) (14)

can be used as an algorithm for the fan beam image reconstruction. The equation (18) can be used

deriving the fan beam sampling theorem as in Ref.[12].

3. The connection between the new integral transform and the

Hankel transform

    In this section the connection between the new integral transform and the Hankel transform is

derived. This connection gives an alternative to obtain the new integral transform pair.

    One algorithm for the parallel-beam image reconstruction is the Hankel transform method [8]




     a } (20) Fm(?)=2

    ; exp(-j

    where Fm(?, is defined by fm(m;,



















    (22) F???,?;;Fm(?,exp(j m;?,; m=-;;

    Jm is the first class Bessel function

     ;;; 1 ?;~,;Jm(r;?,?;



    Hm(r,?, ---> Jm(r;?,;

    Im(r,?, ---> Jm(r;?,;

?his means that the Hankel transform is a special case ( where;? --->0) of the new integral transform.

    The rebinning relation between the fan beam projection and the parallel projection is [11] s (24) u (25)



    fm(r)= exp(j m;,Fm(?Jm(rd???



    exp[+j(m?;r? sin(?,,?d?2;

    First we note that if??D --->0 then

    R?s)=P?u)|u=s D/ D2+s2,?=?tan-1(D)

    P?u)=R?(s)|s=u D/ D2-u2,?=?-sin-1(D)

    combining equation (24) and (25) we obtain Using equations (5) and (21),


     D2+s2 Du u -1 ?Fu{ (26) D D






    D sD S ?s{ F(27) D2+s2 D





    Substituting (16-a) into (26) and comparing with (22), we obtain


    u D uD -1 (28) e } ?Fu{ D2-u2 D2-u2

    Substituting (22) into (27) and comparing with (16-a) , we obtain


    s D sD (29) e } ?Fs{ D2+s2 D2+s2

    Substituting (28) into (19) and comparing with (17-a) we obtain


     s sD -1 (30) } D2+s2 D +

    Substituting (20) to (29) and comparing with (18-a) we obtain


    s D sD (31) } F-:?'{Jm(r?')}|u= ?Fs{e D2+s2u D2+s2

    The two formulas above are connections between the Bessel functions Jm and the new

functions Hm and Im. We have obtained the two function Hm and Im from formulas (6) and (13). This

can also be done directly from (30) and (31). Practically we first obtained (17-b) this way. See


    4. Implementation of this algorithm 1) The implementation of (4)

    calculates the weighted projections:


    sF?'{G'??'}|s= ,???;sin-1D}

    -1 uF?'{F??')}|u= ,???tan-1(D)}

    -j m sin-1(D) sF?'{G'm?'}|s=

    j m tan-1(D) -uF?'{Fm(?')}|u=

    -j m tan-1(D) D2 Fs{e ?2 s2uF?'{Jm(r?')|?'|}|u=

    j m tan-1(D) D1 (32) R'(?i , sk)= R(?i , sk) D12+sk2

    where R(?i,sk) --data of fan-beam projections,;?i=i*;?,;;??;;;??p+1), (Np+1) --number of

projections, which is chosen as a power of 2 so that the FFT algorithm can be easily used. i in [0,Np],

    sk=k;s,;;s=(D1/(D1+D2));t, D1--distance from the source;? to the centre O, D2-- distance from the centre to the detector plane,;;t detector distance, k in [-Ns2,Ns2], Ns2= Ns div 2, div --integer division,

Ns--detector number. See fig 2.


     ; s detector Smax umax plane ?;?;

    D1 D2


    Fig. 2. Object, projections and detector distance.

    2) The implementation of (15-b) calculates the components of the Fourier series. First we define


    (33) MODm-->n


    (34) IMODm-->n

as operators, which change the integer value m to n according to n=m mod N (35) and m if m in [0,N2-1] n={ (36)

    m-N if m in [N2,N-1]

respectively, where N2=N div 2. Define the FFT and IFFT operators for any series Xm,


Report this document

For any questions or suggestions please email