Topics:
Today's homework is subdivided into three parts (and one additional bonus part) and its main topic is the inverse of the forward (Radon) projection transform, that was the topic of last lab. The inverse transform starts with a sinogram, i.e. the projection image $J(\theta, p)$ and reconstructs the input image $f(x, y)$ in Cartesian space.
Naive approach The naive approach (the unfiltered back-projection) is just the reverse process of the projection – we smear the values of each projection back to the reconstructed image and add them together. We can write it as a sum over projection images $f_{\theta}(x, y)$ which we get from the measured projection $J(\theta, p)$ via $$f_{\theta}(x, y) = J(\theta, x \cos \theta + y \sin \theta)$$ The reconstructed image is then $$f_{\text{unfilt}}(x, y) = \mathcal{B}[\mathcal{R}f] = \sum_\theta {f_{\theta}(x, y)}$$
Optimal reconstruction From theoretical point of view, the central slice theorem gives us the optimal solution. The 2D-Fourier Transform $\mathcal{F}_2$ of the input image $f(x, y)$ is the same as the set of 1D-Fourier transformed projections $\mathcal{F}_1 [J(\theta, p)]$. However, since we have only very sparse samples in $\theta, p$, direct application of the inverse transform $\mathcal{F}_2^{-1}$ is impractical.
Nonetheless, we can use the central slice theorem to rewrite the Fourier equality $f(x, y) = \mathcal{F}_2^{-1} [\mathcal{F}_2 f]$ and with help of the Radon transform – $\mathcal{R}f$ to finally arrive at the formula:\begin{equation} f(x, y) = \mathcal{B} \{ \mathcal{F}_1^{-1} [|\omega| \cdot \mathcal{F}_1[\mathcal{R} f] \} \end{equation} i.e. the reconstructed image is obtained by back-projection $\mathcal{B}$ of the measured projections filtered by the ramp filter $|\omega|$ in Fourier domain.
Compared to the naive approach, the filtered back-projection computes $f_{\theta}$ from the filtered projection signal $J^{*}(\theta, p)$. To get $J^{*}(\theta, p)$ from the input projection data $J(\theta, p)$ we take a look at the inner part of equation (1) and write it only for a single projection:
$$\mathcal{B} \{ \mathcal{F}_1^{-1} [|\omega| \cdot \mathcal{F}_1[J(\theta, p)] \}$$
So we get either by multiplication in the Fourier domain or convolution in the spatial domain: $$J^{*}(\theta, p) = \mathcal{F}_1^{-1} [|\omega| \cdot \mathcal{F}_1[J(\theta, p)] = \mathcal{F}_1^{-1} [|\omega|] * J(\theta, p)$$
The ramp filter $|\omega|$ leads to optimal reconstruction, but when applied in Fourier-domain, it also serves as a high-pass filter (since unbounded on both ends) and is therefore magnifying noise frequencies. To minimize this effect, the filter is applied in combination with windows that reduce its frequency span. We will work with Hamming window (hamming
) and the Shepp-Logan, also sinc window (shepp-logan
). The last filter is ram-lak
which does no windowing, but only cuts of high frequencies. Use the provided function H = designFilter(f_type, N, f_d)
to get the frequency-domain filter - plot both $H$ and $\mathcal{F}_1^{-1}[H]$ to get better comprehension of the filtering functions.
Now we can summarize the (filtered) back-projection algorithm:
H = designFilter(f_type, N, f_d)
)
Get the archive with input sinogram (noisy_radon.mat
), reference image (phantom.mat
) and the designFilter.m
script.
myIradon( projim, thetas, f_type, f_d)where
projim
is the projection image of size $[m \times n_{\theta}]$ and theta
is vector of length $n_{\theta}$ containing the projection angles. The last parameter f_type
is optional and if set, it specifies the filtration to be applied during the back-projection, f_d
is a parameter of the cut-off frequency of the filter f_type
.
phantom
image for each of the 7 settings:
hamming
, ram-lak
and shepp-logan
and with filter cut-off frequencies – $d \in {0.7, 1.0}$
X, Y
as last time for positioning the output image $f(x, y)$
P
and[X, Y]
reflecting the spatial relation between $p$ and $(x, y)$, computation of the back-projected $f_{\theta}$ can be done with interp1 function).
designFilter
returns array of even length, if the radon image has odd length (in $p$), first enlarge $J(\theta_0, p)$ and then apply Fourier transform and the filter $H$.
phantom.mat
) in the reconstructed image
phantom.mat
) in the reconstructed image
iradon
function in Matlab (set parameter FILTER
- 'none'
for direct back projection, 'ram-lak'
or different for filtered back projection)