4. Cartesian perodic under-sampling along parallel lines

Here the goal is to illustrate the typical artifacts of standard deterministic regular (or periodic) undersampling along the phase encoding direction (here $k_y$) used in parallel imaging. Below we illustrate the following cases:

1. full Cartesian sampling $R= n/m = 1$ where $n=N^2$ is the image size, $N$ the image dimension and $m$ the number of measurements in k-space:
2. undersampling with a factor $R=2$
3. undersampling with a factor $R=4$
4. undersampling with a factor $R=8$

• Author: Philippe Ciuciu (philippe.ciuciu@cea.fr)
• Date: 06/24/2022
• Date: 04/02/2025 (brainweb compatibility)
• Target: IEEE EMBS-SPS Summer School on Novel acquisition and image reconstruction strategies in accelerated Magnetic Resonance Imaging

#DISPLAY BRAIN PHANTOM
%matplotlib inline

import numpy as np
import os.path as op
import os
import math ; import cmath
import matplotlib.pyplot as plt
import sys

from skimage import data, io, filters

import pywt as pw
import matplotlib.pyplot as plt

import brainweb_dl as bwdl

#Previousmly we used numerical phantoms
if 0:
    cwd = os.getcwd()
    dirimg_2d = op.join(cwd,"..", "data")
    img_size = 512   #256
    FOV = 0.2 #field of view in meters
    pixelSize = FOV/img_size
    #load data file corresponding to the target resolution
    filename = "BrainPhantom" + str(img_size) + ".png"
    mri_filename = op.join(dirimg_2d, filename)
    mri_img = io.imread(mri_filename, as_gray=True)
    plt.figure()
    plt.title("Brain Phantom, size = "+ str(img_size))
    if mri_img.ndim == 2:
        plt.imshow(mri_img, cmap=plt.cm.gray)
    else:
        plt.imshow(mri_img)
    plt.show()

plt.rcParams["image.origin"]="lower"
plt.rcParams["image.cmap"]='Greys_r'

mri_img = bwdl.get_mri(4, "T1")[70, ...].astype(np.float32)
#mri_img = bwdl.get_mri(5, "T2")[150, ...].astype(np.float32)
print(mri_img.shape)
img_size = mri_img.shape[0]
plt.figure()
plt.imshow(abs(mri_img))
plt.title("Original brain image")
plt.show()

(256, 256)

kspace_mask_full = np.ones((img_size, img_size), dtype="float64")

#import numpy.fft as fft
norm = "ortho"
def fft(x):
    return np.fft.fft2(x, norm=norm)

def ifft(x):
    return np.fft.ifft2(x, norm=norm)

# Generate the subsampled kspace with R=2
kspace_data = np.fft.fftshift(fft(mri_img))   # put the 0-freq in the middle of axes as

# Generate the kspace data: first Fourier transform the image
kspace_data = np.fft.fftshift(fft(mri_img))
#add Gaussian complex-valued random noise
signoise = 10
#kspace_data += np.random.randn(*mri_img.shape) * signoise * (1+1j)
# Simulate independent noise realization on the real & imag parts
kspace_data += (np.random.randn(*mri_img.shape) + 1j * np.random.randn(*mri_img.shape)) * signoise
# Mask data to perform subsampling
kspace_data *= kspace_mask_full

# Zero order solution
image_rec0 = ifft(np.fft.ifftshift(kspace_data))

fig, axs = plt.subplots(1, 3, figsize=(8, 8) )
axs[0].imshow(kspace_mask_full, cmap='gray_r')
axs[0].set_title("Full Cartesian mask (R=1)")
axs[1].imshow(np.abs(kspace_data),  cmap='gray_r', vmax=.01*np.abs(kspace_data).max())
axs[1].set_title("Masked data")
axs[2].imshow(np.abs(image_rec0), cmap='gray')
axs[2].set_title("Cartesian recon")

Text(0.5, 1.0, 'Cartesian recon')

