Due date: 30 September 2024

This assignment deals with frequency-domain filtering, image restoration, and image reconstruction from projections. For simplicity, you may use single-component (intensity) images for all problems.

Blur due to an out-of-focus camera is better modelled as a disk rather than a Gaussian. In this problem, we will approximate a disk of radius \(r\) as the following discrete kernel:

\(h(x,y) = \begin{cases} 1 & \text{if } \sqrt{x^2+y^2}\in[0,r-0.5]\\ r+0.5-\sqrt{x^2+y^2} & \text{if } \sqrt{x^2+y^2}\in[r-0.5,r+0.5]\\ 0 & \text{if } \sqrt{x^2+y^2}\in[r+0.5,\infty] \end{cases}\)

For the test image, use a rectangular image whose width and height differ by a ratio of about 2.

Perform blurring using convolution in the spatial domain with mirror padding. You may use a built-in convolution function, or your own implementation. Either way, take care that the kernel is treated as properly centred, with \(x,y=-a,\dots,a\) not \(0,\dots,2a+1\).

Implement the same blurring via multiplication in the frequency domain, including mirror padding. You may use a built-in FFT implementation. Make sure you get nearly the same result as in part (a), up to floating-point rounding error. Also show the intermediate images, as in the last slide of lecture 9; images in the frequency domain should be visualized using the log of the spectrum.

Explain briefly why the (shifted) DFT of \(h\) appears elliptical, even though \(h\) itself has circular symmetry. (Hint: This is related to the role of \(M\) and \(N\) in the discrete Fourier transform.)

Measure the compute time \(T(r)\) of the implementations in (a) and (b) with increasing radii \(r = 1, 2, 4, 8, 16, \dots\) Keep increasing the size until \(T(r)\) > 1 minute or \(r > \min(M, N)/2\). Plot the timings of both methods on the same log-log plot.

Obtain an image \(f\) with some periodic high-frequency content, for example: a close-up of your monitor screen, a photo through a screen door, or a high-resolution scan of a photo in a newspaper. Make sure the periodicity is visible as spikes in the Fourier spectrum, otherwise try a different image.

Rescale the image to a much smaller size by using basic resampling and bilinear interpolation, as in Assignment 1. Choose the scaling factor so that there are visible aliasing artifacts.

Explain in your report how to choose the cutoff frequency \(D_0\) of a low-pass filter to eliminate aliasing, as a function of the scaling factor. (You may get different \(D_0\)’s in the horizontal and vertical directions.) Apply your choice of low-pass filter in the frequency domain, and then perform rescaling as before. Show both the low-pass filtered image in the spatial domain and the rescaled image.

Compute the Fourier spectrum of the original image without padding, and identify the spikes corresponding to the undesired periodic image content. Construct a notch-reject filter to remove the periodic component as best as you can. Describe your procedure in the report, showing all the relevant intermediate images: the original Fourier spectrum \(|F|\), the filter function \(H\) (without log), the filtered spectrum \(|G|\), and the final image \(g\).

To implement an

*optimum*notch filter, we will need the means, variances, and covariances of \(f\) and \(\eta\) in an \(m\times n\) neighbourhood of each pixel. Naively, each of these will take \(O(MNmn)\) total time. Explain how they can be computed in \(O(MN(m+n))\) time instead, using the fact that \({\rm var}(X)\) = \({\rm mean}(X^2)\) − \({\rm mean}^2(X)\) and \({\rm cov}(X,Y)\) = \({\rm mean}(XY)\) − \({\rm mean}(X){\rm mean}(Y)\) for any quantities \(X,Y\). Implement your method and demonstrate it by computing an image \(s(x,y)={\rm var}_{x,y}(f)\) with some moderately large window.Use the same notch filter in (c) and your implementation of (d) to perform optimum notch filtering. Visualize the estimated noise \(\eta\), the weight function \(w\), the modulated noise \(w\eta\), and the restored image \(g = f-w\eta\).

