Image reconstruction with Cartesian GRAPPA
Contents
Image reconstruction with Cartesian GRAPPA¶
In this tutoral we will read in a 2D cine dataset from ocmr (2D+time). This data is a Cartesian bSSFP cine which has been aquired with Partial Fourier acqustion in the kx (readout) direction. We will reconstruct a fully sampled k-space using TensorFlowMRI, and then simulate GRAPPA undersampling
We will then take this undersampled k-space and reconstreuct using GRAPPA in TensorFlow MRI
Set up TensorFlow MRI¶
If you have not yet installed TensorFlow MRI in your environment, you may do so
now using pip:
%pip install --quiet tensorflow-mri
# Upgrade Matplotlib. Versions older than 3.5.x may cause an error below.
%pip install --quiet --upgrade matplotlib
WARNING: You are using pip version 20.2.4; however, version 24.3.1 is available.
You should consider upgrading via the '/usr/local/bin/python -m pip install --upgrade pip' command.
Note: you may need to restart the kernel to use updated packages.
WARNING: You are using pip version 20.2.4; however, version 24.3.1 is available.
You should consider upgrading via the '/usr/local/bin/python -m pip install --upgrade pip' command.
Note: you may need to restart the kernel to use updated packages.
Then, import the package into your program to get started:
import tensorflow_mri as tfmri
print("TensorFlow MRI version:", tfmri.__version__)
TensorFlow MRI version: 0.22.0
We will also need a few additional packages:
from glob import glob
import tensorflow as tf
import matplotlib.pyplot as plt
import numpy as np
%pip install -U matplotlib
import matplotlib.pyplot as plt
Defaulting to user installation because normal site-packages is not writeable
Requirement already up-to-date: matplotlib in /usr/local/lib/python3.8/dist-packages (3.7.5)
Requirement already satisfied, skipping upgrade: pillow>=6.2.0 in /usr/local/lib/python3.8/dist-packages (from matplotlib) (10.4.0)
Requirement already satisfied, skipping upgrade: pyparsing>=2.3.1 in /usr/local/lib/python3.8/dist-packages (from matplotlib) (3.0.9)
Requirement already satisfied, skipping upgrade: fonttools>=4.22.0 in /usr/local/lib/python3.8/dist-packages (from matplotlib) (4.55.3)
Requirement already satisfied, skipping upgrade: python-dateutil>=2.7 in /usr/local/lib/python3.8/dist-packages (from matplotlib) (2.9.0.post0)
Requirement already satisfied, skipping upgrade: cycler>=0.10 in /usr/local/lib/python3.8/dist-packages (from matplotlib) (0.12.1)
Requirement already satisfied, skipping upgrade: contourpy>=1.0.1 in /usr/local/lib/python3.8/dist-packages (from matplotlib) (1.1.1)
Requirement already satisfied, skipping upgrade: numpy<2,>=1.20 in /usr/local/lib/python3.8/dist-packages (from matplotlib) (1.23.2)
Requirement already satisfied, skipping upgrade: importlib-resources>=3.2.0; python_version < "3.10" in /usr/local/lib/python3.8/dist-packages (from matplotlib) (6.4.5)
Requirement already satisfied, skipping upgrade: kiwisolver>=1.0.1 in /usr/local/lib/python3.8/dist-packages (from matplotlib) (1.4.7)
Requirement already satisfied, skipping upgrade: packaging>=20.0 in /usr/local/lib/python3.8/dist-packages (from matplotlib) (21.3)
Requirement already satisfied, skipping upgrade: six>=1.5 in /usr/lib/python3/dist-packages (from python-dateutil>=2.7->matplotlib) (1.14.0)
Requirement already satisfied, skipping upgrade: zipp>=3.1.0; python_version < "3.10" in /usr/local/lib/python3.8/dist-packages (from importlib-resources>=3.2.0; python_version < "3.10"->matplotlib) (3.8.1)
WARNING: You are using pip version 20.2.4; however, version 24.3.1 is available.
You should consider upgrading via the '/usr/local/bin/python -m pip install --upgrade pip' command.
Note: you may need to restart the kernel to use updated packages.
Using a GPU¶
TensorFlow MRI supports CPU and GPU computation. If there is a GPU available in your environment and it is visible to TensorFlow, it will be used automatically.
Tip
In Google Colab, you can enable GPU computation by clicking on Runtime > Change runtime type and selecting GPU under Hardware accelerator.
Tip
You can control whether CPU or GPU is used for a particular operation via
the tf.device
context manager.
Prepare the data¶
We will be using an example brain dataset from the ISMRM Reproducibility Challenge 1. Let’s download it.
cardiac_cine_data_filename = 'fs_0005_1_5T.h5'
cardiac_cine_data_url = "https://ocmr.s3.us-east-2.amazonaws.com/data/fs_0005_1_5T.h5."
!wget --quiet -O {cardiac_cine_data_filename} {cardiac_cine_data_url}
/bin/bash: wget: command not found
You may need to install ‘ismrmrd-python’ and ‘ismrmrd-python-tools’. This can be done like this:
Install ismrmrd-python from here: https://github.com/ismrmrd/ismrmrd-python
%%bash
if [ $(pip list --disable-pip-version-check | grep -c -w 'ismrmrd ') == 0 ]
then
git clone https://github.com/ismrmrd/ismrmrd-python.git /var/tmp/ismrmrd-python
python -m pip install --disable-pip-version-check -r /var/tmp/ismrmrd-python/requirements.txt
python -m pip install --disable-pip-version-check /var/tmp/ismrmrd-python
rm -rf /var/tmp/ismrmrd-python
fi
pip list --disable-pip-version-check | grep -w 'ismrmrd '
Install ocmr reader from here: https://raw.githubusercontent.com/MRIOSU/OCMR/master/Python/read_ocmr.py
%%bash
wget https://raw.githubusercontent.com/MRIOSU/OCMR/master/Python/read_ocmr.py
This dataset contains a fully sampled cartesian dataset. from the OCMR dataset: https://github.com/MRIOSU/OCMR/blob/master/Python/example_ocmr.ipynb The data is stored in a HDF5 file, which we can read using h5py. The downloaded file also has the sampling locations or k-space trajectory, so we do not need to calculate it.
# read_ocmr is downloaded from here:
# https://raw.githubusercontent.com/MRIOSU/OCMR/master/Python/read_ocmr.py
import read_ocmr as read
kData,param = read.read_ocmr(cardiac_cine_data_filename)
print('Dimension of kData: ', kData.shape)
kData = np.squeeze(kData)
print('Dimension of kData: ', kData.shape)
# kx, ky, ch, phase
# Reverse the order of the dimensions.
# [kx, ky, ch, phase] -> [phase, ch, kx, ky ]
kspace = np.transpose(kData, [3,2,0,1])
print(kspace.shape)
#(18, 18, 512, 208)
#[phase, ch, kx, ky ]
Imaging acquisition starts acq 0
Dimension of kData: (512, 208, 1, 18, 18, 1, 1, 1, 1)
Dimension of kData: (512, 208, 18, 18)
(18, 18, 512, 208)
Lets view the k-space data
fig1 = plt.figure(1, figsize=(8, 12)); fig1.suptitle("Original Fully Sampled K-Space and Partial Fourier K-space", fontsize=14);
tmp = plt.imshow(np.abs(np.squeeze(kspace[0,0,:,:])), aspect= 'auto');
plt.xlabel('ky');plt.ylabel('kx'); tmp.set_clim(0.0,0.01*np.max(np.abs(kspace))) # ky by kx
You can see that this k-space was acquired with Partial Foutier in the readout, kx, direction (only 404 of the 512 data points are acquired) We will reconstruct this using the code from the Partial Fourier Tutorial to get a fully sampled k-space
# The tensorflow functions for Partial Fourier Reconstruction takes in only the k-space which was acquired
# kspace should only contain the observed data, without zero-filling of any kind.
# Therefore in this case we need to remove all phase encode lines which are all zeros
# Find rows that contain only zeros
kspace_PF = kspace[:, :, np.any(kspace != 0, axis=(0, 1, 3)),:]
print('kspace_PF.shape: ', kspace_PF.shape)
# (18, 18, 404, 208)
# The tensowflow function needs to know what proportion if k-space was acquired
partial_fourier_factor = kspace_PF.shape[2] / kspace.shape[2]
print('partial_fourier_factor: ', partial_fourier_factor)
# 0.7890625
kspace_PF_flipped = np.flip(kspace_PF, axis=-2)
# Now do a partial fourier reconstruction seperately each time point seperately
# tfmri can do all time point together but it is less likey to run out of memory if we process one time point at a time
Reconstructed_Image = []
for t in range(np.shape(kspace_PF_flipped)[0]):
kspace_PF_single_phase = np.squeeze(kspace_PF_flipped[t,:,:,:])
#print('kspace_PF_single_phase.shape: ', kspace_PF_single_phase.shape)
#(18, 404, 208)
recon_im = tfmri.recon.partial_fourier(
kspace_PF_single_phase, [partial_fourier_factor, 1.0], method="pocs",
preserve_phase=True)
recon_im = np.flip(recon_im, axis=-2)
Reconstructed_Image.append(recon_im)
final_cine_img = np.array(Reconstructed_Image)
fully_sampled_kspace = tfmri.signal.fft(final_cine_img, axes=[-2, -1], norm='ortho', shift=True)
kspace_PF.shape: (18, 18, 404, 208)
partial_fourier_factor: 0.7890625
# Here we simulate the 2x undersampling (by removing every other line of k-space), with 20 ACS lines remaining in the centre of k-space
print('fully_sampled_kspace.shape: ', fully_sampled_kspace.shape)
#(18, 18, 512, 208)
#[phase, ch, kx, ky ]
# Here we work out which lines we will keep and which we will miss out: simulate GRAPPA acquistion
centre=int(fully_sampled_kspace.shape[3]/2)
calib_data = fully_sampled_kspace[:,:,:,(centre-10):(centre+10)]
print('calib.shape: ', calib_data.shape)
# (18, 18, 512, 20)
undersampling_mask = [[1 if ((y % 2 == 0) or (y<(centre+10) and y>(centre-10)) ) else 0 for x in range(fully_sampled_kspace.shape[2])] for y in range(kspace.shape[3])]
undersampling_mask = np.transpose(undersampling_mask)
us_mask = np.array(undersampling_mask)
print('us_mask.shape: ', us_mask.shape)
# [phases, coils, lines, columns]
#(512, 208)
plt.imshow(us_mask, cmap="viridis")
plt.show()
# Number of replicas you want in the third dimension
nTimePts = kspace.shape[0]
# Replicate the 2D array to form a 2D+t array
us_mask_2d_plus_t = np.stack([us_mask] * nTimePts, axis=0)
print(us_mask_2d_plus_t.shape)
#(18, 512, 208)
fully_sampled_kspace.shape: (18, 18, 512, 208)
calib.shape: (18, 18, 512, 20)
us_mask.shape: (512, 208)
(18, 512, 208)
# Now we have created our undersampling mask, we need to multiply the fully sampled data by it to get our synthetially undersampled data
mask_all_coils = np.stack([us_mask_2d_plus_t] * np.shape(kspace)[1], axis=1)
# Now do the undersampling / masking
kspace_US = np.multiply(fully_sampled_kspace, mask_all_coils)
print(np.shape(kspace_US))
#(18, 18, 512, 208)
#[phase, ch, kx, ky ]
fig1 = plt.figure(1, figsize=(8, 12)); fig1.suptitle("US K-space", fontsize=14)
tmp = plt.imshow(np.abs(np.squeeze(kspace_US[0,0,:,:])), aspect= 'auto')
plt.xlabel('ky');plt.ylabel('kx'); tmp.set_clim(0.0,0.01*np.max(np.abs(kspace_US))) # ky by kx
# Check which rows in the 4th dimension are entirely zeros
# axis=3 corresponds to the 4th dimension (rows)
zero_rows = np.all(us_mask_2d_plus_t == 0, axis=(0,1))
print(zero_rows.shape)
# Find rows that contain only zeros
kspace_US_collectedDataOnly = kspace_US[:,:,:,~zero_rows]
print(np.shape(kspace_US_collectedDataOnly))
# (18, 18, 512, 94))
fig1 = plt.figure(1, figsize=(8, 12)); fig1.suptitle("US K-space", fontsize=14)
tmp = plt.imshow(np.abs(np.squeeze(kspace_US_collectedDataOnly[0,0,:,:])), aspect= 'auto')
plt.xlabel('ky');plt.ylabel('kx'); tmp.set_clim(0.0,0.01*np.max(np.abs(kspace_US_collectedDataOnly))) # ky by kx
(18, 18, 512, 208)
(208,)
(18, 18, 512, 114)
# Although tfmri.recon.grappa supports batches, it is quite likely to run out
# of memory. Therefore we perform GRAPPA one phase at a time.
Final_cine = []
for t in range (np.shape(kspace_US_collectedDataOnly)[0]) :
#print('k in shape:', (np.squeeze(kspace_US_collectedDataOnly[t,:,:,:])).shape, ' , mask shape: ', ((np.squeeze(us_mask_2d_plus_t[t,:,:])).astype(np.bool)).shape, ' , calib_data shape: ', calib_data[t,:,:,:].shape)
# k in shape: (18, 512, 208) , mask shape: (512, 208) , calib_data shape: (18, 512, 208)
GRAPPAreconstructedImage = tfmri.recon.grappa((np.squeeze(kspace_US_collectedDataOnly[t,:,:,:])).astype(np.complex64), (np.squeeze(us_mask_2d_plus_t[t,:,:])).astype(np.bool), np.squeeze(calib_data[t,:,:,:]),
return_kspace=False)
Final_cine.append(GRAPPAreconstructedImage)
final_cine = np.array(Reconstructed_Image)
k in shape: (18, 512, 114) , mask shape: (512, 208) , calib_data shape: (18, 512, 20)
k : complex128 , mask : int64 , calib_data : <dtype: 'complex64'>
/tmp/ipykernel_230766/2446214914.py:7: DeprecationWarning: `np.bool` is a deprecated alias for the builtin `bool`. To silence this warning, use `bool` by itself. Doing this will not modify any behavior and is safe. If you specifically wanted the numpy scalar type, use `np.bool_` here.
Deprecated in NumPy 1.20; for more details and guidance: https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations
print('k in shape:', (np.squeeze(kspace_US_collectedDataOnly[t,:,:,:])).shape, ' , mask shape: ', ((np.squeeze(us_mask_2d_plus_t[t,:,:])).astype(np.bool)).shape, ' , calib_data shape: ', calib_data[t,:,:,:].shape)
/tmp/ipykernel_230766/2446214914.py:11: DeprecationWarning: `np.bool` is a deprecated alias for the builtin `bool`. To silence this warning, use `bool` by itself. Doing this will not modify any behavior and is safe. If you specifically wanted the numpy scalar type, use `np.bool_` here.
Deprecated in NumPy 1.20; for more details and guidance: https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations
GRAPPAreconstructedImage = tfmri.recon.grappa((np.squeeze(kspace_US_collectedDataOnly[t,:,:,:])).astype(np.complex64), (np.squeeze(us_mask_2d_plus_t[t,:,:])).astype(np.bool), np.squeeze(calib_data[t,:,:,:]),
k in shape: (18, 512, 114) , mask shape: (512, 208) , calib_data shape: (18, 512, 20)
k : complex128 , mask : int64 , calib_data : <dtype: 'complex64'>
k in shape: (18, 512, 114) , mask shape: (512, 208) , calib_data shape: (18, 512, 20)
k : complex128 , mask : int64 , calib_data : <dtype: 'complex64'>
k in shape: (18, 512, 114) , mask shape: (512, 208) , calib_data shape: (18, 512, 20)
k : complex128 , mask : int64 , calib_data : <dtype: 'complex64'>
k in shape: (18, 512, 114) , mask shape: (512, 208) , calib_data shape: (18, 512, 20)
k : complex128 , mask : int64 , calib_data : <dtype: 'complex64'>
k in shape: (18, 512, 114) , mask shape: (512, 208) , calib_data shape: (18, 512, 20)
k : complex128 , mask : int64 , calib_data : <dtype: 'complex64'>
k in shape: (18, 512, 114) , mask shape: (512, 208) , calib_data shape: (18, 512, 20)
k : complex128 , mask : int64 , calib_data : <dtype: 'complex64'>
k in shape: (18, 512, 114) , mask shape: (512, 208) , calib_data shape: (18, 512, 20)
k : complex128 , mask : int64 , calib_data : <dtype: 'complex64'>
k in shape: (18, 512, 114) , mask shape: (512, 208) , calib_data shape: (18, 512, 20)
k : complex128 , mask : int64 , calib_data : <dtype: 'complex64'>
k in shape: (18, 512, 114) , mask shape: (512, 208) , calib_data shape: (18, 512, 20)
k : complex128 , mask : int64 , calib_data : <dtype: 'complex64'>
k in shape: (18, 512, 114) , mask shape: (512, 208) , calib_data shape: (18, 512, 20)
k : complex128 , mask : int64 , calib_data : <dtype: 'complex64'>
k in shape: (18, 512, 114) , mask shape: (512, 208) , calib_data shape: (18, 512, 20)
k : complex128 , mask : int64 , calib_data : <dtype: 'complex64'>
k in shape: (18, 512, 114) , mask shape: (512, 208) , calib_data shape: (18, 512, 20)
k : complex128 , mask : int64 , calib_data : <dtype: 'complex64'>
k in shape: (18, 512, 114) , mask shape: (512, 208) , calib_data shape: (18, 512, 20)
k : complex128 , mask : int64 , calib_data : <dtype: 'complex64'>
k in shape: (18, 512, 114) , mask shape: (512, 208) , calib_data shape: (18, 512, 20)
k : complex128 , mask : int64 , calib_data : <dtype: 'complex64'>
k in shape: (18, 512, 114) , mask shape: (512, 208) , calib_data shape: (18, 512, 20)
k : complex128 , mask : int64 , calib_data : <dtype: 'complex64'>
k in shape: (18, 512, 114) , mask shape: (512, 208) , calib_data shape: (18, 512, 20)
k : complex128 , mask : int64 , calib_data : <dtype: 'complex64'>
k in shape: (18, 512, 114) , mask shape: (512, 208) , calib_data shape: (18, 512, 20)
k : complex128 , mask : int64 , calib_data : <dtype: 'complex64'>
print(final_cine.shape)
#remove RO oversampling
ROsize = final_cine.shape[2]
final_cine = final_cine[:,:,int(ROsize/4):int(3*ROsize/4),:]
print(final_cine.shape)
#[ph,ch, x, y ]
#(18, 18, 256, 208)
(18, 18, 512, 208)
(18, 18, 256, 208)
# now do coil combination
coil_combined_cine = tfmri.coils.combine_coils(final_cine, maps=None, coil_axis= 1)
print(coil_combined_cine.shape)
# 18, 256, 256
(18, 256, 208)
Finally lets see the final GRAPPA reconstructed image
import matplotlib.animation
plt.rcParams["animation.html"] = "jshtml"
plt.rcParams['figure.dpi'] = 150
plt.ioff()
fig, ax = plt.subplots()
nFrames = len(Reconstructed_Image)
t= np.linspace(0,len(Reconstructed_Image))
def animate(t):
plt.imshow(np.squeeze(tf.math.abs(coil_combined_cine[t,:,:])), cmap = 'gray')
plt.title('Image')
matplotlib.animation.FuncAnimation(fig, animate, frames=nFrames)
Conclusion¶
Congratulations! You performed a Cartesian GRAPPA reconstruction using TensorFlow MRI. The code used in this notebook works for any amount of undersampling/ACS lines. It also works for 3D imaging. Feel free to try with your own data!
For more information about the functions used in this tutorial, check out the API documentation. For more examples of using TensorFlow MRI, check out the tutorials.
Let us know!¶
Please tell us what you think about this tutorial and about TensorFlow MRI. We would like to hear what you liked and how we can improve. You will find us on GitHub.
# Copyright 2022 University College London. All rights reserved.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# https://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.