SOCR ≫ TCIU Website ≫ TCIU GitHub ≫

1 2D Fourier Transform and Image Synthesis

1.1 Functions Used

1.1.1 fftshift()

#' FFT SHIFT
#' This function is useful for visualizing the Fourier transform with the zero-frequency 
#' component in the middle of the spectrum.
#' 
#' @param img_ff A Fourier transform of a 1D signal, 2D image, or 3D volume.
#' @param dim Number of dimensions (-1, 1, 2, 3).
#' @return A properly shifted FT of the array.
#' 
fftshift <- function(img_ff, dim = -1) {

  rows <- dim(img_ff)[1]    
  cols <- dim(img_ff)[2]
  # planes <- dim(img_ff)[3]

  swap_up_down <- function(img_ff) {
    rows_half <- ceiling(rows/2)
    return(rbind(img_ff[((rows_half+1):rows), (1:cols)], img_ff[(1:rows_half), (1:cols)]))
  }

  swap_left_right <- function(img_ff) {
    cols_half <- ceiling(cols/2)
    return(cbind(img_ff[1:rows, ((cols_half+1):cols)], img_ff[1:rows, 1:cols_half]))
  }
  
  #swap_side2side <- function(img_ff) {
  #  planes_half <- ceiling(planes/2)
  #  return(cbind(img_ff[1:rows, 1:cols, ((planes_half+1):planes)], img_ff[1:rows, 1:cols, 1:planes_half]))
  #}

  if (dim == -1) {
    img_ff <- swap_up_down(img_ff)
    return(swap_left_right(img_ff))
  }
  else if (dim == 1) {
    return(swap_up_down(img_ff))
  }
  else if (dim == 2) {
    return(swap_left_right(img_ff))
  }
  else if (dim == 3) {
    # Use the `abind` package to bind along any dimension a pair of multi-dimensional arrays
    # install.packages("abind")
    library(abind)
    
    planes <- dim(img_ff)[3]
    rows_half <- ceiling(rows/2)
    cols_half <- ceiling(cols/2)
    planes_half <- ceiling(planes/2)
    
    img_ff <- abind(img_ff[((rows_half+1):rows), (1:cols), (1:planes)], 
                    img_ff[(1:rows_half), (1:cols), (1:planes)], along=1)
    img_ff <- abind(img_ff[1:rows, ((cols_half+1):cols), (1:planes)], 
                    img_ff[1:rows, 1:cols_half, (1:planes)], along=2)
    img_ff <- abind(img_ff[1:rows, 1:cols, ((planes_half+1):planes)], 
                    img_ff[1:rows, 1:cols, 1:planes_half], along=3)
    return(img_ff)
  }
  else {
    stop("Invalid dimension parameter")
  }
}

1.1.2 ifftshift()

