Download R Code

In this post, I show how to remove some 'noisy' patterns in an image. I will use 'simple' mathematical transformations (Fourier transform) that allow to efficiently process images.

  rm(list=ls()) # clear memory
  if (!require(viridis)) install.packages('viridis')

## Loading required package: viridis

## Loading required package: viridisLite

I recommend to read this lecture (from the Robotics Research Group at Oxford) to know more about FFT in image processing. One should know that performing operations in the Fourier space is generally easier and less computationally demanding. For instance, a convolution is a simple product in Fourier space. The image I have chosen is an 8 bits image (so can have a maximum of 255 and a minimum of 0). It consists of an array that I can easily plot without using additional libraries (just download it here).

  # Load image + Matrix transform + rotation
  image0 <- read.csv("IM1.txt", header = F, sep = '\t')
  # To display this array as an image, I have to convert it to a matrix
  image0 <- data.matrix(image0) 
  # Rotation (not needed)
  image0 <- t(apply(image0, 2, rev))
  
  # Show image (this is a 8 bits image and so there are 2^8 elements)
  image(image0, axes = F, col = grey(seq(0, 1, length = 256)), asp = 1)

We cant then perform the a 2-dimensional discrete Fast Fourier Transform (FFT) using :

  # FFT on the image0
  fft_image0 <- (fft(image0, inverse = F))
  # Display
  # Here the max of Mod(fft_image0) is close to 1E7
  colorBreaks <- 10^(seq(0, 7, by = 0.1)) # set color breaks
  image(Mod(fft_image0), breaks = colorBreaks, col = (cividis(length(colorBreaks) - 1L)),asp=1)

One should keep in mind that the FFT (its magnitude) is not a 8 bits image. In other words, using col = grey(seq(0, 1, length = 256 would result in an image that displays no information at all. Here, one can notice some features (crosses) that most probably orginate from the diagonal stripes. To remove those crosses I will create a 256 x 256 matrix (FDfilter_matrix) with some zero squares (appearing in black in the image below).

  # Create the filter
  N <- ncol(image0)
  inx_filter <- 26
  FDfilter_matrix <- data.matrix(matrix(1:1, nrow = N, ncol = N ))
  FDfilter_matrix[inx_filter:1, (N-inx_filter):N] <- 0 
  FDfilter_matrix[(N-inx_filter):N, inx_filter:1] <- 0
  
  # Display
  image(FDfilter_matrix, axes = F, col = grey(seq(0, 1, length = 255)), asp = 1)

  # Apply the filter and display
  fft_image0_filtered = FDfilter_matrix * fft_image0
  image(Mod(fft_image0_filtered), breaks = colorBreaks, col = rev(cividis(length(colorBreaks) - 1L)),asp =1)

The image above represents the pointwise product of the FFT (magnitude) with the filter. To get back the original (filtered) image, we have to perform an inverse Fourier Transform. As you can see, the result is not perfect but most of the stripes have been removed.

  # Invert the FFT
  inv_fft <- (fft(fft_image0_filtered, inverse = T))
  image(Mod(inv_fft), axes = F, col = grey(seq(0, 1, length = 256)), asp = 1)