home
PDF

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

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 diﬀerent 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 ﬁrst 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 times that of the image resulting from using the all-at-once method. Here is the number of projections.

This report is the result of investigation made to determine the reason for this diﬀerence.

### 1 Radon transform introduction

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

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

Any point on the line with speciﬁc and speciﬁc satisfy

Assuming there exist a function deﬁned over the region shown above. The integral of this function over the line is

where is a diﬀerential element of the line

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

Hence for a speciﬁc the above will integrate over the line . The above is the radon transform of over the line . So the radon transform is really the line integral of a

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

Now suppose we have many parallel lines to and we perform the line integral of over each of these lines (since all these lines are parallel to line , then all of them will have the same , but they will have diﬀerent each This will result in many line of the projection line as follows

So, if we do the above over many parallel lines, we obtain many sample points on the projection line , Notice that the projection line 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 speciﬁc angle we obtain a projection vector as shown in the following diagram

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

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

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

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

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

So, assume the ﬁlter frequency spectrum is given by the vector , (this is a complex vector, since it is the frequency spectrum of the ﬁlter). Hence the ﬁltered backprojection is given by

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

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

Now we have obtain the spatial representation of the projection after being ﬁltered. 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

The above is illustrated by the following diagram

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 ﬁlter ram-lak used to obtain the ﬁltered backprojection

Now that the ﬁltered projection is obtained, then backprojection is done. Notice that the ﬁltering, when we talk about ﬁltered 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 ﬁrst, then apply ﬁltering on the resulting backprojection image.

The ﬁltered projection is , where the tilde indicates this is a ﬁltered projection. Now that we have calculated , we can obtain the backprojection, which will be a 2D function. Assume the original function (which we do not know in practice, was , then let the backprojection function be . Hence

But , hence the above becomes

What the above is saying, is that to ﬁnd at some position, the angle is changed to cover all the angles from zero to , and for each angle in this interval, the function is evaluated at . The sum of all these gives at that values. We do these for each value in the region to obtain all the values of . Assume there are projections made. Hence angles (since each projection corresponds to one angle). Then the discrete version of the above integral becomes

Since there are angles, then we divide by to obtain each speciﬁc angle in the range. Hence

And also we see that So the sum becomes

 (1)

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

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 ﬁlled radially. This diagram from Kak illustrates this:

Then ﬁltering 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 ﬁnal 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 diﬀerent results depending if the call is made in all-at-once vs. one-at-a-time.

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

If we call the sum as then the above becomes

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

 (2)

I am not sure now why Matlab divides by an extra . I think this might be because if the angles given at from to , then there is a double counting involved (since the image can be fully constructed from to , and hence it assumes the angles given have full range of , and so to compensate for doubling the intensity of the image, it divides by an extra 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 angles, and we call iradon 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, , and so we write

Now we add all the above and obtain

 (3)

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

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

Therefore, we see that

Hence the one-at-time result is 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 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.