ECE 401 MP5: Discrete Fourier Transform

In this file, we're going to use the inverse DFT to create images from MRI spectral measurements. Then we'll use the DFT to filter out noise in an image.

import os, h5py
import numpy as np
import matplotlib.figure
import matplotlib.pyplot as plt
%matplotlib inline
import importlib
import submitted

0. Data

The data are provided for you in the file data.hdf5. Let's load it, and see what it contains.

with h5py.File('data.hdf5','r') as f:
    print('data.hdf5 has the keys',f.keys())
    mri_dft = f['mri_dft'][:]
    image = f['image'][:]
    noisy_image = f['noisy_image'][:]
print('mri_dft is an array of shape',mri_dft.shape,'and dtype',mri_dft.dtype)
print('image is an array of shape',image.shape,'and dtype',image.dtype)
print('noisy_image is an array of shape',noisy_image.shape,'and dtype',noisy_image.dtype)

data.hdf5 has the keys <KeysViewHDF5 ['image', 'mri_dft', 'noisy_image']>
mri_dft is an array of shape (1114, 962) and dtype complex128
image is an array of shape (213, 320, 3) and dtype float64
noisy_image is an array of shape (213, 320, 3) and dtype float64

1. MRI

MRI (magnetic resonance images) are acquired in the frequency domain. It works roughly like this:

1. Every atom in your body spins, all the time, emitting weak radio signals. A strong magnetic field will align a few percent of them all in the same direction, so that their radio signals can be measured.
2. If the external magnetic field slightly different in different locations, then the synchronized spin will be at slightly different frequencies in different locations. By holding this difference for a carefully controlled amount of time, we can set it up so that your atoms are spinning with slightly different phases in different locations. The phase at location $(n_1,n_2)$ might be tweaked to be $\phi=-\omega_1n_1-\omega_2n_2$.
3. If the density of hydrogen atoms at location $(n_1,n_2)$ is given by $x[n_1,n_2]$, then the total signal collected is

$$X(\omega_1,\omega_2) = \sum_{n_1}\sum_{n_2} x[n_1,n_2] e^{-j\omega_1 n_1}e^{-j\omega_2 n_2}$$

1. After measuring $X(\omega_1,\omega_2)$ at $N_1$ different values of $\omega_1$, and $N_2$ different values of $\omega_2$, then we compute the inverse DFT in order to calculate the image.

The image you are provided, mri_dft, is a set of complex-valued frequency samples from a real MRI image (a public domain image from Wikipedia, donated by Johns Hopkins University). Let's try taking the inverse DFT, to see if we can see the image.

fig,ax = plt.subplots(figsize=(10,8))
ax.imshow(np.maximum(0,np.real(np.fft.ifft2(mri_dft))), cmap='gray')

<matplotlib.image.AxesImage at 0x7fe5054120d0>

The aliasing we see, here, is caused by capturing only $N_1\times N_2$ frequency samples, but computing the DFT using a spatial range $M_1\times M_2$ such that $M_1>N_1$ and/or $M_2>N_2$. By doing that, we accidentally discover that the DFT is periodic in space, not just in frequency:

$$x[N+n] = \frac{1}{N}\sum_k X[k]e^{j\frac{2\pi k(n+N)}{N}} = \frac{1}{N}\sum_k X[k]e^{j\frac{2\pi kn}{N}} =x[n]$$

In order to solve this problem, we need to figure out how many of the frequency samples in mri_dft are actually nonzero-valued. To get an idea, let's try plotting the first row:

fig, ax = plt.subplots(figsize=(14,4))
ax.stem(np.abs(mri_dft[0,0:200]))

<StemContainer object of 3 artists>

It looks like maybe the odd-numbered spectral samples are zero, and only the even-numbered spectral samples are valid. Let's try taking the inverse fft of only the even-numbered samples.

fig,ax = plt.subplots(figsize=(10,8))
ax.imshow(np.real(np.fft.ifft2(mri_dft[::2,::2])), cmap='gray')

<matplotlib.image.AxesImage at 0x7fe50730dc40>

That's the right size, but now the whole image has been shifted by about $N_1/2$ rows, and by about $N_2/2$ columns! One way to solve this is by multiplying the DFT by a shift operator, before we inverse DFT:

$$X_{\mbox{corrected}}[k_1,k_2] = \begin{cases} X[k_1,k_2] e^{-j\frac{2\pi k_1d_1}{N_1}}e^{-j\frac{2\pi k_2d_2}{N_2}} & 0\le k_1<\frac{N_1}{2},~0\le k_2<\frac{N_2}{2}\\ X[k_1,k_2] e^{-j\frac{2\pi (k_1-N_1)d_1}{N_1}}e^{-j\frac{2\pi (k_2-N_2)d_2}{N_2}} &\frac{N_1}{2}\le k_1<N_1,~\frac{N_2}{2}\le k_2<N_2 \end{cases} $$

where $d_1$ is the number of rows that we want to shift, and $d_2$ is the number of columns that we want to shift. Notice that the second line is actually equal to the first line, in theory (because $e^{-j\frac{2\pi N_1d_1}{N_1}}=1$), but implementing it in the way shown above will help to avoid artifacts caused by floating-point roundoff error. Notice that the formula above only lists two of the four cases; I think you can figure out what the other two cases should be.

importlib.reload(submitted)
help(submitted.downsample_and_shift_dft2)

Help on function downsample_and_shift_dft2 in module submitted:

downsample_and_shift_dft2(oversampled_dft, downsampling_factor, row_shift, col_shift)
    Input: 
      oversampled_dft [M1,M2] - a 2d array containing the oversampled DFT of a grayscale image
      downsampling_factor (scalar) - the factor by which the DFT image is oversampled
      row_shift (scalar)  - the number of rows that the image should be shifted
      col_shift (scalar) - the number of columns that the image should be shifted
    Output: 
      image [M1/downsampling_factor, M2/downsampling_factor] - the real part of the inverse DFT
      of the valid frequency samples, shifted by the specified numbers of rows and columns.

importlib.reload(submitted)
N1, N2 = mri_dft.shape
print(N1,N2)
image = submitted.downsample_and_shift_dft2(mri_dft,2,-N1/4,-N2/4)

fig, ax = plt.subplots(figsize=(10,8))
ax.imshow(np.maximum(0,np.minimum(255,image)),cmap='gray')

1114 962

<matplotlib.image.AxesImage at 0x7fe50730f520>

2. Noisy Image

Filtering out noise is easier in the time domain, if the noise is wideband noise (meaning that there is noise in a wide range of frequencies). Sometimes, though, the noise is in a narrow range of frequencies. For example, consider this image (this is the image "Clifton Beach 5, South Arm, Tasmania, Australia" by JJ Harrison, distributed on Wikimedia under Gnu Free Documentation License):

with h5py.File('data.hdf5','r') as f:
    image = f['image'][:]
    noisy_image = f['noisy_image'][:]
fig, ax = plt.subplots(1,2,figsize=(14,8))
ax[0].imshow(image)
ax[1].imshow(noisy_image)

<matplotlib.image.AxesImage at 0x7fe4c0e4f430>