Introduction

The present invention regards computer tomography and methods of evaluating to¬ mography data. The method of the invention can be applied to all kinds of tomo¬ graphic imagmg, whether registered by X-ray or radioactive absorption, by magnetic resonance (MR) or positron emission tomography (PET). In all those cases are ob¬ tained experimental values in the form of line mtegrals along multitudes of straight lmes The typical CT measurement consists of mounting an X-ray tube on a ro- tatable circular πng and a set of detectors m a row opposite the tube.

The invention therefore regards a computer tomography method of the land recited in the preamble of claim 1

In computerized tomographic imagmg (CT) the mverse problem consists of deter¬ mining a function from its line mtegrals We examine three different types of meas¬ ured data sets, with diffeient scanmng geometries, namely two versions of parallel scanning (smogram. linogram) and fan-beam scanning (even cone-beam scanning) We will study the Filtered Backprojecnon Method (FBM) which is one of the most well-known and widely used methods The reason for the widh use of the filtered backpiojection algorithm is die high image quality.

A. seπous problem m prior an evaluation of tomographic data in order to obtain an image is the heavy computational complexity for the back-projection operator, which is of the order N\ It is the general object of the invention to obtain an im¬ proved method of transforrmng tomographic data mto images, which is less complex and can give improvement either in computanonal speed or in an increased number of voxels/pixels in an image

This object and other objects and advantages are obtained, according to the inven¬ tion, by means of a computer tomography method of the kind recited, which is char¬ acterized by the steps recited in claim 1 Advantageous embodiments are recited in dependent claims

