Load data from the file "oneslice.nii"

The file is located in the folder "data"

In [ ]:
import nibabel as nib
import os

folder_path = './data'
file_name = 'oneslice.nii'

loaded_image = nib.load(os.path.join(folder_path,file_name))
kspace = loaded_image.get_data()

Print the image's dimensions

In [ ]:
Nx = kspace.shape[0]
Ny = kspace.shape[1]

print("kspace x size is " + str(Nx))
print("kspace y size is " + str(Ny))

Visualize the kspace: remember that it is a complex image

In [ ]:
import matplotlib.pyplot as plt
%matplotlib inline

plt.figure(num=None, figsize=(12, 6), dpi=80, facecolor='w', edgecolor='k')
plt.subplot(121)
plt.imshow(kspace.real,cmap='seismic')
plt.colorbar(shrink=0.5)
plt.title('real part',fontsize=20)
head_length = 5.
head_width = 3.
plt.arrow(Nx/2., Ny, 0, -Ny+head_length, head_width=head_width, head_length=head_length, fc='k', ec='k')
plt.arrow(0., Ny/2., Nx-head_length, 0., head_width=head_width, head_length=head_length, fc='k', ec='k')
plt.xlabel('$K_x$',fontsize=20)
plt.ylabel('$K_y$',fontsize=20)
plt.tick_params(
    axis='both',          # changes apply to the x-axis
    which='both',         # both major and minor ticks are affected
    left='off',
    right='off',
    bottom='off',         # ticks along the bottom edge are off
    top='off',            # ticks along the top edge are off
    labelbottom='off',    # labels along the bottom edge are off
    labelleft='off') 

plt.subplot(122)
plt.imshow(kspace.imag,cmap='seismic')
plt.colorbar(shrink=0.5)
plt.title('imaginary part',fontsize=20)
plt.arrow(Nx/2., Ny, 0, -Ny+head_length, head_width=head_width, head_length=head_length, fc='k', ec='k')
plt.arrow(0., Ny/2., Nx-head_length, 0., head_width=head_width, head_length=head_length, fc='k', ec='k')
plt.xlabel('$K_x$',fontsize=20)
plt.ylabel('$K_y$',fontsize=20)
plt.tick_params(
    axis='both',          # changes apply to the x-axis
    which='both',         # both major and minor ticks are affected
    left='off',
    right='off',
    bottom='off',         # ticks along the bottom edge are off
    top='off',            # ticks along the top edge are off
    labelbottom='off',    # labels along the bottom edge are off
    labelleft='off') 

plt.subplots_adjust(left=None, bottom=None, right=None, top=None,
                wspace=0.4, hspace=None)

Visualize the kspace in magnitude and phase

The kspace is a complex image K \begin{equation} K(x,y) = \Re[{K(x,y)}] + i\cdot \Im[{K(x,y)}] \end{equation}

it's magnitude is \begin{equation} |K(x,y)| = \sqrt{\Re[{K(x,y)}]^2 + \Im[{K(x,y)}]^2} \end{equation}

it's phase is \begin{equation} \Phi(x,y) = arctan \left( \frac{\Im[K(x,y)]}{\Re[K(x,y)]} \right) \end{equation}

In [ ]:
import numpy as np

plt.figure(num=None, figsize=(12, 6), dpi=80, facecolor='w', edgecolor='k')
plt.subplot(121)
plt.imshow(np.abs(kspace),cmap='Reds')
plt.colorbar(shrink=0.5)
plt.title('magnitude',fontsize=20)
head_length = 5.
head_width = 3.
plt.arrow(Nx/2., Ny, 0, -Ny+head_length, head_width=head_width, head_length=head_length, fc='k', ec='k')
plt.arrow(0., Ny/2., Nx-head_length, 0., head_width=head_width, head_length=head_length, fc='k', ec='k')
plt.xlabel('$K_x$',fontsize=20)
plt.ylabel('$K_y$',fontsize=20)
plt.tick_params(
    axis='both',          # changes apply to the x-axis
    which='both',         # both major and minor ticks are affected
    left='off',
    right='off',
    bottom='off',         # ticks along the bottom edge are off
    top='off',            # ticks along the top edge are off
    labelbottom='off',    # labels along the bottom edge are off
    labelleft='off') 

