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:
| #required packages: | |
| require(jpeg) | |
| #Read a jpeg from a website and store as matrix "mat_clown" | |
| con <- file("http://4.bp.blogspot.com/-uObQLjPG8Xg/VceaZNyI-NI/AAAAAAAAx0o/ZlcTTB1pomU/s1600/fourier_efteling.jpg", "rb") | |
| #make sure 'n' is big enough to hold your entire file: | |
| jpg_clown <- readBin(con, "raw", n = 3e4) | |
| close(con) | |
| rm(con) | |
| mat_clown <- readJPEG(jpg_clown) | |
| rm(jpg_clown) | |
| #rotate the matrix to its intended position | |
| rotate <- function(x) t(apply(x, 2, rev)) | |
| mat_clown <- rotate(mat_clown) |
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.
| #checker board pattern used to centre the Fourier spectrum: | |
| checkerboard <- (-1)^(row(mat_clown) + col(mat_clown)) | |
| #modify matrix with checkerboard pattern: | |
| mat_checker <- mat_clown*checkerboard | |
| #Discrete fast Fourier transform: | |
| mat_fft <- fft(mat_checker) | |
| #magnitude in frequency domain | |
| #Modulus = length of the vector of the complex number: | |
| magni <- Mod(mat_fft) | |
| #phase in the frequency domain | |
| #Argument = angle of the vector of the complex number: | |
| phase <- Arg(mat_fft) | |
| #Show magnitude in frequency domain (= spectrum): | |
| par(mar = c(0,0,0,0), oma = c(0,0,0,0)) | |
| image(1:nrow(magni), 1:ncol(magni), log(magni + 1), col = grey((0:255)/255)) |
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.
| #select peaks in spectrum to filter out | |
| #click on each peak, then end the locator | |
| peaks <- data.frame(locator()) | |
| #function to create a filter for circle shaped peaks in the spectrum | |
| #dimensions: dimensions of the spectrum matrix | |
| #x: x position of the centre of the circle to be filtered out | |
| #y: y position of ... | |
| #r: radius of circle to be filtered out of the spectrum | |
| #gaussian: flag that indicates whether the circle should have a Gaussian density or not | |
| # when set to TRUE, smoother results are obtained. | |
| filter_circle <- function(dimensions, x, y, r = 2, gaussian = T) | |
| { | |
| result <- matrix(1, dimensions[1], dimensions[2]) | |
| distance <- sqrt((row(result) - x)^2 + (col(result) - y)^2) | |
| if (gaussian) | |
| { | |
| result <- dnorm(distance, sd = r*2) | |
| result <- 1 - result/max(result) | |
| } else | |
| { | |
| result[distance < r] <- 0 | |
| } | |
| return(result) | |
| } | |
| jpeg("fourier_spectrum.jpg", nrow(magni), ncol(magni)) | |
| par(mar = c(0,0,0,0), oma = c(0,0,0,0)) | |
| image(1:nrow(magni), 1:ncol(magni), log(magni + 1), col = grey((0:255)/255), bty = "n", xaxt = "n", yaxt = "n") | |
| points(y~x, peaks, col = "red", cex = 3) | |
| dev.off() | |
| #create filter | |
| filt <- by(peaks, 1:nrow(peaks), function(x) filter_circle(dim(magni), x$x, x$y, 14)) | |
| #function that returns the product of all elements in a list | |
| list_prod <- function(x) | |
| { | |
| if (length(x) > 1) | |
| { | |
| result <- x[2:length(x)] | |
| result[[1]] <- result[[1]]*x[[1]] | |
| return(list_prod(result)) | |
| } else return(x) | |
| } | |
| filt <- list_prod(filt)[[1]] | |
| #filter out the selected peaks in the spectrum: | |
| magni_filter <- magni*filt |
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.
| #reconstruct fourier transformed image from filtered spectrum: | |
| re_construct <- cos(phase)*magni_filter | |
| im_construct <- sin(phase)*magni_filter | |
| mat_fft_construct <- matrix(complex(real = re_construct, imaginary = im_construct), nrow(mat_clown), ncol(mat_clown)) | |
| #Apply inversed Fourier transform to obtain the filtered image: | |
| mat_clown_constr <- Re(fft(mat_fft_construct, inv = T))*checkerboard/prod(dim(mat_clown)) | |
| #normalize the resulting matrix between 0 and 1 | |
| mat_clown_constr <- (mat_clown_constr - min(mat_clown_constr)) / | |
| (max(mat_clown_constr) - min(mat_clown_constr)) | |
| #save image after filtering: | |
| con <- file("clown_filtered.jpg", "wb") | |
| clown_raw <- writeJPEG(rotate(rotate(rotate(mat_clown_constr))), raw(), quality = 0.7) | |
| writeBin(clown_raw, con) | |
| close(con) |
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