July 17, 2008 Compiled on May 26, 2022 at 9:34pm

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 \(2 N\) 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 diﬀerence.

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 \(\theta \) as \[ L\left ( p,\theta \right ) =\left \{ \left ( x,y\right ) :x\cos \theta +y\sin \theta =p\right \} \] Any point \(\left ( x,y\right ) \) on the line \(L\) with speciﬁc \(p\) and speciﬁc \(\theta \) satisfy \(x\cos \theta +y\sin \theta =p\)

Assuming there exist a function \(f\left ( x,y\right ) \) deﬁned over the region shown above. The integral of this function over the line \(L\left ( p,\theta \right ) \) is\[ \int _{L}f\left ( x,y\right ) ds \] where \(ds\) is a diﬀerential element of the line

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\[ {\int \limits _{y=-\infty }^{\infty }\int \limits _{x=-\infty }^{\infty }}f\left ( x,y\right ) \ \delta \left ( x\cos \theta +y\sin \theta -p\right ) \ dxdy \] Hence for a speciﬁc \(p,\theta \) the above will integrate \(f\left ( x,y\right ) \) over the line \(L\left ( p,\theta \right ) \). The above is the radon transform of \(f\left ( x,y\right ) \) over the line \(L\left ( p,\theta \right ) \). So the radon transform is really the line integral of a \(f\left ( x,y\right ) \)

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 \(L\) and we perform the line integral of \(f\left ( x,y\right ) \) over each of these lines (since all these lines are parallel to line \(L\), then all of them will have the same \(\theta \), but they will have diﬀerent \(p\) 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 \(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 speciﬁc angle \(\theta _{i}\) we obtain a projection vector \(\mathbf {g}_{i}\) as shown in the following diagram

The projection shown above is a discrete function. It is a function of \(p\) and \(\theta \). Hence we have \(g\left ( p,\theta \right ) \). But for the same angle \(\theta \), \(g\) is a function of \(p\), So some books write \(g_{\theta }\left ( p\right ) \).

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

Once the projection \(g_{\theta }\left ( p\right ) \) is obtained, then it is converted to the (discrete) Fourier domain using FFT. The discrete Fourier transform is\[ G_{k}={\sum \limits _{n=0}^{N-1}}g_{n}\exp \left ( -\frac {2\pi i}{N}kn\right ) \ \ \ k=0,1,\cdots ,N-1 \] Hence we obtain the vector \(\mathbf {G}\) which is the discrete Fourier transform of the projection \(\mathbf {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.

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 \(\mathbf {G}\) obtained, a
ﬁlter is selected. Consider the ram-Lak ﬁlter. This ﬁlter has a frequency response as
follows^{1}

So, assume the ﬁlter frequency spectrum is given by the vector \(\mathbf {H}\), (this is a complex vector, since it is the frequency spectrum of the ﬁlter). Hence the ﬁltered backprojection is given by\[ \mathbf {\tilde {G}}=\mathbf {GH}\] 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 \(\mathbf {\tilde {g}}\) by inverse discrete Fourier transform of \(\mathbf {\tilde {G}}\) and this is given by\[ \tilde {g}_{n}=\frac {1}{N}{\sum \limits _{k=0}^{N-1}}\tilde {G}_{k}\exp \left ( \frac {2\pi i}{N}kn\right ) \ \ \ \ \ n=0,1,\cdots ,N-1 \] Now we have obtain the spatial representation \(\mathbf {\tilde {g}}\) of the projection \(\mathbf {g}\) 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 \[ \mathbf {\tilde {g}=\operatorname {real}}\left ( \mathbf {\tilde {g}}\right ) \] 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 \(\tilde {g}_{\theta }\left ( p\right ) \), where the tilde indicates this is a ﬁltered projection. Now that we have calculated \(\tilde {g}_{\theta }\left ( p\right ) \), we can obtain the backprojection, which will be a 2D function. Assume the original function (which we do not know in practice, was \(f\left ( x,y\right ) \), then let the backprojection function be \(f_{B}\left ( x,y\right ) \). Hence \[ f_{B}\left ( x,y\right ) ={\int \limits _{\theta =0}^{\pi }}\tilde {g}_{\theta }\left ( p\right ) d\theta \] But \(p=x\cos \theta +y\sin \theta \), hence the above becomes\[ f_{B}\left ( x,y\right ) ={\int \limits _{\theta =0}^{\pi }}\tilde {g}_{\theta }\left ( x\cos \theta +y\sin \theta \right ) d\theta \] What the above is saying, is that to ﬁnd \(f_{B}\left ( x,y\right ) \) at some \(x,y\) position, the angle \(\theta \) is changed to cover all the angles from zero to \(180^{0}\), and for each angle in this interval, the function \(\tilde {g}_{\theta }\left ( p\right ) \) is evaluated at \(p=x\cos \theta +y\sin \theta \). The sum of all these gives \(f_{B}\) at that \(x,y\) values. We do these for each \(x,y\) value in the region to obtain all the values of \(f_{B}\). 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\[ f_{B}\left ( x,y\right ) =\left \{ {\sum \limits _{m=0}^{N-1}}\tilde {g}_{\theta _{m}}\left ( x\cos \theta _{m}+y\sin \theta _{m}\right ) \right \} \Delta \theta \] Since there are \(N\) angles, then we divide \(\pi \) by \(N\) to obtain each speciﬁc angle in the range. Hence \[ \theta _{m}=\frac {\pi }{N}m \] And also we see that \(\Delta \theta =\frac {\pi }{N}\) So the sum becomes\begin {equation} f_{B}\left ( x,y\right ) =\frac {\pi }{N}{\sum \limits _{m=0}^{N-1}}\tilde {g}_{\theta _{m}}\left ( x\cos \left ( \frac {\pi }{N}m\right ) +y\sin \left ( \frac {\pi }{N}m\right ) \right ) \tag {1} \end {equation}

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 \(\theta \) 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.

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\[ f_{B}\left ( x,y\right ) =\frac {\pi }{N}{\sum \limits _{m=0}^{N-1}}\tilde {g}_{\theta _{m}}\left ( x\cos \left ( \frac {\pi }{N}m\right ) +y\sin \left ( \frac {\pi }{N}m\right ) \right ) \] If we call the sum as \(img,\) then the above becomes\[ f_{B}\left ( x,y\right ) =\frac {\pi }{N}img \] Notice however, that matlab iradon.m does the following at the end: (type iradon.m to see)\begin {equation} f_{B}=\frac {\pi }{2N}img \tag {2} \end {equation} 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 \(360^{0}\), then there is a double counting involved (since the image can be fully constructed from \(0\) to \(180^{0}\), and hence it assumes the angles given have full range of \(2\pi \), 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\begin {align*} f_{B_{1}} & =\left ( \frac {\pi }{2\times 2}img_{1}\right ) /2=\frac {\pi }{8}img_{1}\\ f_{B_{2}} & =\left ( \frac {\pi }{2\times 2}img_{2}\right ) /2=\frac {\pi }{8}img_{2}\\ & \vdots \\ f_{B_{N}} & =\left ( \frac {\pi }{2\times 2}img_{N}\right ) /2=\frac {\pi }{8}img_{N} \end {align*}

Now we add all the above and obtain\begin {equation} f_{B_{1}}+f_{B_{2}}+\cdots +f_{B_{N}}=\frac {\pi }{8}\left ( img_{1}+img_{2}+\cdots +img_{N}\right ) \tag {3} \end {equation} Now, \(img_{1}+img_{2}+\cdots +img_{N}\) 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 \(360^{0}\) range and not \(0\) to \(180^{0}\), we will double add the intensity. Hence (3) becomes\begin {align} f_{B_{1}}+f_{B_{2}}+\cdots +f_{B_{N}} & =\frac {\pi }{8}\ \left ( 2img\right ) \nonumber \\ & =\frac {\pi }{4}img\tag {4} \end {align}

Compare (4), which is the one-at-time, with (2), which is the all-at-once, we see that\begin {align*} f_{B_{1}}+f_{B_{2}}+\cdots +f_{B_{N}} & =\frac {\pi }{4}img\\ f_{B} & =\frac {\pi }{2N}img \end {align*}

Therefore, we see that\begin {align*} \frac {f_{B_{1}}+f_{B_{2}}+\cdots +f_{B_{N}}}{f_{B}} & =\frac {\frac {\pi }{4}img}{\frac {\pi }{2N}img}\\ & =\frac {N}{2} \end {align*}

Hence the one-at-time result is \(\frac {N}{2}\) 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 \(\frac {N}{2}\) times the image that would have been obtained if using the all-at-once method.

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