2009 Blair Bohannan

Back to Psych 204 Projects 2009

Project Title - MATLAB image processing tutorial

This page accompanies a MATLAB tutorial which simulates some artifacts that might be found in MR images. Their causes, and in some cases corrections, are demonstrated in the code. One- and two-dimensional examples are provided.
Image artifacts in this tutorial include DC Offset, Quadrature Ghosting, RF Noise, Gradient, and Frequency-Sampling Jitter. Information on the first four topics can be found in Hornak, Joseph P. The Basics of MRI, Chapter 11. More information on the mathematical derivations can be found in Smith, Julius O. Mathematics of the Discrete Fourier Transform (DFT). A full reference list can be found at the end of this page.

About the tutorial

This tutorial consists of a MATLAB (.m) script and some image files. Unzip the files to a folder and either add to path or set the MATLAB workspace to this folder before starting the tutorial.
The tutorial uses cells, and each section below corresponds to a cell in the tutorial. To evaluate a cell in the tutorial, first navigate into it. Ctrl+Enter evaluates the current cell, and Ctrl+Shift+Enter evaluates and moves to the next cell. Variables are re-used throughout the tutorial, so start at the first cell and work through the tutorial in order.
The tutorial package contains a selection of image files. Those titled "stimulus##.tif" (including the one used on this page) are from an image set prepared by Roozbeh Kiani. The associated publication can be found here. Used with permission.

Tutorial

Cell I. Loading images for analysis

The first cell simply loads an image, converts it to grayscale, and sets the size of the image. nPix, the variable setting the number of pixels in the x and y dimensions, should be a power of 2 for faster FFT calculation.
This figure contains the image in grayscale. It has nPix x nPix pixels (in this case, 256 x 256).

Cell II. Using fftshift

The fftshift function is helpful for visualing Fourier transforms in a familiar format. In the 1-dimensional N-point Fourier transform (like the transform of an audio signal), it takes points N/2+1 to N (containing the Nyquist and negative frequency components) and moves them to the beginning of the array. See Oppenheim and Schafer, Chapter 4, if you want to know why this can happen. Thus, the fftshifted function causes the poisitive and negative frequencies to be symmetric about the zero-frequency point. This being the case, it's helpful to keep in mind that the zero-frequency sample is at sample N/2+1.
Here is a 1-dimensional example. The red dot is the zero-frequency component of the transform of the sinusoid in the top subplot.

In the 2-dimensional case, fftshift swaps quadrants 1 and 3, and also 2 and 4 of the Fourier transform. Thus, the low-frequency information is moved to the center of the image, with the [0,0] frequency point moving from matrix index [1,1] to [nPix/2+1,nPix/2+1]. The fftshifted Fourier transform of an image is what we are used to looking at when we talk about k-space. For the purpose of this tutorial, time domain is the 1-dimensional case of image domain, and k-space is specifically the 2-dimensional case of frequency domain.
Here is a 2-d Fourier transform without fftshift. Note the bright spots (high amplitudes) in the outer corners of the image; these are the low frequencies.

Here is the same transform, with fftshift applied to it. Now the low-frequency, high-amplitude information is in the center of the image.
Note: fftshift should be used after performing the FFT. Example:

X=fftshift(fft(x));

Using ifftshift

The ifftshift function reverses the effect of fftshift (use this instead of performing fftshift twice). The ifftshift should occur before taking the ifft:

x=ifft(ifftshift(X));

Cell III. Artifact 1: DC offset

Derivation

The DC offset is a constant, positive shift in amplitude. It results in a large peak at frequency zero. This can be derived from the equation for the DFT. Recall that the DFT is defined as

$X_k=\sum _{n=0}^{N-1}{x}_n{e}^{-\frac{j2\pi }{N}kn}k=0,\dots ,N-1$

where $x$ is the time signal at times $n$ and $X$ is the Fourier transform at frequency indices $k$. The DFT of a signal $x$ plus a constant DC offset $\Delta$ would then be

$X_k=\sum _{n=0}^{N-1}(x_n+\Delta )e^{-\frac{j2\pi }{N}kn}k=0,\dots ,N-1$

Since this is a linear system, we can separate the DC component from the original signal.

$X_k=\sum _{n=0}^{N-1}{x}_n{e}^{-\frac{j2\pi }{N}kn}+\sum _{n=0}^{N-1}{\Delta }_n{e}^{-\frac{j2\pi }{N}kn}k=0,\dots ,N-1$

We are interested in what is going on at frequency $k=0$. Assuming that the first term does not have some other DC artifact, we know that its zero-frequency component should be zero (this has to do with orthogonality of the DFT sinusoids). Focusing on the second term, since we are considering frequency $k=0$, the exponential is equal to $1$ for all $n$. Thus, the amplitude of the transform at this frequency is the sum of every term in the DC signal.

$X_0=N\Delta$

One-dimensional example

Here is a DC signal and its Fourier transform. Note that the DC amplitude is 3 and the FFT length is 32, making the zero-frequency peak 96.

Here is the same DC signal added to the sinusoid from the previous section. Note that while the sinusoid was previously centered around 0, now it is centered around 3, the DC amplitude. The Fourier transform of this signal is the Fourier transform from before, with the added peak at zero.

Two-dimensional example

The same will be the case for 2-dimensional data. The image will have a bright spot in the center (at frequency [0,0]). In the next figure, the zebra image has been combined with a large DC offset and frequency transformed. The original Fourier transform has been subtracted from this artifact Fourier transform, and we are left with only the bright spot.
Figure 7

Cell IV. Removing time-domain DC offset

Two quick solutions come to mind for removing the time-domain (or image-domain) DC offset:

• Subtract the mean of the signal from every sample;
• Use a narrow highpass filter to remove the zero frequency component.

Here is the 1-dimensional example from the previous section. The zero-frequency component has simply been set to zero (so it has effectively been filtered out). The IFFT of this new transform has shifted the original signal back down.

Cell V. Frequency-domain DC offset

Derivation

DC offset can occur in the frequency domain as well, with similar effects. Note, however, that the equation for the IDFT contains a scaling factor $\frac{1}{N}$ before the summation:

$x_n=\frac{1}{N}\sum _{k=0}^{N-1}{X}_k{e}^{\frac{j2\pi }{N}kn}n=0,\dots ,N-1$

(this would be $\frac{1}{N^2}$ for the 2-dimensional case)
Therefore, the peak at time $n=0$ will simply be the DC amplitude $\Delta$, not $N\Delta$.

One-dimensional example

Here is the frequency-domain offset for a 1-dimensional signal. The top subplot shows the sinusoid's Fourier transform, shifted up by a constant amount. The bottom subplot shows the inverse Fourier transform. The signal at time 0 now has the DC amplitude added to what was already present.

As in the time-domain case, this artifact can be removed by subtracting the mean of the Fourier transform from every element (in the Fourier transform). Note in the figure below, though, that doing so removes the entire signal at time zero, so the original signal at the time needs to be recovered (this can be done with some form of interpolation).

Cell VI. Frequency-domain DC offset - two-dimensional case

Here is an example of the effect of a k-space DC offset on a reconstructed image. The original zebra image has been reduced to a smaller pixel set so that the artifact is easier to see. The top two plots show the original image and its Fourier transform. The bottom two plots show the reconstruction of the original Fourier transform and the reconstruction of a Fourier transform that has a large DC offset. Note that the high amplitude of the bright spot has obscured the rest of the image.