import numpy.matlib as mlib

# generate Cartesian lines in a straightforward manner
#a = (np.linspace(0,img_size,img_size+1))/img_size -0.5    # work in normalized frequency
r2 = (int)(img_size/2)
r4 = (int)(img_size/4)
r8 = (int)(img_size/8)
print("2-fold undersampling, m= ", r2)
print("4-fold undersampling, m= ", r4)
print("8-fold undersampling, m= ", r8)

selected_ksp_line = np.ones((1, img_size), dtype="float64")
skipped_ksp_line = np.zeros((1, img_size), dtype="float64")
k_space_pattern_r2 = np.concatenate((selected_ksp_line, skipped_ksp_line), axis=0)
kspace_mask_r2 = np.tile(k_space_pattern_r2, (r2, 1))

#k_space_pattern_r4 = np.concatenate((selected_ksp_line, skipped_ksp_line, skipped_ksp_line,skipped_ksp_line), axis=0)
k_space_pattern_r4 = np.concatenate((selected_ksp_line, np.tile(skipped_ksp_line, (3,1))), axis=0)
kspace_mask_r4 = np.tile(k_space_pattern_r4, (r4, 1))

k_space_pattern_r8 = np.concatenate((selected_ksp_line, np.tile(skipped_ksp_line, (7,1))), axis=0)
kspace_mask_r8 = np.tile(k_space_pattern_r8, (r8, 1))


fig, axs = plt.subplots(1, 3, figsize=(16, 16) )
axs[0].imshow(kspace_mask_r2) #, cmap='Greys_r'
axs[0].set_title("Cartesian regular under-sampling mask (R=2)")
axs[1].imshow(kspace_mask_r4, cmap='Greys_r')
axs[1].set_title("Cartesian regular under-sampling mask (R=4)")
axs[2].imshow(kspace_mask_r8, cmap='Greys_r')
axs[2].set_title("Cartesian regular under-sampling mask (R=8)")

2-fold undersampling, m=  128
4-fold undersampling, m=  64
8-fold undersampling, m=  32

Text(0.5, 1.0, 'Cartesian regular under-sampling mask (R=8)')

Remark 0

$\leadsto$  

R = 2 : under-sample 1 line out of every 2

R = 4 : under-sample 1 line out of every 4

R = 8 : under-sample 1 line out of every 8

• Generate undersampled data for $R=2$ and perform image reconstruction
• What do you observe?

# Generate the kspace data: first Fourier transform the image
kspace_data_r2 = np.fft.fftshift(fft(mri_img))
#add Gaussian complex-valued random noise
signoise = 10
# Simulate independent noise realization on the real & imag parts
kspace_data_r2 += (np.random.randn(*mri_img.shape) + 1j * np.random.randn(*mri_img.shape)) * signoise

# Mask data to perform subsampling
kspace_data_r2 *= kspace_mask_r2

# Zero order image reconstruction
image_rec0_r2 = ifft(np.fft.ifftshift(kspace_data_r2))

fig, axs = plt.subplots(2, 2, figsize=(10, 10) )
axs[0,0].imshow(mri_img, cmap='Greys_r')
axs[0,0].set_title("True image")
axs[0,1].imshow(kspace_mask_r2, cmap='Greys_r')
axs[0,1].set_title("Sampling mask")
axs[1,0].imshow(np.abs(kspace_data_r2),  cmap='gray', vmax=0.01*np.abs(kspace_data_r2).max())
#axs[1].imshow(np.abs(np.fft.ifftshift(kspace_data)), cmap='Greys_r')
axs[1,0].set_title("k-space noisy data (R=2)")
axs[1,1].imshow(np.abs(image_rec0_r2), cmap='gray')
axs[1,1].set_title("Zero-order recon")
plt.show()