In conventional backprojection each pixel value is calculated separately, which means that (9(N) operations are required for each pixel. Then the total number of operations will be of the order (3(N ^J ). Instead of directly doing the N ³ additions the basic idea behind the new algoπthm is to divide the summation into partial sums, which then are used to form new partial sums m a dyadic and hierarchial way, see Figure 1 The reason for this procedure is the benefit in computational reduction due to the fact that a partial sum has approximately the same value for many neighbour¬ ing points Thus each value of a partial sum is calculated once, but used several times when forming a new partial sum from two partial sums in a preceding level of the hierarchy

Although in the present description, the invention is exemplified with forming hier- achial partial sums m a dyadic way, which is presently believed to be the best mode, it is also possible to sum more than two groups, namely, each time a number which is a prime number It is even not necessary to use the same prime number in each stage

In a preferred embodiment, families of lines are combmed two and two, indicating that the total number of line families is a power of two. With a number of pix¬ els/voxels of 512 x 512 an an equal number of families, i.e. angles of registration, this condition is fulfilled. However, without using the full advantages of the inven¬ tion, it is possible m the partial sum calculation to combine three (or five, or any prime number) space distributions by making sums for those points where the lines for the largest angular difference mtersect, and adding interpolation values for the further distributions.

It is also possible to add line integrals for intermediate lines, for which no measure¬ ment was made, i.e by sernng intermediate (non-measured) line mtegrals to the mean value of the adjacent lmes, for which a measurement was made.

Although the method is most easily understood for a case where all the line families belong to the same plane, which is the most common configuration in practice, and which is exemplified as xhe fan-type, the smogram and the linogram, it may also be adapted to the cone-beam case, where radianon from a source fills a cone, and de¬ tection is made in a two-dimensional array of detectors. Each detector men corre¬ sponds to a line through the field swept, and the contraption and the object are mu¬ tually rotated As shall be shown, a considerable saving in calculations may be ob¬ tained also this case

In a special case, where each family of lmes consists of Lines which are mutually - parallel, it is possible, in a very favourable way, to obtain values for points which are equidistant in a rectangular screen pattern. This is obtained by making me two last line families, for which sums at intersections are calculated, be at mutually right angles Thereby , in principle, the last repeated second step is performed by two combined gro ps, symmetrized in relaπon to a direction at right angles to an angle belonging to a family of lmes which is at an extreme end of said enumerable order, thus obtaining a set of sum values which correspond to pixels in a cartesian system of screen pattern

List of Figures

1 The sumπαac on process for N = 3 21

2 Division of the (.-: ^* , H-plane ^f ° ^r I'{ , ι. k'. f) 21

3 Division of che O-plane for I'(l. ι, k', l') 22

4 Division of che (k", l ")-plaae for I'{j, ι, k l") 23

5 Cone-beam geometry 24

6 Points of mcerseccion for the partial sum [(j. i. k. l. m) . . . 2 ..->

2 The filtered backprojection algorithm for

sinogram

In the parallel scanning geometry a cross-section of the human body is scanned by a set of equally spaced parallel X-ray beams for a number of equally distributed directions. Thus the recorded detector values correspond to a discrete set of values from the Radon transform. The (n-dimensional) Radon transform of a function / on R ⁿ is defined by

(Rf)(θ, s) = I f(X) dX,

X θ=s which is the integral of / over the hyperplane whose normal is θ <E S ^{n~ l} , the unit sphere in R ⁿ , and with (signed) distance s € R ^l from the origin. The term sinogram is chosen since an attenuation function consisting of a (5-function, Dirac function, gives a Radon transform supported on a sinusoid in parameter space, i.e., (θ, s)-space.

Before giving the inversion formula we will describe the backprojection. The (n-dimensional) backprojection operator R ⁱ maps a function g on S ^{n~ i} x R ¹ into the set of functions on R ⁿ . The backprojection operator is defined by

(Rtg) (X) = jf ^ g(θ, X ■ θ) dθ, X E R ⁿ .

R ⁱ is the dual operator of R. Hence

/ (Rf)(θ, s) g(θ, s) dsdθ = I f(X) (R ⁱ g)(X) dX.

The filtered backprojection algorithm is a numerical implementation of

the inversion formula

/ = I (2T) ¹ -" R ^u H ⁿ - ^l {Rf) ⁱⁿ - ^y) .

where H denotes the Hubert transform and (-) ^(n_1) denotes the (n-l)st derivative both with respect to the variable s.

3 Fast backprojection for sinogram

In this section we will develope a method for backprojection. requiring only 0(N ² \ogN) operations. This method can be used on the (n-dimensional) backprojection operator but in order to describe the main idea it suffices to study the case n — 2. Then

(R ⁱ g)(X)= f g(Θ.X-θ)dθ,

where X = 6 S ^l , or. with slightly different notation,

(R ⁱ g)(x,y)= / g(υ,xcosυ + ysinυ) dυ.

Strictly we should write g(θ(υ),xcosυ -f ysinv). We assume that g{θ,s) satisfies the condition g(—θ, —s) = g(θ, s) which implies that g(v -f- π, — s) = g(v, s). Then

( Rg) (x, y) = 2 / g(v, x cos v + y sin v) dv.

Let g(υ, s) be sampled at the points g(v _m , s _n ) where

N — 2 ^J , J positive integer. In order to compute the backprojection, we need to calculate the sums

for A: ² + - ² < ^ _j -, where A; and I are integers.

Since ø is defined on a rectangular grid we might use linear interpolation for the variable s in the sum of formula 3.1. We choose instead the following procedure.

First extend the definition of g by increasing the number of sampling points in the s-coordinate with a factor e, using linear interpolation, and introduce the notation

for m = 0, ..., N - 1 and n = -≤f , ..., . Then, take

7 ' ] ) ~ Σ G(m, [e(fc cos u _m + . sin u _m )]).,

where [α] denotes α rounded to the nearest integer.

Now, instead of directly doing the N ³ additions we will divide the sum¬ mation into partial sums, which then are used to form new partial sums in a hierarchial way. The reason for this procedure is the benefit in computational

reduction, due to the fact that a partial sum has the same value for many neighbouring points. First we define i-2->-l

I(j.i,k,l)= ∑ G(m,[e(kcosυ _m ÷lsi v _m )]) m=(ι-l)-2^ for 0 < j < J, 1 < i < 2 ^J~j . Then we have the following simple relation

I(j,i,k.l) = I(j-l,2i-l,kJ) + I(j-l,2i,k,l).

In the following description we will show that the number of points (k,l) where the value of the partial sum I(j, i. k. I) has to be calculated is approx¬ imately twice as large as the number of points where the value of a partial sum I(j — l,i,k,l) has to be calculated. Hence, since the number of partial sums with different i- values at level j is N/2 ³ , the total number of summa¬ tion points for all partial sums at level j is approximately the same as for level j — I and since there are J = ² log N j-levels we conclude that the total number of operations for all levels will be proportional to N ² lo N.

The variable i just stands for a rotation about the origin. Thus it will be convenient to make the following definition

2'-l

I'(j,i,k .l') = ∑ G(m + (i- 1) ^■ 2> .[e(k ^* cosυ _m + T s v _m ))). m=0

Then we have

I(j,i,k,l)

(3.2)

= I ^* (j,i,kcosv(i-.i)2i + - ^' sinι/( _i _ _i ), _I Jcosu(j_ι.)2J — A ^; sinι/(,_ι)2 ) and

By definition we have

We now see that / ^* (0. i. k{, t[) = I"(0. i, A-..], l ₂ ^* ) if the points (k{. l\) and (JfcJ, l ₂ *) are located in the same strip, see Figure 2. Now consider the next level, =l. We have

r(l, ι. k', I') = C(2ι - 2, [ek'\) + G(2ι - 1. [e(k ^* cos υ _x ÷ I' sin _Vi )}),

and we conclude that 7 ^* (1. i. k ^* . / ^* ) = 7 ^* (l, ι, k*_. I*) if (k{, / ^* ) and (Jfc.]. .J) are located in the same region, see Figure 3. Let

e

A ;- « cos vι ^■ +- 1 r . sm vi = - ^ e

According to Figure 3 it will be necessary to calculate 7 ^* (1, i. k* , I*) only for a and b integers. Therefore we mtroduce, general A; ^* and I ^* will not be integers. ϊ(l, ι, a, b) = I'(l. ι, k ^* , r) where a — ek" b = e(k' cos υ - / ^* sin Vι)

Hence

1(1, i, a. b) = G(2ι - 2. a) + G(2ι - 1, b) Can we do a similar division of the (A; ^* , / ^* )-plane for I"(j, ι, k*, l") , > 1?

Let

* ^* =

k _. ^• cosuoj.i + , 1 ;• sm ^• t^-i = - ^b e If two points (k{. l ) and (k ₂ , I ^* ,) are in the same region, see Figure 4, we want r(j,i,kuι) = ι*(j,ι.ti.r ₂ ). .

\(k{ - k;) cos^ _m + (l[ - .;) smt/ _TO | < - (3.4) e for = 0, ..., 2 ^J - 1. In Figure 4 we see that all regions approximately satisfy condition (3.4) if u ₂₎ _ι < , which corresponds to j < J — 1. In the same way as for j = l we introduce

ϊ(j,ι,a,b) = r(j.i,k',n (3.5)

where a = ek ^* b = e( ; ^* cosu ₂ j_i -)- I'sinvj-i) By equations (3.3) and (3.5) we see that we have the following relation

- l,2t-l,a,b _l ) + ϊ(3-l,2ι,a _l ,b) (3.6)

where

and

By symmetry we have

ϊ{j,i,b,a) « /( - l,2i - 1,6,0 + 7 - 1.2ϊ,δι,o)

/(j, i, -α, -6) « /(j - 1, 2t - 1, -α, -0ι) + /(j/ - 1, 2ι, -oi, -6)

and ϊ(j, i. a. a) z ϊ(j — l, 2i — l, a, aι) -r ϊ(j — 1.2? ^' , a _x . a)

where

which means that we only have to compute αi and b for 0 < a < \b\.

The procedure suggested by this discussion of the partial sums leads to the following algorithm. We start by adding "extra" sample values by linear interpolation. Then we form partial sums of adjacent pairwise projections ( A projection m consists of all sample points G(m, n) where rn is fixed. ). The next step will be to form partial sums of four projections by adding two adjacent sums from the previous step. Thereafter we will subsequently form partial sums of 8, 16,32, ...projections. The iteration stops when full backprojection is obtained, that is, when the sum consists of all projections. A partial sum is only determined in points where two lines intersect whose angles are v ₀ and v _2> -ι, with respect to the (A; ^* , / ^* )-plane, and whose per¬ pendicular distances to the origin are ^ and ^ respectively, n _t and π ₂ integers. If e > 2 the process is stopped before we obtain /(J — 1, i, a, b). If the total number of points where the partial sums, at level j. are calculated is of the same order as the number of pixels, the process is stopped and the pixels are interpolated directly from the values of the partial sums at level j — 1. The pixel values are then calculated from equation (3.2). When e = 1 we calculate all the partial sums up to /(J — 1, i, α, b). Although the number of points where I(J — l, i, a, b) is calculated is of the same order as the num¬ ber of pixels they do not coincide. This is caused by U j-i --, = f - f •

Therefore we have to calculate the pixel values by equation (3.2) which is unnecessary if we instead do the following. We want to calculate

.V/2

∑ G(m, [k cos v _m + I sin v _m ]) (3.9) m=0 and

;V

53 G(m. [k cos v _m + l smv _m ])- (3.10) m=N/2

But now v _N / — VQ = V _N — v^/ ₂ = so the points where we calculate the sums are in fact the pixels. One problem is that each sum, (3.9) and (3.10), consist of 2 ^{;_ i} -I- 1 elements. However in order for the summation process to work, the number of elements in a sum has to be a power of two. Therefore we calculate ∑^Q ¹ , ∑^ ^' _N μ ₊ ∑mϊv/2 ^and ∑ΪLΛN/ +I - ^Then

and

Note that since the partial sums ∑ i ₀ and ∑ ,= / ₂ both contain the pro¬ jection rn = N/2 and the projection m = 0 ~ m = N, we do the following. We assign the projections the values G(N/2, π) = 0 and G(0, π) = 0 in one of the partial sums.

The total number of different values for ^ or - ^l = , k ² + I ² < ^, is

If and i are fixed then the number of points where we calculate ϊ(j, i, a, b) is approximately the area of the circle, ^E -, divided by the area of one region

( ^rhomb), . c_. sin' V.n _j _ , , i.e., e ²²¹ 4 - sin ϊΛv, -ι - Hence the total number of points a jt. i level l j ^■ , i 1 = ^e 2 τ ² V ² "" -TI _' When j is small, ( =l, 2, 3.), the overlapping of the circle by rhombs is not so efficient but then -^ = ^, f , . Hence the total number of points for all levels is e"τr

< -— A- log N ÷ N ² (e - l)

for e > 2 and

7T

< Λ , ^' -2 log N

for e = 1.

4 The filtered backprojection algorithm for

linogram

In 1987 a new reconstruction algorithm was developed by Edholm [Edh87]. It is called the linogram method and is similar to the direct Fourier method but it is based on a different scanning geometry. In the parallel scanning geometry for sinogram, the projection lines. X-rays, are taken for a num¬ ber of equally distributed directions which requires the detectors to rotate around the object. The linogram projections actually consist of two sets of projections. In the first linogram the detectors and the transmitters lie equally spaced on fixed lines parallel to the rr-axis and the projection angles are —45° < θ < 45°. Instead of equally distributed directions, Δ0=constant, we now get Δ tan f?=constant. In the second linogram the

detectors and transmitters lie equally spaced on fixed lines parallel to the y-axis and 45° < θ < 135°,(Δcot(9=constant).

Instead of the Radon transform it now will be natural to study the two linogram operators Ri and R ₂ .

and

The term linogram is chosen since an attenuation function consisting of a (5-function gives a i? _r and β ₂ -transform supported on a line in parameter space, i.e.. (u, t/)-space.

There exist dual operators, backprojection operators, R\ and R ₂ iven by

and

f∞ -OO

I

(R _l f)(u,υ)g(u.v)dudv= l f(x.y)(R[g](x,y)dxdy -OO J —OC '

_! _ ^{R2f U} ' ^V)9{U ' ^{V) dUdV =} S-oo i-oa ^(X ' ^V) (^) ^(X ' ^{V) dXdy} -

The filtered backprojection algorithm for linogram is a numerical imple¬ mentation of the inversion formula

where H denotes the Hubert transform and ( ■ ) ⁽¹⁾ denotes the derivative both with respect to the variable u.

5 Fast backprojection for linogram

Here the procedure of summation for the backprojection is analogous to the one used for sinogram. It suffices to study the backprojection R and projection data from Ri since R ₂ and R corresponds to a rotation about the origin of 90 ^c degrees. We want to determine (x, y) in the pixels (x, y) = (τ ) ) . N = 2 ^J where J is a positive integer and A:, I integers belonging to the disk k ² + I ² < ^. Then the function g(u. v) has to be sampled in the points g{u _m , v _n ), v _n — n- n — — and u _m = m-— , m integer, for all g(u _m , υ _n ) where the corresponding projection line intersect the unit circle. Therefore the number of parallel projection lines depends on the projection angle. For instance if n = 0 we have \m\ < y but if n = y we have \m\ < 2^-. We will divide the summation.

into partial sums in the same way as we did for sinogram. In order for the summation process to work, the number of elements in a sum has to be a power of two. Therefore we start by splitting the summation in the following

ay.

+ v _n -i) _t v _n )

(5.11)

N/2 -i

Then we divide the sum ∑ into partial sums ( The procedure for ∑ is n=l n=-N/2 similar ). Define

for 0 < j < J— 1.1 < i < 2 ^J~J~1 . The procedure suggested by this discussion of the partial sums is analogous to the method described in section 3 on page 8. A partial sum I(j, i. k, I) is only determined in points where two projection lines intersect, one projection line corresponds to and the other to υ _{t 2} , . If the projection line, which corresponds to (u _a , U ₍ ,_ _I) . ₂ J+ _I ). intersects the projection line, which corresponds to (uj,, v,. ₂₁ ) , at the point ( ι,a,b _j .ι, . _b ) in the (A .)-plane we will use the notation

ϊ(j, i. a, b) = I( , i, k _Jtl<aιb , l _3ll , _a ,b)-

The projection line x = U _m — υ _n ^■ y is the straight line through the points (m + n, —N/2) and (rn — n, N/2) in the (k, .)-plane. We have r jfc ^» j,ι,a,b , - — _Ω ^u . ₂ ^x ; ^V -l b ^{U f} '- 2 ¹ J ⁾ -- ^J 1^ ¹

' ₇ ,t,α,6 — 2J-1

Since we only are interested in points on the disc k ² 4- 1 ² < ^ N- ² we get

\b-a\ <2> -1.

We obtain a similar relation for linogram as relation (3.6) for sinogram We have

I(j, i, a, b) « ϊ(j - 1, % - 1. , 6 -- /(j - 1, 2z. α _λ , 6)

where

and

2^- ¹ 2 ^J-1 — 1 2- ^7-1 — 1

6ι = ~ — - a — ^- — 6 = a - ¹ - - (6 — )

2^ — 1 2^-1 2^ — 1

One major advantage of linogram compared to sinogram is that in the equa¬ tions for αi and 6ι the terms which not are integers depend only on the difference b — α Then if we use linear interpolation the interpolation factors will be constant if b = b — a =constant Hence if we introduce the notation

I(j,ι,a.b) = ϊ{j,ι,a,a-rb)

we finallv obtain the simple and very effective summation algorithm

« (1 - c)T(j - 12% - 1 a b - A) -*- cϊ(j - 12; - 1, , 6 - .4 - 1)

+ (1 - c)I(j -l.2ι,a- ^L -A b-A) + cl(j - 1, 2ι, α ^J - + 1, b - A - 1)

where 6- ₃ 7 ₃ 7 = A+c A integer and 0 < c < 1 Since c = 0 for b = 0, ±(2 ^J -1) and by using that we have the same interpolation factors for ±6, we only have to compute for 1 < b < 2 ^J -2 ( Let 1 < b < 2 ³ -2 and b f _j = A + c, A integer and 0 < c < 1 Then -b • ^ = - A - 1 + (1 - c) ) For sinogram,

however, we had to compute a _\ and b _\ for all combinations of (a, b), a > b. Therefore the total number of times αi and 61 have to be computed, at level j, is approximately 2(2 ^J - 1)/V for sinogram. Hence although the number of summation points (a. b) is approximately the same for linogram and smogram the total computational complexity for a,ι and 61 _. is O(N) for linogram but 0(N ² ) for sinogram.

Finally we calculate according to equation (5.11). Let

7ϊ(a, b) = ϊ(J - 1, 1. a, b)

and I2(a, b) for the corresponding partial sum

∑ 9 (k + V-n - l). V- _n ) ,

where a corresponds to the direction ι/_ι and b = a + b corresponds to the direction υ_,v/2- Then

for 0 < I < N/2.

The total computational complexity for the filtered fast backprojection algorithm for the linogram geometry is the same as for the sinogram geome¬ try, i.e., N ² log N.

6 Cone-beam and Fan-beam geometry

The most common scanning geometries in medical tomography are the cone- beam geometry in the three-dimensional case, and the analogous fan-beam geometry in the two-dimensional case. In cone-beam scanning an X-ray source is rotating around the object along a circle while emitting cones of X-rays, which are recorded by a two-dimensional detector. In fan-beam scanning the cone-beams are screened off to a single plane and consequently recorded by a one-dimensional detector. Feldkamp, Davis and Kreiss [Fel84] have developed an approximate 3D-reconstruction algorithm for cone-beam projections. The computational complexity is of the order 0(N ⁴ ) for an image / of N x N x Nx pixels. By using the same technique as we did for the previous backprojection algorithms we will reduce the number of operations to C(Λ ^τ3 lo N).

Next we will describe the cone-beam geometry. Let the position of the X-ray source be denoted by S(λ),

where 0 < λ < 2π and r constant. Hence the X-ray source rotates around the object along a circle in the zy-plane. The detector plane is always per¬ pendicular to the vector (λ), see Figure 5. Hence it will be practical to introduce the plane spanned by the normed vectors p(λ) and z, where p(λ) is the vector in the xy-p\a,ne which is perpendicular to S(λ), see Figure 5.

Then the measured data can be represented by a discrete set of values from

-7( ,p,g)= J f(x)dx.

-.(A.p.f) where L(\.p,q) is the straight line through the points S(X) and D(p,q) == pp(X) + qz. Now the approximate inversion formula can be written as [Fel84]

1

/(£) = A(λ, x) ² G{X, A(λ, x)x ^■ p(X), A(X, x)z) dλ, (6.12)

4τr ² Jo w here

A(X,x) r ² -x-5(λ) and

{l, F.I) ii,F- ^l .i)

G{λ,p.q) = \σ\

_, \S(X)-D(p,q)\ 9( P,Q) [X. σ, q) (λ,P,<?),

i.e., we multiply the projection dataρ(λ,p, q) with a weight function .^. ^, — ΓT and then perform a rampfiltering in the variable p. The final backprojection in (6.12) is computationally the most heavy part of the inversion algorithm. The major difference between this backprojection and all the previous is the weight function A(X,x). Now let g(X,p.q) be sampled at the points g(X _a ,Pb, Qc) where, r > 1 and N = 2 ^J . with J integer,

and we want to estimate f(x) in the points Xk.i.m = (%k, yι, Zm) — ( , j, η^j , k,l and m integers, for k ² + I ² ÷ τn ² < ^- . Thus

ι 2ιV-l ~ -r— _r - ∑ A(X _a , Xk. m) ² G(X _a , (λ _α , XkΛ,m)Xk.l,m ^■ P(λ _β ), A(X _a , Xk,l,m)Zm),

4ττ^ _α -^o w here

A{X _a ,Xk,ι,m) = = A(a,k,l) and

X _k ,ι _,m -p(λ _β ) = Xfcsinλ _α - yι cos λ _α = (x-p){a,k,l).

If we now divide the sum into partial sums in the same way as we did in section 3 we get

53 A(a. k, l G(X _a , A(a, k, l)(x ^■ p)(a, k, l),A(a, k. l)z _m ). α=(t-l)-2J

Now if the interval is small, i.e. j small, we can make the approximation

A(a.k,l) vA(a _JιX .k.l) (6.13)

where a. _jιX = % • 2 ^J — (2 ^j + l)/2. the average angle. Thus t-2->-l

-4(α _jιt , k, I) ² 5Z ^G ( ^λ <" ^A (H ^k > (* ^■ P)( ^Q _; k, I), A(a _jιtl k. l)z _m ). a={x-ιyv

Define i-2-»-l

I(j,i,k,l,m)= G(λ _α ,( -p)(α,A ^; , ,- )- α=(t-l)-2>

Then as before

J( , t, A;, Z, m) = I(j - 1, 2t - 1, A, , ) + /( - 1, 2i, k, I, m).

A partial sum is only determined in points where two projection lines, pro¬ jected on the xy-pl&ne. intersect one projection line coming from S X ₍ ,_ιy _2] ) and the other from -?(λ,. _2J _, ), see Figure 6. Now since .4(λ, x) obviously can't be approximated by a constant for all λ, except x = 0, the iteration process is stopped after a few steps, at level JO _> and the pixels are interpolated directly from the values of the partial sums at level j _Q . We have

~ k, 1)1, A(a _n , _t , k. l)m).

At a first glance the method seems to be far less efficient for this backpro¬ jection than for the previous ones. But if we consider that the level j — I reduces the number of summations by half and j = 2 to approximately a fourth of the original number, and so on, we see that it is almost as efficient as the other. The accuracy of the approximation (6.13) depends strongly on the distance between the point (k, I) and the origin, i.e., But since the inversion formula (6.12) is correct only for points in the plane z = 0 and the artifacts in the other planes z = m increase with the absolute value of rn the region of interest usually is k + I ² -r- m ² « r. For (k, l) fixed point the accuracy of (6.13) depends on the angle difference -I — λ ₍ ,_ _{1) 2} j not on the number of elements. 2 ^] , in the partial sum.

2 ^] - 1 1 -I — λ(j_ι).2r = TΓ — — —

Hence if N — 2 ^J corresponds to that the iteration process is stopped at level jo then 2N corresponds to that the iteration process is stopped at level j ₀ -f- 1

in order to have the same accuracy for the approximation (6.13). Then if the region of interest consists of n ³ pixels the number of operations for the back- projection will be of the order 0(n ² N log N). In the two-dimensional case, fan-beam scanning, the procedure is identical, m = 0. Hence if the region of interest is of size n ² the computational complexity will be 0(nN log N).

References

[Edh87j P.R. Edholm. G.T. Herman. Linograms in image reconstruction from projections IEEE Transactions on medical imaging, Vol. MI- 6. pp.301-307. 1987.

[Fel84] L.A. Feldkamp. L.C. Davis, J.W. Kress, Practical cone beam algo¬ rithms Journal of optical soc. Am., Vcl. A6. pp.612-619, 1984.

[Νat86] F. Natterer, The mathematics of computerized tomography John Wiley & Sons Ltd, Chichester, 1986.

[Radl7] J. Radon. Uber die bestimmung von funktionen durch ihre integral- werte langs gewisser manmgfaltigkeiten Ber. Verh. Sachs. Akad. 69, p.262. 1917.

Instead of using the approximation (6.13) we can even perform the multi¬ plication of the weight function in an iterative way. In order to describe the main idea it suffices to study the two-dimensional case.2 = 0. the extension

to the case z φ 0 is obvious. We introduce the notations x = R-, y )

L{X.p) = L(X,p,Q) and i 2-»-l

I(j,i,x)= ∑ A(X _a ,x) ² G(X _a .A(X _a ,x)x-p{X _a )). α=(»-l)-2J

Let x(j,ι,pι,p ₂ ) denote the point of intersection between the lines 7-(λ _(l _ ₁₎ . ₂ j,2pι/ ^' ) and L(X _l . _2} _ _i ,2p ₂ /N). Then we will use the notation

We obtain a similar relation as (3.6)

Hj,i,Pi,P2) ~Ki ~ l,2i- - 1.2i,d ₂ ,p ₂ ),

where d _\ and d are calculated in a similar way as equations (3.7) and (3.8) for aj and b _x . In general d _x and d will not be integers. In order to obtain the correct value of the weight function we will apply the following procedure. We will only consider the partial sum I{j — l,2ι — l, ι.d ) the procedure for ϊ(j — 1,2., d , p ) is analogous. We start by only consider the term which corresponds to a = (i — 1) • 2° -- 2 ^}~1 — 1, i.e.,

