home
PDF

Note on radon and iradon transforms and Matlab's iradon on the all-at-once call vs. the one-at-time call

Nasser M. Abbasi

July 17, 2008 page compiled on June 30, 2015 at 11:43pm
This note was motivated by a strange result found when using Matlab iradon.

In matlab, when using iradon to obtain a backprojection image from a set of projections, a different result is obtained depending if one calls iradon() passing it the set of projections all at once vs. if one calls iradon() once for each projection, and then sum the resulting set of backprojections.

Calling the first method above all-at-once method and the method as the one-at-time method, then it is found that using the one-at-time method resulted in an image whose intensity levels is 2N  times that of the image resulting from using the all-at-once method. Here N  is the number of projections.

This report is the result of investigation made to determine the reason for this difference.

1 Radon transform introduction

PIC

The equation of the line L  can be written in 2 ways. The standard way is

y = mx  + b

Where m  is the slope and b  is the intercept. It can also be written in terms of the parameters p  and    𝜃  as

L (p,𝜃) = {(x,y ) : x cos𝜃 + y sin 𝜃 = p}

Any point (x, y)  on the line L  with specific p  and specific 𝜃  satisfy xcos 𝜃 + ysin𝜃 = p

Assuming there exist a function f (x,y)  defined over the region shown above. The integral of this function over the line L (p,𝜃)  is

∫

   f (x,y)ds
 L

where ds  is a differential element of the line

PIC

It is simpler to express the above integral in terms of x  and y  . To do that, a trick is used with the help of the delta function. The above integral can be written as

  ∫∞   ∫∞
          f (x,y) δ (xcos 𝜃 + ysin𝜃 − p) dxdy

y=−∞ x=−∞

Hence for a specific p,𝜃  the above will integrate f (x, y)  over the line L (p,𝜃)  . The above is the radon transform of f (x,y)  over the line L (p,𝜃)  . So the radon transform is really the line integral of a f (x,y)

The result of the above radon transform is one numerical value. It is the line integral value. We can imagine a projection line into which we accumulate the result of these line integrals as follows

PIC

Now suppose we have many parallel lines to L  and we perform the line integral of f (x, y)  over each of these lines (since all these lines are parallel to line L  , then all of them will have the same 𝜃  , but they will have different p  each).  This will result in many line of the projection line as follows

PIC

So, if we do the above over many parallel lines, we obtain many sample points on the projection line z  , Notice that the projection line z  in the above diagram is some arbitrary line drawn just to collect the result of the line integrals into. It represents a detector which collects the results of each line integral. If we collect many line integrals to cover the whole region. Hence for each specific angle 𝜃i  we obtain a projection vector gi  as shown in the following diagram

PIC

The projection shown above is a discrete function. It is a function of p  and 𝜃  . Hence we have g(p,𝜃 )  . But for the same angle 𝜃  , g  is a function of p  , So some books write g𝜃 (p)  .

Notice that z  is parallel to p  , and can be called the axis of the projection.

Once the projection g𝜃 (p)  is obtained, then it is converted to the (discrete) Fourier domain using FFT. The discrete Fourier transform is

      N∑− 1      (  2 πi   )
Gk  =    gn exp  − ----kn     k = 0,1,⋅⋅⋅,N  − 1
      n=0           N

Hence we obtain the vector G  which is the discrete Fourier transform of the projection g  , This is illustrated in the following diagram. Notice that the numbers G
  k  are complex numbers and hence have phase and magnitude. In the implementation of the discrete Fourier transform, the FFT is used for performance.

PIC

So, why do we do the above? The reason is to filter the projection data. Filtering the projection produces better backprojection (sharper) than without filtering (more blurred). It is easier to apply filtering in the frequency domain than in the spatial domain (multiplication vs convolution). Now that the FFT is done and G  obtained, a filter is selected. Consider the ram-Lak filter. This filter has a frequency response as follows1

PIC

So, assume the filter frequency spectrum is given by the vector H  , (this is a complex vector, since it is the frequency spectrum of the filter). Hence the filtered backprojection is given by

˜G =  GH

Where I use the tilde symbol to represent a filtered frequency response.

Now that we filtered the projection, we need to return back to the spatial domain. Hence we obtain  ˜g  by inverse discrete fourier transform of ˜G  and this is given by

        N−1       (       )
     -1-∑   ˜       2πi-
˜gn = N     Gk exp    N  kn      n = 0,1,⋅⋅⋅ ,N − 1
        k=0

Now we have obtain the spatial representation ˜g  of the projection g  after being filtered. However, this contains both complex and real components (since it is complex valued as result of the IFFT). Then we need to take its real part

˜g = real(˜g)

The above is illustrated by the following diagram

PIC

The above can also be represented using plots of the spatial and frequency spectrum of the above vectors. For these plots, I used a projection signal made up from some simple function of x. The plots shown are the actual FFT and IFFT and the filter ram-lak used to obtain the filtered backprojection

PIC

Now that the filtered projection is obtained, then backprojection is done. Notice that the filtering, when we talk about filtered backprojection, is carried on the projection itself, and not on the backprojection. I am not sure if it is possible to do backprojection on the projection first, then apply filtering on the resulting backprojection image.

