Description
The goal of this lab is to help you to • get an understanding of frequency analysis in image data by using the Fourier transform, expressed spatially or in the Fourier domain, • become familiar with the two-dimensional Fourier transform, by studying and testing its properties in practice,
• learn the relation between the continuous and discrete Fourier transform, and design a smoothing filter for Gaussian convolution,
• understand the effect of different smoothing operations on various kinds of noise, • understand the effect of smoothing on different image resolutions generated by subsampling, • get practical experience of differences between synthetic and real data.
The purpose of exercise 1 is for you to become familiar with the characteristics of the Fourier transform. The purpose of exercise 2 is for you to design a Gaussian smoothing filter using the Fourier transform and to understand the relation between continuous and discrete cases.
The purpose of exercise 3 is to study different smoothing filters and their effects on noise and subsampling. Remaining exercises are more problem-oriented and experimental in nature and formulations are thus more concise. If you happen to have difficulties with Matlab or the environment in general contact the lab assistants during the scheduled lab hours.
Reporting: For you to efficiently report your labs, you need to create a script file that reproduces the experimental results. You should also summarize results and conclusions from your experiments, writing down answers to explicitly stated questions in the answer sheet. For many exercises it is recommended that you create illustrations with multiple images simultaneously shown on screen. For this purpose use the embedded Matlab command subplot (see “help subplot”).
As prerequisites to this lab exercise, you should have studied the course material regarding spatial filtering and image restoration. To get started with the lab you need to set the path, in order to access the image samples and functions of the course library.
1 Properties of the discrete Fourier transform
1.1 The continuous and discrete Fourier transform From the definitions used during the lectures the continuous Fourier transform ˆf : R 2 → C of a twodimensional signal f : R 2 → R is ˆf(ω) = FC(f)(ω) = Z x∈R2 f(x) e −iωT x dx (1) and the Fourier inversion theorem states that f(x) = F −1 C ( ˆf)(x) = 1 (2π) 2 Z ω∈R2 ˆf(ω) e +iωT x dω.
(2) Correspondingly, in the discrete case the Fourier transform Fˆ : [0..N − 1]2 → C of a quadratic image F : [0..N − 1]2 → R of N2 image elements is defined as Fˆ(u) = FD(F)(u) = 1 N X x∈[0..N−1]2 F(x) e −2πiuT x N (3) with the corresponding inversion theorem F(x) = F −1 D (Fˆ)(x) = 1 N X u∈[0..N−1]2 Fˆ(u) e +2πiuT x N .
(4) There expressions are symmetric with respect to the factor 1/N in the Fourier transform and inverse respectively. Note that in Matlab the fast Fourier transform (FFT) is implemented using a factor 1 in the FFT-routine, and factor 1/N2 in the inverse.
1.2 Definition domains In the continuous case the spatial variable x = (x1, x2) and the frequency variable ω = (ω1, ω2) are defined in the whole two-dimensional plane, i.e. x, ω ∈ R 2 , while the corresponding discrete variables x = (x1, x2) and u = (u1, u2) are defined in the interval [0..N − 1]2 , i.e. x1, x2, u1, u2 ∈ [0..N − 1] = {0, 1, . . . , N − 1}.
From the periodicity of the basis functions in the discrete case, thus follows that the discrete Fourier transform will become periodic with a period of N. To maintain the consistency with the inversion theorem we similarly have to consider the original signal as periodic with the same period. By comparing equations (1) and (3) we can relate the continuous angular frequency variable ω to the discrete frequency variable u: ωD = 2π u N . (5) Since u ∈ [0..N − 1]2 we may from this relation treat the frequency variable ωD as being defined in the interval [0, 2π] 2 .
Furthermore, the discrete Fourier transform may similarly be treated as periodic with a period of 2π, and we may change the definition domain to the interval [−π, π] 2 without losing in generality. This corresponds to a centering operation in the discrete Fourier transform.
3 1.3 Basis functions The Fourier transform can be regarded as a change of basis functions, from compact (discrete) delta functions defined on the Cartesian image domain to complex exponential function with maximum spatial extension. In the discrete case this is even more obvious noting that the discrete Fourier transform can be found by pre- and post-multiplying the image with orthogonal matrices.
Doing so we may consider the complex pixel values in the Fourier transform Fˆ, as components of the discrete image F, with respect to the new basis. Specifically, the basis element corresponding to an image point at coordinates (p, q) in Fˆ is proportional to e i2π p·0 N . . . e i2π p·(N−1) N e i2π q·0 N · · · e i2π q·(N−1) N
These basis vectors are the discrete counterparts to the complex exponential functions ew(x) = e iωT x You may convince yourself that this is true by calculating the inverse discrete Fourier transform of an image Fˆ, that is zero everywhere except for a single point (p, q) at which the value is one. For such an image the expansion of Fˆ only contains one term, i.e. the basis vector given by the index (p, q). To visualize this, define a 128 × 128 pixel image, with elements defined as Fhat = zeros(128, 128); Fhat(p , q ) = 1; with (p, q) = (5, 9).
Display this image on screen using showgrey. This function is presented in Lab 0B: Elementary Image Operations. Type help showgrey in the command prompt of Matlab for a short description of the function and its parameters.
Thereafter, compute its inverse discrete Fourier transform with F = ifft2(Fhat); and look at its real and imaginary parts, as well as its magnitude and phase, by writing Fabsmax = max(abs(F(:))); showgrey(real(F), 64, -Fabsmax, Fabsmax) showgrey(imag(F), 64, -Fabsmax, Fabsmax) showgrey(abs(F), 64, -Fabsmax, Fabsmax) showgrey(angle(F), 64, -pi, pi) A clever way of organizing the experiment is by writing a small procedure fftwave, that assembles the resulting images into a single figure in Matlab, with the help of the command subplot and illustrative headings.
1 Please use this procedure as of prototype example of how the experimentation in later exercises can be arranged. 1A template of such a procedure can be found in fftwavetemplate.m in the course library.
function fftwave(u, v, sz) if (nargin < 2) error(’Requires at least two input arguments.’) end if (nargin == 2) sz = 128; end Fhat = zeros(sz); Fhat(u, v) = 1; F = ifft2(Fhat); Fabsmax = max(abs(F(:))); subplot(3, 2, 1); showgrey(Fhat); title(sprintf(’Fhat: (u, v) = (%d, %d)’, u, v)) % What is done by these instructions? if (u <= sz/2) uc = u – 1; else uc = u – 1 – sz; end if (v <= sz/2) vc = v – 1; else vc = v – 1 – sz; end wavelength = 0.0; % Replace by correct expression amplitude = 0.0; % Replace by correct expression subplot(3, 2, 2);
showgrey(fftshift(Fhat)); title(sprintf(’centered Fhat: (uc, vc) = (%d, %d)’, uc, vc)) subplot(3, 2, 3); showgrey(real(F), 64, -Fabsmax, Fabsmax); title(’real(F)’) subplot(3, 2, 4); showgrey(imag(F), 64, -Fabsmax, Fabsmax); title(’imag(F)’) subplot(3, 2, 5); showgrey(abs(F), 64, -Fabsmax, Fabsmax); title(sprintf(’abs(F) (amplitude %f)’, amplitude)) subplot(3, 2, 6); showgrey(angle(F), 64, -pi, pi); title(sprintf(’angle(F) (wavelength %f)’, wavelength))
5 Question 1: Repeat this exercise with the coordinates p and q set to (5, 9), (9, 5), (17, 9), (17, 121), (5, 1) and (125, 1) respectively. What do you observe?
Question 2: Explain how a position (p, q) in the Fourier domain will be projected as a sine wave in the spatial domain. Illustrate with a Matlab figure.
Question 3: How large is the amplitude? Write down the expression derived from Equation (4) in these notes. Complement the code (variable amplitude) accordingly.
Question 4: How does the direction and length of the sine wave depend on p and q? Write down the explicit expression that can be found in the lecture notes. Complement the code (variable wavelength) accordingly.
Question 5: What happens when we pass the point in the center and either p or q exceeds half the image size? Explain and illustrate graphically with Matlab!
Question 6: What is the purpose of the instructions following the question What is done by these instructions? in the code? 1.4 Linearity Define some rectangular shaped 128 × 128 pixel test images by F = [ zeros(56, 128); ones(16, 128); zeros(56, 128)]; G = F’; H = F + 2 * G; and display them with showgrey.
Then compute the discrete Fourier transform of the images by writing Fhat = fft2(F); Ghat = fft2(G); Hhat = fft2(H); and show their Fourier spectra with showgrey(log(1 + abs(Fhat))); showgrey(log(1 + abs(Ghat))); showgrey(log(1 + abs(Hhat))); Also try to run the following showgrey(log(1 + abs(fftshift(Hhat)))); and explain why the fftshift and the log commands are used here. Another way of achieving the same effect is using the Matlab function showfs located in the course library. To better understand the role of the logarithm function, see also Lab 0B ”Elementary image operations”.
Question 7: Why are these Fourier spectra concentrated to the borders of the images? Can you give a mathematical interpretation? Hint: think of the frequencies in the source image and consider the resulting image as a Fourier transform applied to a 2D function. It might be easier to analyze each dimension separately!
Question 8: Why is the logarithm function applied?
Question 9: What conclusions can be drawn regarding linearity? From your observations can you derive a mathematical expression in the general case? 1.5 Multiplication With F and G as previously defined F = [ zeros(56, 128); ones(16, 128); zeros(56, 128)]; G = F’; Try the following commands showgrey(F .* G); showfs(fft2(F .* G)); and explain the results. (The notation F .* G in Matlab refers to point-wise multiplication of corresponding matrix elements.)
Question 10: Are there any other ways to compute the last image? Remember what multiplication in Fourier domain equals to in the spatial domain! Perform these alternative computations in practice.
1.6 Scaling Define a test image F = [zeros(60, 128); ones(8, 128); zeros(60, 128)] .* … [zeros(128, 48) ones(128, 32) zeros(128, 48)]; and display it with showgrey. Determine the discrete Fourier transform and look at the magnitude with showfs.
Question 11: What conclusions can be drawn from comparing the results with those in the previous exercise? See how the source images have changed and analyze the effects of scaling. 1.7 Rotation Activate a new figure with figure, but keep the Fourier spectrum of F in the first figure. Rotate F by for example alpha = 30◦ and show the results with G = rot(F, alpha ); showgrey(G) axis on Then calculate the discrete Fourier transform of the rotated image with Ghat = fft2(G); and display the results with showfs.
Do you recognize it? Finally, rotate the spectrum back by the same angle with Lab 1: Filtering operations 7 Hhat = rot(fftshift(Ghat), -alpha ); and show the results by writing showgrey(log(1 + abs(Hhat))) Assemble the original images and their Fourier spectra into the same figure with the subplot command. Vary the angle with different values of alpha, for example alpha = 30◦ , 45◦ , 60◦ , 90◦ and illustrate with at least one case.
Question 12: What can be said about possible similarities and differences? Hint: think of the frequencies and how they are affected by the rotation.
1.8 Information in Fourier phase and magnitude Above we have primarily used the magnitude of the Fourier transform when we have visualized the transform as an image. This visualization technique is also what dominates in literature.
As an illustration of the limitations of this way of visualizing, we will in this section devote ourselves to a bit of image manipulation, where we simply replace the power spectrum for a given image f with a power spectrum of the form 2 | ˆf(ω)| 2 = 1 a + |ω| 2 (6) In the course library there is a function pow2image that performs this (nonlinear) function. (Take a look at the code by typing “type pow2image”.) Apply this function on for example the following images phonecalc128, few128, nallo128 and study the results on screen (for very small values of a ≈ 10−10).
An image can be loaded like this: img = phonecalc128. Conclusions? As a comparison you should apply the function randphaseimage that keeps the magnitude of the Fourier transform, but replaces the phase information with a random distribution.
Question 13: What information is contained in the phase and in the magnitude of the Fourier transform? 2 Gaussian convolution implemented via FFT Here we study how a smoothing filter can be designed using the Fourier transform, in order to perform a 2D Gaussian convolution on images.
2.1 Continuous Gaussian convolution Gaussian convolution means that we convolve a given image fin with a Gaussian kernel. In the continuous case the output image fut is given by fut(x, y) = Z ∞ ξ=−∞ Z ∞ η=−∞ fin(x − ξ, y − η) g(ξ, η; t) dξ dη (7) where g(x, y; t) = 1 2πt e −(x 2+y 2 )/(2t) . (8) 2For continuous signals (and in the case a = 0) this form of Fourier spectrum can be derived as an idealized model of images that contain different types of image structures distributed on all scales.
(The purpose of parameter a is simply to avoid division by zero.)
Correspondingly, in the Fourier domain ˆfut(ω1, ω2) = ˆg(ω1, ω2; t) ˆfin(ω1, ω2) (9) where ˆfut and ˆfin are the Fourier transforms of fut and fin respectively, and ˆg(·, ·; t) denotes the Fourier transform of the Gaussian function, that is gˆ(ω1, ω2; t) = Z ∞ ω1=−∞ Z ∞ ω2=−∞ g(x, y; t) e −i(ω1x+ω2y) dx dy = e −(ω 2 1+ω 2 2 ) t/2 . (10)
2.2 Discretization in the spatial domain To discretize the Gaussian convolution, we may either choose to discretize the convolution operation (7) in the spatial domain, or proceed by discretizing the Fourier transform of the Gaussian function (10). In this exercise we will try the first method that is the spatial domain. If we discretize the integral in (7) with the trapezoidal rule and set the step length to one (corresponding to a unit distance between neighboring image elements), we get fut(x, y) = X∞ m=−∞ X∞ n=−∞ fin(x − m, y − m) g(m, n; t) (11) which corresponds to a discrete convolution with a sampled version of the Gaussian function.
To perform this operation in practice, we may thus use either one of the following methods: With spatial discretization and spatial convolution: 1. Generate a filter based on a sampled version of the Gaussian function. 2. Convolve the image with this filter using the embedded Matlab-function conv2. This method is probably the most straight-forward one. Such explicit convolution in the spatial domain generally works very well, if one exploits the fact that the Gaussian function is separable.
This method is particularly suitable when the variance of the Gaussian kernel is small, and the kernel is truncated to a small filter of compact support. However, since the computational work grows linearly with the number of non-zero elements in the filter, this method can be computationally expensive for large values of t, in particular if the separability is not exploited.
In practice, this is not the case with the Matlab-function conv2. With spatial discretization and convolution via FFT: 1. Generate a filter based on a sampled version of the Gaussian function. 2. Fourier transform the original image and the Gaussian filter. 3. Multiply the Fourier transforms. 4. Invert the resulting Fourier transform. This method has the characteristic that it is essentially independent of the size of the Gaussian kernel. In practice, you always generate a Gaussian kernel of the same size as the original image to be filtered; we assume here that this is a power of two and large enough for the shape of the Gaussian kernel to be kept within the given image size.
Specifically, with the image processing system we use, the second method is facilitated by an efficient implementation of FFT available in Matlab. For the above mentioned reasons the last method is in this case the preferable one, even if it involves more steps.
2.3 Filtering procedure Write a function that performs Gaussian filtering using FFT! Write one Matlab-procedure gaussfft(pic, t) that, using the fast Fourier transform, convolves the image pic with a two-dimensional Gaussian function of arbitrary variance t via a discretization of the Gaussian function in the spatial domain, as described in steps 1-4 above.
In Matlab you may simply generate a coordinate system with the command meshgrid and move the position of the origin with the command fftshift. Proposed test procedure: When you perform this operation you might want to make sure that your convolution procedure has an approximatively correct behavior by analyzing its impulse response.
Hence, generate the impulse response explicitly psf = gaussfft(deltafcn(128, 128), t); for different values of t (e.g. t = 0.1, 0.3, 1.0, 10.0 and 100.0). Look at the result and compute also the (spatial) covariance matrix of the Gaussian function with variance(psf) Compare your results to the ideal continuous case, for which the covariance matrix is t multiplied by the identity matrix. C(g(·, ·; t)) = t 1 0 0 1 (12) If you find it difficult to get this to properly work, a routine discgaussfft, that convolves an image with the discrete version of the Gaussian kernel, is provided in the course library. Your procedure should be different and your results are not expected to be exactly the same! This function might help you to understand the effect of Gaussian filtering by visualizing the results.
Question 14: Show the impulse response and variance for the above mentioned t-values. What are the variances of your discretized Gaussian kernel for t = 0.1, 0.3, 1.0, 10.0 and 100.0?
Question 15: Are the results different from or similar to the estimated variance? How does the result correspond to the ideal continuous case? Lead: think of the relation between spatial and Fourier domains for different values of t.
Question 16: Convolve a couple of images with Gaussian functions of different variances (like t = 1.0, 4.0, 16.0, 64.0 and 256.0) and present your results. What effects can you observe?
3 Smoothing
3.1 Smoothing of noisy data Load the image office256 from the image database of the course office = office256; and create two noisy images using the Matlab-functions gaussnoise and sapnoise. (These functions exist in the course function library): add = gaussnoise(office, 16); sap = sapnoise(office, 0.1, 255); Show the images on screen with showgrey.
What kind of noise do they represent? (Use the help command, if you are unable to guess.) Try to reduce the noise in the images add and sap with • Gaussian smoothing (using function gaussfft that you wrote in previous exercise, alternatively the function discgaussfft available in the course library) • Median filtering (using the function medfilt in the course library) • Ideal low-pass filtering (using the function ideal in the course library) Try suitable values of the parameters of each filter respectively (standard deviation for the Gaussian filter, window size for the median filter, cut-off frequency for the ideal low-pass filter).
Question 17: What are the positive and negative effects for each type of filter? Describe what you observe and name the effects that you recognize. How do the results depend on the filter parameters? Illustrate with Matlab figure(s).
Question 18: What conclusions can you draw from comparing the results of the respective methods?
3.2 Smoothing and subsampling We will now study the effect of smoothing when applied to images sampled at lower resolutions. There are many applications to subsampling motivated by performance or storage. Generate a series of subsampled versions of phonecalc256, keeping 1 pixel out of 2 (using the rawsubsample function in the course library).
Alternatively, smooth the original image before subsampling, using the Gaussian filter (gaussfft) or the ideal low-pass filter (ideal). To display the results you can use this template: img = phonecalc256; smoothimg = img; N=5; for i=1:N if i>1 % generate subsampled versions img = rawsubsample(img); smoothimg = % (smoothimg, ); smoothimg = rawsubsample(smoothimg); end subplot(2, N, i) showgrey(img) subplot(2, N, i+N) showgrey(smoothimg) end As the resolution of the image decreases, observe the changes in the original image and its smoothed variants.
Try the Gaussian and ideal low-pass filters and find suitable values for their parameters, that preserve most of the characteristics/quality with respect to the original image (iteration i = 1).
Question 19: What effects do you observe when subsampling the original image and the smoothed variants? Illustrate both filters with the best results found for iteration i = 4.
Question 20: What conclusions can you draw regarding the effects of smoothing when combined with subsampling? Hint: think in terms of frequencies and side effects.