A(X _i ,x _l ) ² G(X _i ,2d ₁ /N)

« A(Xi, Xi) ² ( c _t G(Xά, 2^3/N) -r c ₂ G(X _ά , 2([d ₁ ] + 1)/N))

+c ₂ ^A , ^Xa~ ' A(X _a -, x ₃ ) ² C(λ ₅ , 2{[d _x ] + 1)/N),

^ ^■ \ ^Λ ά^Xi) ^" where £-, = x(j — 1, 2£ — l,pι,dι), x ₂ = x(j - l,2i- l,pι, [c-ι]) _; x ₃ = x(j - 1, 2i - l,pι, [dι] + 1), d _x = [r,] -f c ₂ , 0 < c ₂ < 1 and c _x + c ₂ = 1.

Therefore we will use the following approximation

ϊ(j -l,2ι- l.p^i;) « _Cl ^ ^ /(j - 1,2* - l, _Pl) ft])

+c ₂ ^TTT ^ ^ - 1- 2z - l _)Pll [d _x ] + 1),

Λ^Λύ- 3J where ά = (i — 1) ^• 2 ^J - 2 ^]~1 — 1. In order for the approximation to be accurate for all terms in the partial sum ϊ(j — l,2ι — l,p _x .d ) all weight functions which belong to the partial sum have to satisfy the condition