The filtered projection is ˜g𝜃 (p)  , where the tilde indicates this is a filtered projection. Now that we have calculated ˜g𝜃 (p)  , we can obtain the backprojection, which will be a 2D function. Assume the original function (which we do not know in practice, was f (x, y)  , then let the backprojection function be f  (x, y)
 B  . Hence

            ∫π
fB (x,y ) =   ˜g𝜃 (p)d𝜃
           𝜃=0

But p =  xcos 𝜃 + y sin𝜃  , hence the above becomes

           ∫π

fB (x, y) =   ˜g𝜃 (xcos 𝜃 + ysin𝜃) d𝜃
          𝜃=0

What the above is saying, is that to find fB (x,y)  at some x,y  position, the angle 𝜃  is changed to cover all the angles from zero to 1800   , and for each angle in this interval, the function ˜g (p)
 𝜃  is evaluated at p = x cos𝜃 + y sin 𝜃  . The sum of all these gives fB  at that x,y  values. We do these for each x,y  value in the region to obtain all the values of fB  . Assume there are N  projections made. Hence N  angles (since each projection corresponds to one angle). Then the discrete version of the above integral becomes

           {                           }
             N∑− 1
fB (x,y) =      ˜g𝜃m (x cos𝜃m + y sin 𝜃m)   Δ𝜃
             m=0

Since there are N  angles, then we divide π  by N  to obtain each specific angle in the range. Hence

𝜃  =  π-m
 m    N

And also we see that Δ 𝜃 = Nπ  So the sum becomes

            π N∑−1    (     ( π   )        ( π   ))
fB (x,y) = ---   ˜g𝜃m  x cos  --m   + y sin  ---m
           N  m=0            N             N
(1)

And the above is the equation used for the backprojection. The following diagram illustrates the above.

PIC

Notice that in the above, backprojection is carried out in the spatial domain. This is how the Matlab iradon does it. Other way to do backprojection is to stay in the frequency domain by using the central slice theorem. This is done as follows: One the fourier transform of the projection is found, it is moved to the 2D plane which represents the 2D fourier transform of the image being reconstructed by the backprojection. It goes into a radial line through the center of the 2D fourier transform at the same angle 𝜃  of the projection. This is done for each projection. The result is a 2D fourier transform of the image. However, it is filled radially. This diagram from Kak illustrates this:

PIC

Then filtering is done now in the frequency domain. Next, gridding is carried out to transform the above to a Cartesian coordinates since IFFT works on these only (by interpolation). Then IFFT is done to obtain the final backprojection image.

But since the purpose of this note is to talk about Matlab iradon(), then we will use (1) as the method of backprojection.

2 Matlab iradon and the all-at-once vs the one-at-time

Now we start the investigation on why when using Matlab to determine backprojection from a set of projections it gives different results depending if the call is made in all-at-once vs. one-at-a-time.

Let consider the all-at-once first. From (1) we see that the backprojection image is

            π N∑−1    (     ( π   )        ( π   ))
fB (x,y) = N--   ˜g𝜃m  x cos  N-m   + y sin  N--m
              m=0

If we call the sum as img,  then the above becomes

           π
fB (x, y) =--img
           N

Notice however, that matlab iradon.m does the following at the end: (type iradon.m to see)

      π
fB = ----img
     2N
(2)

I am not sure now why Matlab divides by an extra 2  . I think this might be because if the angles given at from 0  to     0
360   , then there is a double counting involved (since the image can be fully constructed from 0  to     0
180   , and hence it assumes the angles given have full range of 2π  , and so to compensate for doubling the intensity of the image, it divides by an extra 2  at the end.

So, now that (2) gives the backprojection in the all-at-once method, let consider the one-at-time method. Supposed again we have N  angles, and we call iradon N  times. But each time we call iradon with one angle, we have to pass the same angle twice, and the same projection twice, and divide the result by 2. From help:

  Compute the backprojection of a single projection vector. The IRADON
      syntax does not allow you to do this directly, because if THETA is a
      scalar it is treated as an increment.  You can accomplish the task by
      passing in two copies of the projection vector and then dividing the
      result by 2.

So, in each call now, N  = 2  , and so we write

pict

Now we add all the above and obtain

                         π
fB1 + fB2 + ⋅⋅⋅ + fBN =  --(img1  + img2 + ⋅⋅ ⋅ + imgN )
                         8
(3)

Now, img1  + img2 + ⋅⋅⋅ + imgN  must be twice the img  from the call to iradon in the once-at-time method, since using the same assumption as Matlab, where it assumed 0  to 3600   range and not 0  to    1800   , we will double add the intensity. Hence (3) becomes

pict

Compare (4), which is the one-at-time, with (2), which is the all-at-once, we see that

pict

Therefore, we see that

pict

Hence the one-at-time result is N2-   times the all-at-once. And this is what was found.

Therefore, when using iradon() to obtain a backprojection in the one-at-time method, we obtain N2-   times the image that would have been obtained if using the all-at-once method.

3 References

1.
http://www.comap.com/product/samples/UMM794.pdf
2.
Principles of computerized tomographic imaging by Kak and Slaney
3.
Matlab radon and iradon help.