$\leadsto$  Since R = 2 (under-sampling 1 line out of every 2), we can observe that in the reconstructed image np.abs(image_rec0_r4), the noyau spatially repeats 2 times. The image contains artifacts, where the artifact at the top of the reconstructed image is the lower half of the noyau, while the artifact at the bottom is the upper half of the noyau.

If we vertically concatenate two copies of the reconstructed image, we can observe spatial continuity, revealing the periodic nature of the artifacts. This is essentially because kspace_data_r2 is multiplied by kspace_mask_r2 in k-space, which corresponds to a convolution in the spatial domain.

• Generate undersampled data for $R=4$ and perform image reconstruction
• What do you observe?

# Generate the kspace data: first Fourier transform the image
kspace_data_r4 = np.fft.fftshift(fft(mri_img))
#add Gaussian complex-valued random noise
signoise = 10
# Simulate independent noise realization on the real & imag parts
kspace_data_r4 += (np.random.randn(*mri_img.shape) + 1j * np.random.randn(*mri_img.shape)) * signoise
# Mask data to perform subsampling
kspace_data_r4 *= kspace_mask_r4

# Zero order image reconstruction
image_rec0_r4 = ifft(np.fft.ifftshift(kspace_data_r4))

fig, axs = plt.subplots(2, 2, figsize=(10, 10) )
axs[0,0].imshow(mri_img, cmap='Greys_r')
axs[0,0].set_title("True image")
axs[0,1].imshow(kspace_mask_r4, cmap='Greys_r')
axs[0,1].set_title("Sampling mask")
axs[1,0].imshow(np.abs(kspace_data_r4),  cmap='gray', vmax=0.01*np.abs(kspace_data_r4).max())
axs[1,0].set_title("k-space noisy data (USF=2)")
axs[1,1].imshow(np.abs(image_rec0_r4), cmap='Greys_r')
axs[1,1].set_title("Zero-order recon")
plt.show()

$\leadsto$  Since R = 4 (under-sampling 1 line out of every 4), we can observe that in the reconstructed image np.abs(image_rec0_r2), the noyau spatially repeats 4 times. The image contains artifacts, where the artifact at the top of the reconstructed image is the lower half of the noyau, while the artifact at the bottom is the upper half of the noyau.

If we vertically concatenate two copies of the reconstructed image, we can observe spatial continuity, revealing the periodic nature of the artifacts. This is essentially because kspace_data_r4 is multiplied by kspace_mask_r4 in k-space, which corresponds to a convolution in the spatial domain.

• Generate undersampled data for $R=8$ and perform image reconstruction
• What do you observe?

# Generate the kspace data: first Fourier transform the image
kspace_data_r8 = np.fft.fftshift(fft(mri_img))
#add Gaussian complex-valued random noise
signoise = 10
# Simulate independent noise realization on the real & imag parts
kspace_data_r8 += (np.random.randn(*mri_img.shape) + 1j * np.random.randn(*mri_img.shape)) * signoise
# Mask data to perform subsampling
kspace_data_r8 *= kspace_mask_r8

# Zero order image reconstruction
image_rec0_r8 = ifft(np.fft.ifftshift(kspace_data_r8))

fig, axs = plt.subplots(2, 2, figsize=(10, 10) )
axs[0,0].imshow(mri_img, cmap='Greys_r')
axs[0,0].set_title("True image")
axs[0,1].imshow(kspace_mask_r8, cmap='Greys_r')
axs[0,1].set_title("Sampling mask")
axs[1,0].imshow(np.abs(kspace_data_r8),  cmap='gray', vmax=0.01*np.abs(kspace_data_r4).max())
axs[1,0].set_title("k-space noisy data (USF=2)")
axs[1,1].imshow(np.abs(image_rec0_r8), cmap='Greys_r')
axs[1,1].set_title("Zero-order recon")
plt.show()

