The file is located in the folder "data"
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
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
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)
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}
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)
Before that there is a pre-proocessing we have to do.
kspace_shifted = np.fft.fftshift(kspace)
plt.imshow(np.abs(kspace_shifted),cmap='Reds',clim=(0,15e5))
plt.colorbar()
plt.title('Shifted (magn) k-space',fontsize=20)
img = np.fft.fft2(kspace_shifted / np.sqrt(Nx*Ny),norm="ortho")
Visualize the image (remember, it is again a complex image)
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)
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
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);
the corresponding pixels coordinates are
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)
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.
plt.imshow(np.abs(kspace_low_frequencies))
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
typical Field Of View (FOV) is 220mm..
Each pixel side is 220/112 mm ..
# 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')
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)
kspace_lf_shifted = np.fft.fftshift(kspace_low_frequencies)
img_lf = np.fft.fft2(kspace_lf_shifted / np.sqrt(Nx*Ny),norm="ortho")
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 "
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..
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')
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))
kspace_hf_shifted = np.fft.fftshift(kspace_high_frequencies)
img_hf = np.fft.fft2(kspace_hf_shifted / np.sqrt(Nx*Ny),norm="ortho")
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)
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)
We keep the whole k-space, but we set a point of the k-space to a MUCH MUCH bigger value...
# 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)
# # 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)
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
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')
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.
# 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)