CG-SENSE
This tutorial is based on the RRSG_Challenge_01 github repository.
pip install
The cgsense2023 package was uploaded to PyPI and can be easily installed using the below command.
pip install cgsense2023
How to use
Plot Module
# !pip install ../.from cgsense2023.plot import *
import numpy as np# Sequential - full - spoke radial trajectory
traj_full = gen_radial_traj(golden=False, full_spoke=True)
show_trajectory(traj_full, golden=True)
# golden - full - spoke radial trajectory
traj_full_golden = gen_radial_traj(golden=True, full_spoke=True)
show_trajectory(traj_full_golden, golden=True)
Nim, Nx, Ny = 12, 128, 256
x1 = np.random.randn(Nim, Nx, Ny)
x2 = np.random.randn(1, Nx, Ny)show_image_grid(x1)Warning: number of images (12) is larger than number of panels (2x5)!
show_image_grid(x2)
Metrics Module
from cgsense2023.metrics import *
import numpy as np# same dimensions
Nim, Nx, Ny = 11, 128, 256
x1 = np.random.randn(1, Nx, Ny)
x2 = np.random.randn(1, Nx, Ny)
print_metrics(x1, x2)+---------+-----------+
| Metrics | Values |
+---------+-----------+
| MSE | 2.002e+00 |
| NMSE | 1.998e+00 |
| RMSE | 1.415e+00 |
| NRMSE | 1.414e+00 |
| PSNR | 15.966 |
| SSIM | 0.010322 |
+---------+-----------+
Gridding Reconstruction
from cgsense2023.math import *
from cgsense2023.io import *
import cgsense2023 as cgs2003
from cgsense2023.mri import *
from cgsense2023.optimizer import *
import h5py
import numpy as np
import matplotlib.pyplot as pltData Paths
# path to the data file
fpath = '../testdata/rawdata_brain.h5'
# path to the reference results
rpath = '../testdata/CG_reco_inscale_True_denscor_True_reduction_1.h5'with h5py.File(rpath, 'r') as rf:
print(list(rf.keys()))
ref_cg = np.squeeze(rf['CG_reco'][()][-1])
ref_grid = np.squeeze(rf['Coil_images'][()])['CG_reco', 'Coil_images']
ref_cg.shape, ref_grid.shape((300, 300), (12, 300, 300))
Setup Parameters
# one stop shop for all Parameters and Data
params = setup_params(fpath, R=1)['Coils', 'InScale', 'rawdata', 'trajectory']
Setup MRI Operator
# Set up MRI operator
mrimodel = MriImagingModel(params)
mriop = MriOperator(data_par=params["Data"],optimizer_par=params["Optimizer"])
mriop.set_operator(mrimodel)
# Single Coil images after FFT
my_grid = mriop.operator.NuFFT.adjoint(params['Data']['rawdata_density_cor'])my_grid.shape, ref_grid.shape((12, 300, 300), (12, 300, 300))
# test gridding recon results
np.allclose(ref_grid, my_grid)True
# test gridding recon results
np.array_equal(ref_grid, my_grid)False
print_metrics(np.abs(ref_grid[0]), np.abs(my_grid[0]))+---------+-----------+
| Metrics | Values |
+---------+-----------+
| MSE | 1.472e-24 |
| NMSE | 5.905e-15 |
| RMSE | 1.213e-12 |
| NRMSE | 7.684e-08 |
| PSNR | 154.87 |
| SSIM | 1.0 |
+---------+-----------+
show_image_grid(my_grid, figsize=(10,10), rows=3, cols=4)
Gradient Descent
guess = np.zeros((params['Data']['image_dim'],params['Data']['image_dim']))
SD_result, SD_residuals, SD_ref_res = steepest_descent(mriop, guess,
params['Data']['rawdata_density_cor'],
iters=50,
ref=ref_cg)Residuum at iter 50 : 6.553379e-06
show_compared_images(np.abs(ref_cg), np.abs(SD_result), diff_fac=10,
labels=['Reference', 'Steepest Descent', 'diff'])
np.allclose(ref_cg, SD_result)False
print_metrics(np.abs(ref_cg), np.abs(SD_result))+---------+-----------+
| Metrics | Values |
+---------+-----------+
| MSE | 5.633e-14 |
| NMSE | 7.888e-05 |
| RMSE | 2.373e-07 |
| NRMSE | 8.882e-03 |
| PSNR | 55.242 |
| SSIM | 0.9988 |
+---------+-----------+
CG’s Semi-Convergence Behavior
CG_result, CG_residuals, CG_ref_res = conjugate_gradient(mriop, guess,
params['Data']['rawdata_density_cor'],
iters=50,
ref=ref_cg)Residuum at iter 50 : 2.993229e-06
plt.plot(np.log10(SD_ref_res),'*--', label='SD reference_norms');
plt.plot(np.log10(SD_residuals),'*--', label='SD residual_norms');
plt.plot(np.log10(CG_ref_res),'*--', label='CG reference_norms');
plt.plot(np.log10(CG_residuals),'*--', label='CG residual_norms');
plt.grid();
plt.xlabel("# iteration")
plt.ylabel("residuals (log10)")
plt.legend();
Conjugate Gradient vs. REF
Based on the “semi-convergence” plot above, the optimal number of iterations for CG and SD are around 10 and 28, respectively.
CG_result_vs_REF, _, _ = conjugate_gradient(mriop, guess,
params['Data']['rawdata_density_cor'],
iters=10,
ref=ref_cg)Residuum at iter 10 : 1.826174e-05
show_compared_images(np.abs(ref_cg), np.abs(CG_result_vs_REF), diff_fac=200000,
labels=['Reference', 'Conjugate Gradient', 'diff'])
np.allclose(ref_cg, CG_result_vs_REF)True
print_metrics(np.abs(ref_cg), np.abs(CG_result_vs_REF))+---------+-----------+
| Metrics | Values |
+---------+-----------+
| MSE | 1.196e-22 |
| NMSE | 1.675e-13 |
| RMSE | 1.094e-11 |
| NRMSE | 4.093e-07 |
| PSNR | 141.97 |
| SSIM | 1.0 |
+---------+-----------+
Developer install
If you want to develop cgsense2023 yourself, please use an editable
installation of cgsense2023.
git clone https://github.com/hdocmsu/cgsense2023.git
pip install -e "cgsense2023[dev]"
You also need to use an editable installation of nbdev, fastcore, and execnb.
Happy Coding!!!
References
This template was created based on the below references. The list is not exhaustive so if you notice any missing references, please report an issue. I will promptly correct it.
-
RRSG_Challenge_01 github repository
-
fastMRI github repository
-
CG-SENSE revisited: Results from the first ISMRM reproducibility challenge
-
Prof. James V. Burke’s Lecture on “The Conjugate Gradient Algorithm”
-
Magnetic Resonance Image Reconstruction: Theory, Methods, and Applications, (2022)
Invitation to Contributions
If you would like to contribute to improve this repository please feel free to propose here.
License
Apache License 2.0