$\leadsto$  Since R = 8 (under-sampling 1 line out of every 8), we can observe that in the reconstructed image np.abs(image_rec0_r8), the noyau spatially repeats 8 times. The image contains artifacts, where the artifact at the top of the reconstructed image is the lower half of the noyau, while the artifact at the bottom is the upper half of the noyau.

If we vertically concatenate two copies of the reconstructed image, we can observe spatial continuity, revealing the periodic nature of the artifacts. This is essentially because kspace_data_r8 is multiplied by kspace_mask_r8 in k-space, which corresponds to a convolution in the spatial domain.

QUESTION:

• Do you know what key ingredient may help to recover the reference image pretty well while still using these regular under-sampling patterns?

$\leadsto$ We should try to find additional information (sparsity) or add some regularizations (constraints), like the method Compressed Sensing .

Q1  

• For the different under-sampling scenarios (factor R, cf. Block 3: r2, r4 and r8), comment the effect of increasing R.

$\leadsto$ As R increases, the noyau repeats more frequently in the spatial domain, and the quality of the reconstructed image degrades.

Q2  

• Vary the noise level and measure the impact on image quality (SSIM score).

def Reconstruction (mask, signoise) :
    '''
    vector_mask = [2, 4, 8]
    vector_signoise = [0, 10, 20, 100]
    '''
    if type(mask) == int or float :
        r = (int)(img_size/mask)

        selected_ksp_line = np.ones((1, img_size), dtype="float64")
        skipped_ksp_line = np.zeros((1, img_size), dtype="float64")

        k_space_pattern_r = np.concatenate((selected_ksp_line, np.tile(skipped_ksp_line, (mask-1,1))), axis=0)
        kspace_mask_r = np.tile(k_space_pattern_r, (r, 1))
        # k_space_pattern_r8 = np.concatenate((selected_ksp_line, np.tile(skipped_ksp_line, (7,1))), axis=0)
        # kspace_mask_r8 = np.tile(k_space_pattern_r8, (r8, 1))

        # Generate the kspace data: first Fourier transform the image
        kspace_data_r = np.fft.fftshift(fft(mri_img))
        #add Gaussian complex-valued random noise
        # signoise = 10
        # Simulate independent noise realization on the real & imag parts
        kspace_data_r += (np.random.randn(*mri_img.shape) + 1j * np.random.randn(*mri_img.shape)) * signoise
        # Mask data to perform subsampling
        kspace_data_r *= kspace_mask_r

        # Zero order image reconstruction
        image_rec0_r = ifft(np.fft.ifftshift(kspace_data_r))

        return kspace_data_r, image_rec0_r
        

$\leadsto$   Here we write a function using a mask based on standard deterministic regular (or periodic) undersampling along the phase encoding direction (here $k_y$).

We under-sample 1 line out of every R .

We apply this mask to the k-space data for subsampling. Finally, the reconstructed image is obtained by applying the inverse Fourier transform to the masked k-space data.

vector_mask = [2, 4, 8]
vector_signoise = [0, 10, 20, 100]

'''
vector_mask = [2, 4, 8]
vector_signoise = [0, 10, 20, 100]
'''
fig, axs = plt.subplots(len(vector_signoise), 2, figsize=(12, 12))

for i in range(len(vector_signoise)) :
    for j in range(2):

        mask = vector_mask[0]
        signoise = vector_signoise[i]

        kspace_data, image_rec0 = Reconstruction (mask, signoise)

        # 显示被噪声影响的 k-space 数据
        axs[i, 0].imshow(np.abs(kspace_data), cmap='gray', vmax=0.005*np.abs(kspace_data).max())
        
        # 显示重建的图像
        axs[i, 1].imshow(np.abs(image_rec0), cmap='Greys_r')

        axs[i, j].set_xlabel(f"signoise = {signoise}")
        axs[i, j].xaxis.set_label_position('top')  # 移动 xlabel 到顶部
        axs[i, j].set_ylabel(f"mask R = {mask}")
    

axs[0, 0].set_title("k-space noisy data")
axs[0, 1].set_title("Zero-order recon")

