Fourier transformation can be a useful instrument in the analyses, processing and compression of data such as images or sound. There is a tremendous amount of information on two-dimensional discrete fast Fourier transformation available online. However, clearcut examples in R are scarce. It took my quite a while to piece together all information and figure out how to get it done in R. With this post I neither intend to explain the principles of fast Fourier transformation nor the mathematics behind it. However, I couldn't avoid talking about some technical stuff, without explaining it too much. Therefore, it is advisable to look around for some background info on the topic, before delving into this post. For starters you could have a look at Steven Lehar's page. Hopefully you can use the example presented here to your advantage.
For those readers I didn't scare off yet: let's get started. With a Fourier transform, information from discrete space (e.g. the brightnes of your image) is transformed to a frequency domain. This sounds very technical, which it is. If you took the time Googling a bit on the topic, you will understand that sinusoidal patterns in an image will show up as peaks in the frequency domain spectrum generated with a Fourier transform. One application of the Fourier transform is to detect and remove regular distortion patterns in images. In the example provided here, I will guide you through the following steps:
- Read an image and store as a matrix in the working memory
- Transform the image to the frequency domain
- Do some filtering magic
- Back transform the image to normal space and save the filtered image
Let's start with the image below. I've taken this picture years ago in the wonderful Dutch amusement park De Efteling. I've added some distortion, the kind that you will typically get when scanning a picture from a newspaper.
First things first. We first need to load the jpeg in memory such that we can process it as a matrix. Let's do this with the jpeg package. We also need to rotate the matrix as the jpeg uses a different origin for plotting:
Transforming the image to the frequency domain in R is actually ridiculously simple. Just one call to the fft function is enough. The multiplication with the ‘checkerboard’ is to ensure that the frequency spectrum is centred in the middle. The call to the fft function results in a matrix with the same dimensions as the original image. However, the matrix contains complex numbers. A complex number is composed of a real part and an imaginary part, which can be represented as a vector on two axes. The length of this vector (the modulus, calculated with the Mod function) gives the magnitude of frequencies in the spectrum. The angle of the vector (the argument, calculated with the Arg function) gives the phase shift of the sinusoidal waves. The latter is not very interesting for filtering of the image, but crucial in the reconstruction of the image after filtering.
When the magnitude of frequencies is plotted, the spectrum as shown below is obtained. Peaks in the spectrum are shown as white areas. The white blur in the centre of the spectrum contains much of the information of the original image. Other peaks circled in red are the result of the regular distortions in the image. Using the locator function we can now manually select the circled peaks with mouse clicks, so we have their location in the spectrum.
Next step is to filter out the selected (red circled) peaks from the spectrum. Basically, what we need to do is set the magnitude to zero at the location of the peaks. I've written a fancy function that does that, only slightly more complicated. The function multiplies the selected location with a bell-shaped curve (with the bell upside down). If you don't like this fancy approach, you set the ‘gaussian’ argument in the function to FALSE. When you do that, all values with a radius of r around the selected location in the spectrum are set to zero. We're all done filtering the unwanted peaks from the spectrum.
Note that we modified the magnitudes in the frequency spectrum. We left the phase shifts unchanged. We need to reconstruct the complex values from the modified magnitudes and the phase. This is not complex at all using the complex function (pardon the awful wordplay). Next we can apply the inverse Fourier transform, by simply setting the parameter ‘inv’ to TRUE in the fft function.
Just like that we have our filtered image which we can save with the writeJPEG function:
|
|
 |
| Before filter | After filter |
I know, the image is a bit blurry, noisy and a bit pale, but at least we got rid of the nasty distortion (which I put there myself in the first place). As you can imagine, there are plenty of other applications possible with Fourier transform. Maybe I'll save those for a later blog...
No comments:
Post a Comment