#' IFFT SHIFT
#' This function is useful for moving back the zero-frequency component in the middle of the spectrum
#' back to (0,0,0).  It rearranges in reverse (relative to fftshift()) the indices appropriately,
#' so that the image can be correctly reconstructed by the IFT in spacetime
#' 
#' @param img_ff An Inverse Fourier transform of a 1D signal, 2D image, or 3D volume.
#' @param dim Number of dimensions (-1, 1, 2, 3).
#' @return A properly shifted IFT of the input array.
#' 
ifftshift <- function(img_ff, dim = -1) {

  rows <- dim(img_ff)[1]    
  cols <- dim(img_ff)[2]    

  swap_up_down <- function(img_ff) {
    rows_half <- floor(rows/2)
    return(rbind(img_ff[((rows_half+1):rows), (1:cols)], img_ff[(1:rows_half), (1:cols)]))
  }

  swap_left_right <- function(img_ff) {
    cols_half <- floor(cols/2)
    return(cbind(img_ff[1:rows, ((cols_half+1):cols)], img_ff[1:rows, 1:cols_half]))
  }

  if (dim == -1) {
    img_ff <- swap_left_right(img_ff)
    return(swap_up_down(img_ff))
  }
  else if (dim == 1) {
    return(swap_up_down(img_ff))
  }
  else if (dim == 2) {
    return(swap_left_right(img_ff))
  }
  else if (dim == 3) {
    # Use the `abind` package to bind along any dimension a pair of multi-dimensional arrays
    # install.packages("abind")
    library(abind)
    
    planes <- dim(img_ff)[3]
    rows_half <- floor(rows/2)
    cols_half <- floor(cols/2)
    planes_half <- floor(planes/2)
    
    img_ff <- abind(img_ff[1:rows, 1:cols, ((planes_half+1):planes)], 
                    img_ff[1:rows, 1:cols, 1:planes_half], along=3)
    img_ff <- abind(img_ff[1:rows, ((cols_half+1):cols), (1:planes)], 
                    img_ff[1:rows, 1:cols_half, (1:planes)], along=2)
    img_ff <- abind(img_ff[((rows_half+1):rows), (1:cols), (1:planes)], 
                    img_ff[(1:rows_half), (1:cols), (1:planes)], along=1)
    return(img_ff)
  }
  else {
    stop("Invalid dimension parameter")
  }
}

1.1.3 fftinv()

#' Implicitly Invert the FT (IFT)
#' This function does the IFT and scales appropriately the  result to ensure IFT(FT()) = I()
#' 
#' @param x FT of a dataset.
#' @return The IFT of the input array.
#'
fftinv <- function( x ) { 
  fft( x, inverse=TRUE ) / length( x ) 
  }

1.2 Figure 6.2

Let’s try to perform the same image synthesis to reconstruct the Cyrillic and English alphabet images.

img_CyrAlpha <- load.image("BulgAlpha.jpg") 
# plot(img_CyrAlpha, axes = F)  
# Grayscaled # img_gray <- im[ , , 1, 1]
img_CyrAlpha <- matrix(img_CyrAlpha, nrow = dim(img_CyrAlpha)[1], ncol = dim(img_CyrAlpha)[2])
# EBImage::display(img_CyrAlpha, title='Image', method = "raster")
img_EngAlpha <- load.image("EngAlpha.jpg")
img_EngAlpha <- matrix(img_EngAlpha, nrow = dim(img_EngAlpha)[1], ncol = dim(img_EngAlpha)[2])
# add an extra zero column at the end to make the 2D Alphabet images homologous: 411 353
img_EngAlpha <- cbind(img_EngAlpha, rep(0, dim(img_EngAlpha)[1])) 
# dim(img_CyrAlpha); dim(img_EngAlpha)

# FFT
ft_img_CyrAlpha <- fft(img_CyrAlpha)  # fftw2d # Display Re(FT): display(fftshift(ft_img_CyrAlpha))
ft_img_EngAlpha <- fft(img_EngAlpha)  # display(fftshift(ft_img_EngAlpha))

# Magnitude and Phase
mag_ft_img_CyrAlpha <- sqrt(Re(ft_img_CyrAlpha)^2+Im(ft_img_CyrAlpha)^2)
mag_ft_img_EngAlpha <- sqrt(Re(ft_img_EngAlpha)^2+Im(ft_img_EngAlpha)^2)

# Phase  <- atan(Im(img_ff)/Re(img_ff))
phase_ft_img_CyrAlpha  <- atan2(Im(ft_img_CyrAlpha), Re(ft_img_CyrAlpha))
phase_ft_img_EngAlpha  <- atan2(Im(ft_img_EngAlpha), Re(ft_img_EngAlpha))

# FFT SHIFT
shift_ft_img_CyrAlpha <- fftshift(mag_ft_img_CyrAlpha)
shift_ft_img_EngAlpha <- fftshift(mag_ft_img_EngAlpha)
# Display FT
EBImage::display(log(shift_ft_img_CyrAlpha), title="FT Magnitude (Cyrillic Alphabet)")