plt.subplot(122)
plt.imshow(np.angle(kspace),cmap='seismic',clim=(-np.pi,np.pi))
plt.colorbar(shrink=0.5)
plt.title('phase',fontsize=20)
plt.arrow(Nx/2., Ny, 0, -Ny+head_length, head_width=head_width, head_length=head_length, fc='k', ec='k')
plt.arrow(0., Ny/2., Nx-head_length, 0., head_width=head_width, head_length=head_length, fc='k', ec='k')
plt.xlabel('$K_x$',fontsize=20)
plt.ylabel('$K_y$',fontsize=20)
plt.tick_params(
    axis='both',          # changes apply to the x-axis
    which='both',         # both major and minor ticks are affected
    left='off',
    right='off',
    bottom='off',         # ticks along the bottom edge are off
    top='off',            # ticks along the top edge are off
    labelbottom='off',    # labels along the bottom edge are off
    labelleft='off') 

plt.subplots_adjust(left=None, bottom=None, right=None, top=None,
                wspace=0.4, hspace=None)

Apply the 2D Fourier transform to pass into "image space"

Before that there is a pre-proocessing we have to do.

  • use the function fftshift (np.fft.fftshift?) to shift the zero frequency to the corners of the image (this step is required because of the implementation of the Fourier transform
In [ ]:
kspace_shifted = np.fft.fftshift(kspace)
  • visualize the shifted k-space to understand
In [ ]:
plt.imshow(np.abs(kspace_shifted),cmap='Reds',clim=(0,15e5))
plt.colorbar()
plt.title('Shifted (magn) k-space',fontsize=20)
  • apply the 2D Fourier Transform to the shifted k-space (Note, we have to multiply the k-space by the root square of the number of samples, again because of implementation issues of the Fast Fourier Transform)
In [ ]:
img = np.fft.fft2(kspace_shifted / np.sqrt(Nx*Ny),norm="ortho")

Visualize the image (remember, it is again a complex image)

In [ ]:
plt.figure(num=None, figsize=(12, 6), dpi=80, facecolor='w', edgecolor='k')
plt.subplot(121)
plt.imshow(img.real,cmap='seismic',clim=(-20000,20000))
plt.colorbar(shrink=0.5)
plt.title('Image real part',fontsize=20)
plt.xlabel('x pixels',fontsize=20)
plt.ylabel('y pixels',fontsize=20)

plt.subplot(122)
plt.imshow(img.imag,cmap='seismic',clim=(-20000,20000))
plt.colorbar(shrink=0.5)
plt.title('Image imaginary part',fontsize=20)
plt.xlabel('x pixels',fontsize=20)
plt.ylabel('y pixels',fontsize=20)

plt.subplots_adjust(left=None, bottom=None, right=None, top=None,
                wspace=0.4, hspace=None)

Visualize the Image as magnitude and phase

In [ ]:
plt.figure(num=None, figsize=(12, 6), dpi=80, facecolor='w', edgecolor='k')
plt.subplot(121)
plt.imshow(np.abs(img),cmap='gray',clim=(0,22000))
plt.colorbar(shrink=0.5)
plt.title('Image\'s magnitude',fontsize=20)
plt.xlabel('x pixels',fontsize=20)
plt.ylabel('y pixels',fontsize=20)

plt.subplot(122)
plt.imshow(np.angle(img),cmap='seismic',clim=(-np.pi,np.pi))
plt.colorbar(shrink=0.5)
plt.title('Image\'s phase',fontsize=20)
plt.xlabel('x pixels',fontsize=20)
plt.ylabel('y pixels',fontsize=20)

plt.subplots_adjust(left=None, bottom=None, right=None, top=None,
                wspace=0.4, hspace=None)

The PIXEL intensity in the magnitude image IS the AMPLITUDE (or MAGNITUDE) of the magnetization vector for that location

Select a region in the k-space corresponding to the LOW frequencies

In [ ]:
import matplotlib.patches as patches

## This is the side of the the selected rectangle
x_half_side = 2
y_half_side = 2

fig, ax = plt.subplots(num=None,figsize=(8, 8))
ax.imshow(np.abs(kspace))
rect = patches.Rectangle((Nx/2-x_half_side,Ny/2-y_half_side),
                         2*x_half_side,
                         2*y_half_side,
                         linewidth=2,
                         edgecolor='r',
                         facecolor='none')
ax.add_patch(rect)
ax.set_title('Selected area of the k-space',fontsize=20);

Set to ZERO all the frequencies OUTSIDE the selected region

the corresponding pixels coordinates are

In [ ]:
x_start = int(Nx/2 -   x_half_side)
x_stop  = int(Nx/2 + 2*x_half_side)

y_start = int(Ny/2 -   y_half_side)
y_stop  = int(Ny/2 + 2*y_half_side)
In [ ]:
low_frequency_mask = np.zeros((Nx,Ny),dtype=bool)
low_frequency_mask[y_start:y_stop,x_start:x_stop] = True

kspace_low_frequencies = kspace.copy()
kspace_low_frequencies[~low_frequency_mask] = 0. + 1j*0.
In [ ]:
plt.imshow(np.abs(kspace_low_frequencies))

Plot the frequencies corresponding to the selected area of the k-space

In [ ]:
def calculate_frequency(dimension, wx, wy, fx, fy, phix, phiy, mask=None):
    '''
    dimension: tuple specifying size (size x,size y)
    wx:    x period length
    wy:    y period length
    fx:    x frequency
    phix:  x initial phase
    fy:    y frequency
    phiy:  y initial phase
    mask:  brain mask (optional)

    returns a n-array of size dimension with a sinusoidal phase between -pi and
    +pi
    '''
    if len(dimension) != 2:
        raise ValueError('dimension must be a tuple of length 2')

    sin_xy = np.zeros(dimension)
    for y in range(sin_xy.shape[0]): #use xrange if Python2
        for x in range(sin_xy.shape[1]): #use xrange if Python2
            if mask is None:
                sin_xy[x,y] = np.pi * np.sin(2*np.pi*fx*x/wx+phix+2*np.pi*fy*y/wy+phiy)
            else:
                if mask[x,y]==1:
                    sin_xy[x,y] = np.pi * np.sin(2*np.pi*fx*x/wx+phix+2*np.pi*fy*y/wy+phiy)
    return sin_xy

What is a spatial frequency?

typical Field Of View (FOV) is 220mm..

Each pixel side is 220/112 mm ..

In [ ]:
# The samallest resolution is
resolution = 0.22/112.   # 22 cm divided by the number of pixels
# The highest frequency in the image is the reciprocal of the resolution
K_max = 1. / resolution
print('the resolution is ' + str(resolution) + ' meters')
print('the maximum spatial frequency is ' + str(K_max) + ' 1/meters')
In [ ]:
frequency_mask = low_frequency_mask.copy()
frequency_to_show = np.zeros((112,112))
# frequency_mask = np.ones((112,112),dtype=bool)

plt.figure(num=None, figsize=(16, 16), dpi=80, facecolor='w', edgecolor='k')
image_number = 0
for y in range(112): #use xrange if Python2
    for x in range(112): #use xrange if Python2
        if frequency_mask[y,x] == True:
            image_number = image_number + 1
            K_y = K_max/2. * (y-Ny/2.) / (Ny/2.) 
            K_x = K_max/2. * (x-Nx/2.) / (Nx/2.)
           
            sin_xy = calculate_frequency((112,112), K_max, K_max, K_y, K_x, 0., 0., mask=None)
            
            if image_number>20:
                break
            
            plt.subplot(7,5,image_number)
            plt.imshow(sin_xy,origin=True)

Regenerate the images

In [ ]:
kspace_lf_shifted = np.fft.fftshift(kspace_low_frequencies)
img_lf = np.fft.fft2(kspace_lf_shifted / np.sqrt(Nx*Ny),norm="ortho")
In [ ]:
plt.figure(num=None, figsize=(12, 6), dpi=80, facecolor='w', edgecolor='k')
plt.subplot(121)
plt.imshow(np.abs(img),cmap='gray',clim=(0,22000))
plt.colorbar(shrink=0.5)
plt.title('Original magnitude',fontsize=20)
plt.xlabel('x pixels',fontsize=20)
plt.ylabel('y pixels',fontsize=20)

plt.subplot(122)
plt.imshow(np.abs(img_lf),cmap='gray',clim=(0,22000))
plt.colorbar(shrink=0.5)
plt.title('Low frequency magnitude',fontsize=20)
plt.xlabel('x pixels',fontsize=20)
plt.ylabel('y pixels',fontsize=20)

plt.subplots_adjust(left=None, bottom=None, right=None, top=None,
                wspace=0.4, hspace=None)

Now try to change the selected frequencies.. In particular try to select a larger area.. Example: before

" x_half_side = 2 y_half_side = 2 " now do

" x_half_side = 5 y_half_side = 5 "

Select region with HIGH (DIRECTIONAL) frequencies

Now, we use the variables "posx" and "posy" to move the center of the selected region far from the center of the k-space. Then play with these variables to see the effects..

In [ ]:
import matplotlib.patches as patches

x_half_side = 5
y_half_side = 5

# POSTITIONS WITH RESPECT TO THE CENTER OF THE K-SPACE
posx = 15
posy = 5

fig, ax = plt.subplots(num=None,figsize=(8, 8))
ax.imshow(np.abs(kspace))

rect = patches.Rectangle((Nx/2+posx,Ny/2-posy),
                         2*x_half_side,
                         2*y_half_side,
                         linewidth=2,
                         edgecolor='r',
                         facecolor='none')
ax.add_patch(rect)
ax.set_title('Selected area of the k-space',fontsize=20);
ax.set_xlabel('x')
ax.set_ylabel('y')

set to zero the other frequencies

In [ ]:
x_start = int(Nx/2+posx)
x_stop  = int(x_start + 2*x_half_side)

y_start = int(Ny/2-posy)
y_stop  = int(y_start + 2*y_half_side)

high_frequency_mask = np.zeros((Nx,Ny),dtype=bool)
high_frequency_mask[y_start:y_stop,x_start:x_stop] = True

kspace_high_frequencies = kspace.copy()
kspace_high_frequencies[~high_frequency_mask] = 0. + 1j*0.
plt.imshow(np.abs(kspace_high_frequencies))

Go in image space

In [ ]:
kspace_hf_shifted = np.fft.fftshift(kspace_high_frequencies)
img_hf = np.fft.fft2(kspace_hf_shifted / np.sqrt(Nx*Ny),norm="ortho")
In [ ]:
plt.figure(num=None, figsize=(12, 6), dpi=80, facecolor='w', edgecolor='k')
plt.subplot(121)
plt.imshow(np.abs(img),cmap='gray',clim=(0,22000))
plt.colorbar(shrink=0.5)
plt.title('Original magnitude',fontsize=20)
plt.xlabel('x pixels',fontsize=20)
plt.ylabel('y pixels',fontsize=20)

plt.subplot(122)
plt.imshow(np.abs(img_hf),cmap='gray')
plt.colorbar(shrink=0.5)
plt.title('High (directional) frequency magnitude',fontsize=20)
plt.xlabel('x pixels',fontsize=20)
plt.ylabel('y pixels',fontsize=20)

plt.subplots_adjust(left=None, bottom=None, right=None, top=None,
                wspace=0.4, hspace=None)
In [ ]:
frequency_mask = high_frequency_mask.copy()
frequency_to_show = np.zeros((112,112))
# frequency_mask = np.ones((112,112),dtype=bool)

plt.figure(num=None, figsize=(16, 16), dpi=80, facecolor='w', edgecolor='k')
image_number = 0
for y in range(112): #use xrange if Python2
    for x in range(112): #use xrange if Python2
        if frequency_mask[y,x] == True:
            image_number = image_number + 1
            K_y = K_max/2. * (y-Ny/2.) / (Ny/2.) 
            K_x = K_max/2. * (x-Nx/2.) / (Nx/2.)
           
            sin_xy = calculate_frequency((112,112), K_max, K_max, K_y, K_x, 0., 0., mask=None)
            
            if image_number>35:
                break
            
            plt.subplot(7,5,image_number)
            plt.imshow(sin_xy,origin=True)

To make it more clear.. let's do a different experiment now..

We keep the whole k-space, but we set a point of the k-space to a MUCH MUCH bigger value...

In [ ]:
# POSTITION, with respect to k-space center
posy = 10
posx = 0

# AMPLITUDE OF THE MODIFIED FREQUENCY
amplification_factor = 40000000.


kspace_modified = kspace.copy()

posy = int(posy+Ny/2.)
posx = int(posx+Nx/2.)

kspace_modified[posy,posx] = amplification_factor*(1. + 1j*1)
In [ ]:
# # POSTITION, with respect to k-space center
# posy = 23
# posx = 0

# # AMPLITUDE OF THE MODIFIED FREQUENCY
# amplification_factor = 40000000.

# posy = int(posy+Ny/2.)
# posx = int(posx+Nx/2.)

# kspace_modified[posy,posx] = amplification_factor*(1. + 1j*1)
In [ ]:
kspace_modified_shifted = np.fft.fftshift(kspace_modified)
img_modified = np.fft.fft2(kspace_modified_shifted / np.sqrt(Nx*Ny),norm="ortho")

visualize the frequency we are amplifying

In [ ]:
K_y = K_max/2. *(posy-Ny/2.) / (Ny/2.) 
K_x = K_max/2. *(posx-Nx/2.) / (Nx/2.)
           
sin_xy = calculate_frequency((112,112), K_max, K_max, K_y, K_x, 0., 0., mask=None)
plt.imshow((sin_xy),cmap='gray')
In [ ]:
plt.figure(num=None, figsize=(12, 6), dpi=80, facecolor='w', edgecolor='k')
plt.subplot(121)
plt.imshow(np.abs(img),cmap='gray',clim=(0,22000))
plt.colorbar(shrink=0.5)
plt.title('Original magnitude',fontsize=20)
plt.xlabel('x pixels',fontsize=20)
plt.ylabel('y pixels',fontsize=20)

plt.subplot(122)
plt.imshow(np.abs(img_modified),cmap='gray',clim=(0,22000))
plt.colorbar(shrink=0.5)
plt.title('Modified magnitude',fontsize=20)
plt.xlabel('x pixels',fontsize=20)
plt.ylabel('y pixels',fontsize=20)

plt.subplots_adjust(left=None, bottom=None, right=None, top=None,
                wspace=0.4, hspace=None)

Play with the parameters

posy = 10 posx = 10

amplification_factor = 20000000.

Visualize the difference between the 2 images

In [ ]:
# The true difference
difference = np.abs(img)-np.abs(img_modified)
max_diff = np.max(np.abs(difference))

# the difference scaled to pi
difference_pi = difference / max_diff * np.pi

plt.figure(num=None, figsize=(8, 8), dpi=80, facecolor='w', edgecolor='k')
plt.imshow(difference)
plt.colorbar(shrink=0.5)
# plt.title('Difference scaled to $\pi$',fontsize=20)
plt.title('Difference',fontsize=20)
plt.xlabel('x pixels',fontsize=20)
plt.ylabel('y pixels',fontsize=20)
In [ ]: