SOCR ≫ TCIU Website ≫ TCIU GitHub ≫

1 Fourier Transform of images

Let’s demonstrate the use of R’s fft() to calculate the 2D FT of images. Information on data acquisition frequency and block length (in sec) cannot be included into the fft() call.

R generates a single 2D vector of the same dimensions as the data containing a list of complex numbers. The function Mod(fft()) is used to extract the magnitudes of the Fourier coefficients, which are computed by: \[magnitude = \sqrt{(real * real + imaginary * imaginary)},\] where the real and imaginary components are extracted by functions Re(fft()) and Im(fft()), respectively.

The method fft() generates only meaningful frequency up to half the sampling frequency. The FT returns values of the discrete Fourier transform for both positive and negative frequencies. Although, as sampling a signal in discrete time intervals causes aliasing problems, R yields all frequencies up to the sampling frequency. For instance, sampling a 50 Hz sine wave and 950 Hz sine wave with 1000 Hz will generate identical results, as the FT cannot distinguish between the two frequencies. Hence, the sampling frequency must always be at least twice as high as the expected signal frequency. For each actual frequency in the signal, the FT will give 2 peaks (one at the “actual” frequency and one at sampling frequency minus “actual” frequency). This will make the second half of the magnitude vector a mirror image of the first half.

As long as the sampling frequency is at least twice as high as the expected signal frequency, all meaningful information will be contained in the the first half of the magnitude vector. However, a peak in the low frequency range might be present when high “noise” frequency is present in the signal (or image). At this point, the vector of extracted magnitudes is only indexed by the frequencies but has no associated frequencies. To calculate the corresponding frequencies, the FT simply takes (or generates) the index vector (1, 2, 3, …, length(magnitude vector)) and divides it by the length of the data block (in sec).

In 1D, the phases would represent a vector of the same length as the magnitude vector with the phases (0 to \(2\pi\) or \(-\pi\) to \(+\pi\)) of each frequency. Phase shifts are translations in space (e.g., x-axis) for a given wave component that are measured in angles (radians). For instance, shifting a wave \(f(x)=0.5sin(3wt)+0.25sin(10wt)\) by \(\frac{\pi}{2}\) would produce the following Fourier series:

\[f(t)=0.5sin\left (3wt+\frac{\pi}{2}\right )+0.25sin\left (10wt+\frac{\pi}{2}\right ).\]

In 2D, The Fourier transform (FT/IFT) for images is defined by: \[\hat{f}(u,v)=F(u,v)=\int_{-\infty}^{\infty}{\int_{-\infty}^{\infty}{f(x,y)e^{-i2\pi(ux+vy)}dxdy}},\] \[f(x,y)=\hat{\hat{f}}(x,y)=\hat{F}(x,y)=\int_{-\infty}^{\infty}{\int_{-\infty}^{\infty}{F(u,v)e^{i2\pi(ux+vy)}dudv}},\] where \(u\) and \(v\) are the spatial frequencies, \(F(u,v)=F_R(u,v)+iF_I(u,v)\) is a complex number for each pair of arguments, \[|F(u,v)|=\sqrt{F_R^2(u,v)+F_I^2(u,v)}\] is the magnitude spectrum, and \[arctan\left (\frac{F_I(u,v)}{F_R(u,v)}\right )\] is the phase angle spectrum.

The complex exponential \[e^{-i2\pi(ux+vy)}=cos(2\pi(ux+vy)) +i\ sin(2\pi(ux+vy))\] represents the real and imaginary (complex) sinusoidal terms in the 2D plane. The extrema of its real part (\(cos(2\pi(ux+vy))\)) occur at \(2\pi(ux+vy)=n\pi\). Using vector notation, \[2\pi(ux+vy)=2\pi \langle U, X \rangle =n\pi,\] where the extrema points \(U=(u,v)^T\) and \(X=(x,y)^T\) represent sets of equally spaced parallel lines with normal \(U\) and wavelength \(\frac{1}{\sqrt{u^2+v^2}}\).

Let’s define the index shifting paradigm associated with the discrete FT, which is simply used for convenience and better visualizaiton. It has no other relevance to the actual calculation of the FT and its inverse, IFT.

When applying the forward or reverse generalized discrete FT it is possible to shift the transform sampling in time and frequency domain by some real offset values, \(a,b\). Symbolically,

\[\hat{f}(k) = \sum_{n=0}^{N-1} f(n) e^{-\frac{2 \pi i}{N} (k+b) (n+a)} \quad \quad k = 0, \dots, N-1.\]

Note: Remember that in R, the array indices start with 1, not 0, as in some other languages.

The function fftshift() is useful for visualizing the Fourier transform with the zero-frequency component in the middle of the spectrum. Its inverse counterpart, ifftshift(), is needed to rearrange again the indices appropriately after the IFT is employed, so that the image is correctly reconstructed in spacetime. The FT only computes half of the frequency spectrum corresponding to the non-negative (positive and zero if the length(f) is odd) frequencies in order to save computation time. To preserve the dimensions of the output \(\hat{f}=FT(f)\), the second half of the frequency spectrum (the complex conjugate of the first half) is just added at the end of this vector. In a 1D setting, the result of fft() is:

\(0\ 1\ 2\ 3\ ...\ (freq\ bins > 0)\ ... {\frac{Fs}{2}}\) and \(-{\frac{Fs}{2}}\ ... \ (freq\ bins < 0)\ ...\ -3\ -2\ -1.\)

where \(F_s\) is the frequency sample. The fftshift method sets the zero-frequency component in the center of the array, i.e., it just shifts (offsets) the second part with the negative frequency bins to the beginning and the first part to the end of the resulting FT vector, or matrix. Thus, the shifted discrete FT can be nicely plotted in the center covering the frequency spectrum from \(-{\frac{Fs}{2}}\) on the left to \({\frac{Fs}{2}}\) on the right. This is not necessary, but is used for better visualization aesthetics. To synthesize back the correct image, after using fftshift on the FT signal, we always have to undo that re-indexing by using ifftshift() on the inverse-FT.

# FFT SHIFT
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")
  }
}

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")
  }
}

Now we will present the FT and IFT of a pair of synthetic 2D images (a square and a disk).