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 \(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 difference.
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 specific \(p\) and specific \(\theta \) satisfy \(x\cos \theta +y\sin \theta =p\)
Assuming there exist a function \(f\left ( x,y\right ) \) defined 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 differential 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 specific \(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 different \(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 specific 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 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 \(\mathbf {G}\) obtained, a filter is selected. Consider the ram-Lak filter. This filter has a frequency response as follows1
So, assume the filter frequency spectrum is given by the vector \(\mathbf {H}\), (this is a complex vector, since it is the frequency spectrum of the filter). Hence the filtered backprojection is given by\[ \mathbf {\tilde {G}}=\mathbf {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 \(\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 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 \[ \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 filter ram-lak used to obtain the filtered backprojection
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 \(\tilde {g}_{\theta }\left ( p\right ) \), where the tilde indicates this is a filtered 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 find \(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 specific 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 filled radially. This diagram from Kak illustrates this:
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.
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\[ 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.