A(Λ _αl ,x(j - l,2i-l.p _x ,d _γ )) ² _ A(X _a2 .x(j-1.2ι-l,p _x ,d _x )) ² A(X _aX ,x(j- 1.2ι-l.p _x ,d) - ~ A(X _a2 ,x( - l,2ι-l.p d)) ' iov(ι-l)-2 ^] <ol,o2< {ι-l)-2 ⁱ ÷2 ^J - ^l -l<md[d _x ] <d< [d _x ]÷l. The error obtains its maximum when l = (1 — 1) ■ 2 ^J and α2 = (1 — 1) ■ 2 ^J ÷ 2 ^J~l — 1. It suffices to consider the case

A(0,x( - 1, pχ.d _x )) ² _ A(X.x(j - 1, l.pi.di)) ² A(0,x(j-l,l, _Pl ,d))* ~ ^(λ^U- i .p _! ,-.)) ² ' for i = 1 and λ = 2 ^3~l - 1. The points x( - 1. l.pi.di) and xl - 1, l,pι,fi) are located on the fine (0, 2p _x /N) which can be written

Let x(j - 1, l,pι,dι) = x{t _x ) and £(ji - 1, l,p _x ,d) = x(f). Then

(0,x(r ) = ^ and A(0,x(£)) = 7

e have

^ t l + (^)(l-cosA)-^sinΛ _ t ₌ A(0,x( )) ~ t _x i + (!=£_.) (i _ cosλ) - 2*. sin λ ~ *ι (0, (*)) ^'

The number of points of intersection increase with λ. In fact we have l*ι - *l ^κ / - ^Hence 1 1--tt for some constant c. At the first steps of the iteration the error will be extremely small oc 1/N ² and then increase to o 1/N at the final step.

Now we return to the three-dimensional case, z ψ 0. We introduce the notations i-2>-l

I(j, i, x, z) = Z ⁱ⁴ ( ^λ "- ^) ² G{X _a , A(Xa, x) ^■ p(λ _α ), A(X _a , x)z) α=(t-l)-2-' and ϊ(J, P P2,z) = I(j,i,x(j,i,pι,P2),z).

We have the relation

HJ,I,P P ₂ ,Z) ~ ϊ(j - l,2z- l,p _x ,d _x ,z) + ϊ(j - l,2i,d ₂ ,p ,z),

We will use the same approximation as before the only difference is that we have to adjust the value of the variable z.

ϊ(j-l,2i-l,p _x ,[d _x ],z )

where ά = (i - 1) • V + 2 ^3~l - 1, A = (i - 1) • 2 ^] + (2 ^j ~ ^l - l)/2 the average angle and z ^■ A(X,x _x ) = ri ^■ i4(λ,x ₂ ) = 2 ₂ • A(λ,x ₃ ).

In general z _x and z ₂ are not integers but since the weight function and the points of intersection ₂ and x ₃ are independent of z we can perform an interpolation of the variable z in the usual way. The accuracy of the approx¬ imation in the .r-direction depends on the previous condition, i.e.,

A(X,ϊι) _ .4(A _B ,Jι) _and A(X,xι) _ A(X _a ,x ) A(X,x ₂ ) ~ (λ _α , ₂ ) ⁿ A(X,x ₃ ) ~ A(X _a ,x ₃ ) for (i - 1) ^■ 2 ^j < a < (i - 1) ^• 2 ^j + 2 ^"1 - 1.