In this problem on denoising, choose a test image \(f\) that has some regions with smoothly varying intensities, and some regions with lots of detail / texture.

Implement operations that add Gaussian noise of a specified standard deviation \(\sigma\), or salt-and-pepper noise with specified probabilities \(p_s,p_p\), to a given image. Create two noisy images, \(g_1\) which has Gaussian noise with \(\sigma=0.4I_\max\), and \(g_2\) which has salt-and-pepper noise with \(p_s=p_p=0.2\).

Implement an arithmetic mean filter and a median filter with user-specified width \(w=2a+1\). Obtain restored images \(\hat f\) using a mean filter on \(g_1\) with varying sizes \(w=3,5,7,\dots\); plot the mean-squared error (MSE) \(\tfrac1{MN}\|\hat f-f\|^2\) as a function of \(w\) and show the image with the lowest MSE. Repeat using a median filter with the same sizes, add its MSE to the same plot, and show the best resulting image. Finally, repeat with both mean and median filters on \(g_2\).

You should see the MSE go down and then up again as \(w\) increases. To see why, perform the following additional experiment: obtain the Gaussian noise image \(\eta_1 = g_1 - f\), and apply the mean filter (let’s call it \(h\)) to \(f\) and \(\eta_1\) separately. Plot the MSE due to blurring \(\frac1{MN}\|h*f-f\|^2\), the MSE due to noise \(\frac1{MN}\|h*\eta_1\|^2\), and the total MSE \(\frac1{MN}\|h*g_1-f\|^2\) on the same plot as a function of \(w\). Does this explain the observed behaviour?

Now let us perform deconvolution. Choose the test image from problem 1 or problem 3, apply a disk blur from problem 1 using a moderately-sized radius, and add a small amount of Gaussian noise (say \(\sigma\) between \(0.05I_\max\) and \(0.1I_\max\)) to get the degraded image \(g\).

Perform inverse filtering and verify that it does not work well. To avoid division by zero, replace \(H(u,v)\) by some small \(\epsilon\sim10^{-6}\) if \(|H(u,v)|<\epsilon\).

In your report, derive the formula for \(S_\eta(u,v)=E\{|N(u,v)|^2\}\) for the case where \(\eta\) is zero-mean Gaussian noise of standard deviation \(\sigma\), using the fact that the mean and variance are additive for independent random variables. Verify your formula numerically by taking the DFT of \(\eta\) and checking its values.

To perform Wiener filtering, we know \(S_\eta\) from part (b), but we still need \(S_f(u,v)=E\{|F(u,v)|^2\}\). Implement this in two ways. First, use the actual value of \(|F(u,v)|^2\) from the known original image (this is of course not possible outside of synthetic experiments). Second, take \(S_\eta/S_f\) to be some constant \(k\) independent of \((u,v)\). Try a reasonable range of values, plot the MSE as a function of \(k\), and show the plot as well as the restored image corresponding to the best \(k\).

Let us see what happens if we use the wrong guess of the degradation transfer function. Apply the Wiener filter with \(H\) corresponding to a disk radius that is too small, for example \(0.7r\), and experimentally choose a \(k\) that gives the best results. Similarly, choose a large disk radius e.g. \(1.5r\) and repeat the experiment. Discuss what you observe and informally explain why.

Submit a zip file that contains (i) all your code for the assignment, and (ii) a PDF report that includes your results for each of the assignment problems. In the report, each output image must be accompanied with a brief description of the procedure used to create it, and the values of any relevant parameters.

All images should ideally be saved in PNG format, which is lossless and so does not cause any information loss (other than quantization if the intensities are not already in 8-bit format). JPEG is permitted if your submission PDF is becoming too big for the upload limit.

Your assignment should be submitted on Moodle before midnight on the due date. Only one person in a group needs to submit. Late days are counted with a quantization of 0.5 days: if you cannot finish the assignment by midnight, get some sleep and submit by noon the following day.

Separately, each of you will individually submit a short response in a separate form regarding how much you and your partner contributed to the work in this assignment.