1. Basic deterministic under-sampling

In this first notebook, we start with simple deterministic under-sampling using either a low-frequency or a high-frequency mask. Hence we want to generate a sampling mask \(\Omega\in\{0,1\}^N\) such that \(y = F_{\Omega} x\) where \(x\in R^{N}\) is the original image. The data \(y\) is then corrupted by additive noise \(b\):

\[ y = F_{\Omega} x + b \]

where \(b\) is a complex-valued circular Gaussian noise of variance \(\sigma^2\) and \(F_{\Omega} = \Omega F\) in this Cartesian setting (i.e. sampling on the grid). We consider:

1. a low-frequency Cartesian mask \(\Omega\) is defined by the central lines of k-space.

2. a low-frequency Cartesian mask \(\Omega\) is defined by a square box in k-space centered at \((k_x,k_y)=(0,0)\).

3. a high-frequency Cartesian mask \(\Omega\) is defined as the complementary set of the lines defined in to 1)

4. a high-frequency Cartesian mask \(\Omega\) is defined as the complementary set of the box defined in 2)

Based on these masks we generate the corresponding measured data from a reference MR image and finally perform Cartesian image reconstruction \(\widehat{x}\) from the under-sampled data \(y\) as follows:

\[ \widehat{x} = \tilde{F}^* y \]

where \(F^*\) is the zero-filled inverse Fourier transform (ifft).

• Author: Philippe Ciuciu (philippe.ciuciu@cea.fr)

• Date: 02/03/2025

• Target: IEEE EMBS-SPS Summer School on Novel acquisition and image reconstruction strategies in accelerated Magnetic Resonance Imaging

#DISPLAY BRAIN PHANTOM
%matplotlib inline
import sys
import os.path as op
import os
import math
import cmath

import numpy as np

import matplotlib.pyplot as plt

from skimage import data, io
import brainweb_dl as bwdl

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(4, "T2")[120, ...].astype(np.float32)
print(mri_img.shape)
img_size = mri_img.shape[0]
plt.figure()
plt.imshow(abs(mri_img))
#plt.imshow(abs(mri_2D),origin="lower")
#plt.imshow(abs(mri_2D),origin="upper")
plt.title("Original brain image")
plt.show()

/volatile/github-ci-mind-inria/gpu_mind_runner/_work/mri-acq-recon-book/mri-acq-recon-book/venv/lib/python3.10/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm

(256, 256)

• To get started, the objective is to build a low-frequency sampling mask that consists of the central k-space lines

• Fourier transform the reference MR image to compute k-space data, add zero-mean Gaussian complex-valued noise

• Mask the data using the above defined mask

• Next, perform zero-filled MR image reconstruction from masked k-space data

• Study the impact of the number of lines in the mask on the final image quality

# Generate trivial Cartesian sampling masks
#img_size = 512

mask="low_res"
factor_ctr = 4

#Objective: construct a k-space mask that consists of the central lines
# 0 entries: points not sampled, 1 entries: points to be sampled
# Initialize k-space with 0 everywhere and then place the "1" appropriately
kspace_masklines = np.zeros((img_size,img_size), dtype="float64")

