In the last post, many moons ago, I introduced the 2-D FFT and discussed the magnitude and phase components of the spatial Fourier domain. In the next few posts, I would like to describe a concrete application of the 2-D FFT, namely blurring. Blurring and deblurring are useful image processing techniques that can be used in a variety of fields. I will explain what blur is mathematically and how it is performed artificially. In future posts, I’ll go into more depth about what happens in the spatial domain, different types of blur, and some current deblurring technology.
When we look at a blurred image, we see that the colors from various objects are spread out in the image over an area around the objects. Mathematically, this is exactly what is happening. Instead of the light from one point in the scene corresponding with the light at one point in the image, the light from each point in the scene is spread out over multiple points in the image. The function that describes this spread is called the point spread function (PSF). The PSF is also the impulse response of the imaging system, and in applications dealing with blur, is sometimes referred to as the “blur kernel”.
Let’s use an example to demonstrate the FFT method of blurring. Here is a picture of Toni Morrison, an African-American author and winner of the Pulitzer and Nobel prizes. We will use a Gaussian blur kernel, a common type of blur, to blur this image.
The following code shows how to load the image and plot it. The imread function is used to load the jpeg file and is the same function used in the last post. However, in this post, we use imagesc to display the image. After loading the image into the Matlab workspace, we delete the second and third indices of the third dimension of the image, since they are nearly identical to the first index. The rest of the code after the plotting command just format the figure and axes so that the aspect ratio of the image is preserved and the image looks nice.
%Blur Demo %Import image origimage = imread('tonimorrison','jpg'); %Reduce image to 2-D origimage = origimage(:,:,1); %Plot image figure, imagesc(origimage) axis square colormap gray title('Original Image') set(gca, 'XTick', , 'YTick', )
The following code shows how to create a Gaussian blur kernel. First, the size of the blur kernel is selected. The size of the kernel determines the amount of blur. Second, the variance of the Gaussian function is chosen. A larger variance will result in more blur, and a smaller variance will result less blur, as the color from each point is spread out over a smaller area. Third, the mean of the Gaussian is chosen. The mean should just be in the center of the kernel image. Next, the values of the blur kernel matrix are filled using the meshgrid function and the equation for a two-dimensional Gaussian with mean m and variance s^2. Finally, the blur kernel is plotted the same way as the image.
%Blur Kernel ksize = 31; kernel = zeros(ksize); %Gaussian Blur s = 3; m = ksize/2; [X, Y] = meshgrid(1:ksize); kernel = (1/(2*pi*s^2))*exp(-((X-m).^2 + (Y-m).^2)/(2*s^2)); %Display Kernel figure, imagesc(kernel) axis square title('Blur Kernel') colormap gray
Now that we have the image and blur kernel, we can blur the image. But first, let’s take a step back and talk about the theory behind the code. As we said before, blurring is just spreading out the information from each point into the surrounding points. The function that accomplishes this in the image domain is called convolution. Convolution involves sliding a window over an image, then setting each point at the center of the window to the sum of the pointwise multiplications of the values of the window and the values of the points in the image that overlap each other. This website gives a good visual and mathematical description of how this works. Although convolution will correctly blur an image, there exists another method that is faster, called the Fast Fourier Transform (FFT). According to the convolution theorem, convolution in the time (or image) domain is equivalent to multiplication in the frequency (or spatial) domain. The image domain is the 2-D equivalent of the time domain, and the spatial domain is the 2-D equivalent of the Fourier domain. Thus, we have:
Let’s look through the code, now. In order for the FFT of the kernel to be multiplied with the FFT of the image, the kernel must be placed into an image of the same size as the original image. The kernel can be placed anywhere in the image; it doesn’t really matter. After the kernel is placed into a new image, 2-D FFTs of the kernel image and original image are taken. The FFT of the kernel is then modified so that there are no zero values. This is not important now because we are performing a multiplication with the FFT of the kernel, but later, when we deblur, we will need to perform a division, and any zero values in the FFT of the kernel will prompt a warning. To ensure that we are performing symmetric operations, we should remove the zero values from the FFT of the kernel first. After the zeros are removed from the FFT of the kernel, a multiplication is performed with the FFT of the original image. This is a pointwise multiplication, not a matrix multiplication, which would give much different results and probably an error, in this case. Next, the FFT of this product is transformed back into the image domain with an inverse 2-D FFT. This procedure follows the second equation mentioned earlier. Finally, the blurred image is displayed using imagesc. Scroll down below the code to see the results.
%Embed kernel in image that is size of original image [h, w] = size(origimage); kernelimage = zeros(h,w); kernelimage(1:ksize, 1:ksize) = kernel; %Perform 2D FFTs fftimage = fft2(double(origimage)); fftkernel = fft2(kernelimage); %Set all zero values to minimum value fftkernel(find(fftkernel == 0)) = 1e-6; %Multiply FFTs fftblurimage = fftimage.*fftkernel; %Perform Inverse 2D FFT blurimage = ifft2(fftblurimage); %Display Blurred Image figure, imagesc(blurimage) axis square title('Blurred Image') colormap gray set(gca, 'XTick', , 'YTick', )
Ok, so this looks really bad. Although the image is blurred, it also seems to be aliased in both dimensions. So, where did we go wrong? The answer is that we need to use padding before we do the FFT. There are several good ways to pad the edges of the image. The method I use below sets each point in the padded region around the image to the closest edge pixel. To visualize the type of padding I use, imagine vertical lines in the padded regions above and below the image, horizontal lines to the right and left of the image, and squares equal to the corner pixels in the corners of the padded region. This method is good enough for our purposes, but other methods could also be used, such as setting the padded areas to mirror images of the regions of the image close to the edges. The essential requirement of the padding is that it must have approximately the same values as the regions close to the edges of the image.
To implement the padding function, create the padimage function, as written below, in a separate file or within the main file, and modify the previous block of code in the following way.
%Pad image origimagepad = padimage(origimage, ksize); %Embed kernel in image that is size of original image + padding [h1, w1] = size(origimagepad); kernelimagepad = zeros(h1,w1); kernelimagepad(1:ksize, 1:ksize) = kernel; %Perform 2D FFTs fftimagepad = fft2(origimagepad); fftkernelpad = fft2(kernelimagepad); fftkernelpad(find(fftkernelpad == 0)) = 1e-6; %Multiply FFTs fftblurimagepad = fftimagepad.*fftkernelpad; %Perform Reverse 2D FFT blurimagepad = ifft2(fftblurimagepad); %Remove Padding blurimageunpad = blurimagepad(ksize+1:ksize+h,ksize+1:ksize+w); %Display Blurred Image figure, imagesc(blurimageunpad) axis square title('Blurred Image - with Padding') colormap gray set(gca, 'XTick', , 'YTick', )
function Ipad = padimage(I, p) %This function pads the edges of an image to minimize edge effects %during convolutions and Fourier transforms. %Inputs %I - image to pad %p - size of padding around image %Output %Ipad - padded image %Find size of image [h, w] = size(I); %Pad edges Ipad = zeros(h+2*p, w+2*p); %Middle Ipad(p+1:p+h, p+1:p+w) = I; %Top and Bottom Ipad(1:p, p+1:p+w) = repmat(I(1,1:end), p, 1); Ipad(p+h+1:end, p+1:p+w) = repmat(I(end,1:end), p, 1); %Left and Right Ipad(p+1:p+h, 1:p) = repmat(I(1:end,1), 1, p); Ipad(p+1:p+h, p+w+1:end) = repmat(I(1:end,end), 1, p); %Corners Ipad(1:p, 1:p) = I(1,1); %Top-left Ipad(1:p, p+w+1:end) = I(1,end); %Top-right Ipad(p+h+1:end, 1:p) = I(end,1); %Bottom-left Ipad(p+h+1:end,p+w+1:end) = I(end,end); %Bottom-right
As you can see, the image is now blurred with no strange edge effects. In future posts, I’ll discuss the effects of blurring in the spatial domain, different types of blur, and some current deblurring technology.
*In the original post, there were several errors in the section of the code that performed blurring with padding, which have been corrected. Additionally, there was a mismatch between the code that created the final deblurred figure and the code that was presented in the post, which caused some confusion for the readers. The deblurred figure in the post was actually created using the blurred figure before its padding was removed, which made a large difference in the final product. Deblurring is actually more complicated that what I presented in the original post, and I will address deblurring methods in later posts in this series. I apologize for any confusion. In addition, I would like to thank several people who pointed out these errors in the comments section.