# Add a main title (suptitle) to the figure
fig.suptitle(f"k-space noisy data and its Reconstruction Results for Mask: R = {mask}", fontsize=16)

# Display the figure with adjusted layout
plt.tight_layout(rect=[0, 0, 1, 0.96])  # Adjust layout to fit suptitle
plt.show()

$\leadsto$   In the case where mask = 2 (R=2), we plot 4 reconstructed images with signoise varying in vector_signoise = [0, 10, 20, 100] :

• Initially, as signoise increases, the quality of the reconstructed image degrades.

• To quantitatively evaluate the reconstruction accuracy, we can calculate the Structural Similarity Index between the true image and the reconstructed image.

$\leadsto$   We can compute the Structural Similarity Index to evaluate the quality of reconstruction :

$$ SSIM(x,y) = \frac{(2m_x m_y + C_1)(2\Sigma_{xy} + C_2)}{(m_x^2 + m_y^2 + C_1)(\Sigma_{xx}^2 + \Sigma_{yy}^2 + C_2)} $$

where $m$ represents the mean, and $\Sigma$ denotes the covariance, and $C_1$ and $C_2$ are constants.

The result represents the overall similarity between the two images, ranging from [0, 1], where 1 indicates they are identical, and 0 means there is no similarity at all.

In Python, import the module :

from skimage.metrics import structural_similarity as ssim

Then, simply use the following line to compute the SSIM score and map :

ssim_score, ssim_map = ssim(mri_img, np.abs(image_rec0), full=True)

from skimage.metrics import structural_similarity as ssim

vector_mask = [2, 4, 8]
vector_signoise = np.linspace(0, 100, 100)

'''
vector_mask = [2, 4, 8]
vector_signoise = [0, 10, 20, 100]
'''

vector_ssim_score = []
vector_image_rec0 = []


for i in range(len(vector_signoise)) :
    

        mask = vector_mask[0]
        signoise = vector_signoise[i]

        kspace_data, image_rec0 = Reconstruction (mask, signoise)

        vector_image_rec0.append(image_rec0)
        ssim_score, ssim_map = ssim(mri_img, np.abs(image_rec0), full=True)
        vector_ssim_score.append(ssim_score)
    
    
plt.figure(figsize=(10, 6))
plt.plot(vector_signoise, vector_ssim_score, linestyle='-', color='red', marker='.', markersize=8, markerfacecolor='red', label='SSIM Score')

plt.xlabel('Signoise')
plt.ylabel('SSIM Score')
plt.title('SSIM Score vs Signoise')
plt.grid(True)
plt.legend()
plt.show()

$\leadsto$   Remember, we add a Gaussian complex-valued random noise by kspace_data_r += (np.random.randn(*mri_img.shape) + 1j * np.random.randn(*mri_img.shape)) * signoise .

The results show that the maximum SSIM Score occurs at signoise = 0.00 . Below, we display the reconstructed images for it :

• signoise = 0.00 (corresponding to the maximum SSIM Score)

In this scenario, we fixed R = 2 .

# Find the index and value of the ；maximum Frobenius norm.
max_index = np.argmax(vector_ssim_score)
image_corresponding_max_index = vector_image_rec0[max_index]
signoise_corresponding_max_index = vector_signoise[max_index]
# print(signoise_corresponding_max_index)

fig, axes = plt.subplots(1, 2, figsize=(8, 4))

axes[0].imshow(np.abs(image_corresponding_max_index), cmap='Greys_r')
axes[0].set_title(f'Signoise (for Max SSIM Score) : {signoise_corresponding_max_index:.2f}')
axes[0].axis('off')

axes[1].imshow(mri_img, cmap='Greys_r')
axes[1].set_title("True image")
axes[1].axis('off')

plt.tight_layout()
plt.show()

• When we select signoise = 0.00 , the best reconstructed image is obtained in terms of the SSIM Score.