low_res_size = img_size // factor_ctr + 1
idx_vec = np.linspace(img_size // 2 - low_res_size // 2, img_size // 2 + low_res_size // 2, low_res_size)
idx_vec_ =  idx_vec.astype("int")
#print(idx_vec_)

#Use fancy indexing to select low frequency lines only
kspace_masklines[idx_vec_, ] = 1


if 0:
    plt.figure()
    plt.imshow(kspace_masklines, cmap='Greys_r')
    plt.title("Sampling mask")
    plt.show()

norm = "ortho"
#norm = None

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

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

# 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_masklines

# Zero order solution
image_rec0 = ifft(np.fft.ifftshift(kspace_data))            # get back to the convention

fig, axs = plt.subplots(2, 2, figsize=(10, 10))
axs[0, 0].imshow(mri_img)
axs[0, 0].set_title("True image")
axs[0, 1].imshow(kspace_masklines)
axs[0, 1].set_title("Low frequency sampling mask")
axs[1, 0].imshow(np.abs(kspace_data), vmax=0.01*np.abs(kspace_data).max())
axs[1, 0].set_title("k-space noisy data")
axs[1, 1].imshow(np.abs(image_rec0))
axs[1, 1].set_title("Zero-order recon")
plt.show()

• Question: based on the previous example, try to construct a low-freqency sampling mask defined by a square box centered around \((k_x, k_y) = (0,0)\).

• Then, replicate the same steps:

  1. Generate noisy masked data

  2. Perform zero-filled MR image reconstruction

  3. Visualize results and study the impact of both the noise level and the mask size

#Objective: construct a k-space mask consisting of a central box
# Initialize k-space with 0 everywhere and then place the "1" appropriately
kspace_maskbox = np.zeros((img_size,img_size), dtype="float64")
# use fancy indexing along rows
kspace_maskbox[idx_vec_, ] = 1

list_img_size = np.arange(0., img_size).tolist()
filtered_center = [x for x in list_img_size if x not in idx_vec_]
array_idx_center = np.array(filtered_center)
array_idx_center_ = array_idx_center.astype("int")
# use fancy indexing along cols
kspace_maskbox[:, array_idx_center_] = 0

# 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
# 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_maskbox

fig, axs = plt.subplots(2, 2, figsize=(10, 10))
axs[0, 0].imshow(mri_img)
axs[0, 0].set_title("True image")
axs[0, 1].imshow(kspace_maskbox)
axs[0, 1].set_title("Low frequency sampling mask")
axs[1, 0].imshow(np.abs(kspace_data),  vmax=0.01*np.abs(kspace_data).max())
axs[1, 0].set_title("k-space noisy data")
axs[1, 1].imshow(np.abs(image_rec0))
axs[1, 1].set_title("Zero-order recon")
plt.show()

• Next question: Construct a k-space sampling mask that consists of the high-frequency lines, just discard the central lines

• Study the impact of the number of lines removed

# 0 entries: points not sampled, 1 entries: points to be sampled
# Initialize k-space with "1" everywhere and then place the "0" appropriately

kspace_masklines = np.ones((img_size,img_size), dtype="float64")
low_res_size = img_size // factor_ctr + 1
idx_vec_upper = np.linspace(img_size // 2 - low_res_size // 2, img_size // 2 + low_res_size // 2, low_res_size)
idx_vec_upper_ =  idx_vec.astype("int")

# set central k-space lines to 0: use fancy indexing along rows
kspace_masklines[idx_vec_upper_, ] = 0

# 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_masklines

# Zero order solution
image_rec0 = ifft(np.fft.ifftshift(kspace_data))            # get back to the convention

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

Question: Based on the previous example, try to construct a k-space mask that consists of the removing the central box centered in \((𝑘_𝑥,𝑘_𝑦)=(0,0)\).

• Then, replicate the same steps:

  1. Generate noisy masked data

  2. Perform zero-filled MR image recon

  3. Visualize results and study the impact of both the noise level and the mask size

# 0 entries: points not sampled, 1 entries: points to be sampled
# Initialize k-space with "1" everywhere and then place the "0" appropriately
kspace_maskbox = np.ones((img_size,img_size), dtype="float64")
# use fancy indexing along rows
kspace_maskbox[idx_vec_, ] = 0

list_img_size = np.arange(0., img_size).tolist()
filtered_center = [x for x in list_img_size if x not in idx_vec_]
array_idx_center = np.array(filtered_center)
array_idx_center_ = array_idx_center.astype("int")
# use fancy indexing along cols
kspace_maskbox[:, array_idx_center_] = 1

# 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_maskbox

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

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_maskbox, cmap='Greys_r')
axs[0,1].set_title("High frequency box sampling mask")
axs[1,0].imshow(np.abs(kspace_data),  cmap='gray', vmax=1*np.abs(kspace_data).max())
#axs[1].imshow(np.abs(np.fft.ifftshift(kspace_data)), cmap='Greys_r')
axs[1,0].set_title("k-space noisy data")
axs[1,1].imshow(np.abs(image_rec0), cmap='Greys_r')
axs[1,1].set_title("Zero-order recon")
plt.show()