cgsense2023 1.0.2

Conjugate Gradient SENSE (CG-SENSE) MRI Reconstruction

CG-SENSE HOME

Step-by-step tutorial of Conjugate Gradient SENSE (CG-SENSE) reconstruction by Hung P. Do.

Citing this work

DOI

Do, H. P. (2026). CG-SENSE Tutorial version 1.0.0 (1.0.0). Zenodo. https://doi.org/10.5281/zenodo.18526383

Organization of the tutorial

This tutorial was written by Hung Do based on the RRSG_Challenge_01 github repository.

• Original tutorial notebook

• Original python code:

  • helper_fun.py
  • linop.py
  • solver.py
  • recon.py

• Hung Do’s modifications:

  • Modified tutorial notebook
  • Development notebook for the cgsense2023 PyPI package

• cgsense2023 PyPI Package re-written: Details can be seen in the P4_HPD_cgsense section.

My modifications are:

1. Most importantly, I organized the code into coherent modules and repackaged the code into a Python package named cgsense2023 for easy installation and use. You can install the package by running pip install cgsense2023 in your terminal.

2. For teaching purposes, I created this tutorial for anyone who is interested in learning about CG-SENSE and the artifacts that can arise from it.

3. The gridding part was untouched, but I retouched the I/O interface such that parameters are consolidated and grouped in one and only one central location.

4. For my own learning purposes, I re-implemented the CG-SENSE algorithm from scratch without referring to the original code. There are several variations of the CG-SENSE algorithm, but the results are generally similar.

5. Implemented gradient descent (GD) and steepest descent (SD) algorithms for comparison with CG-SENSE.

6. Added code to create k-space trajectories and visualize them.

7. Added quantitative metrics for evaluating the reconstructed images, such as RMSE and SSIM, etc.

8. Added code to visualize the artifacts in the reconstructed images and assess convergence of the CG-SENSE algorithm.

W.I.P Tasks

Contributions to the below tasks are welcome!

• Document the standardization of the raw h5-file, otherwise read_data() has to be modified every time a deviation in h5-file is encountered.

• Refactor, simplify, modularize the gridding reconstructions. Currently, it is difficult to understand, especially for new learner.

• Add dstore flag in gradient descent or conjugate gradient optimization functions to store intermediate reconstruction in each iteration.

• Add theoretical descriptions on MRI, MRI recons, and gradient Descent, and gonjugate gradient algorithms.

• etc.

How it works

pip install

The cgsense2023 package was uploaded to PyPI and can be easily installed using the below command.

pip install cgsense2023

Plot Module

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, figsize=(10,8))

# golden - full - spoke radial trajectory
traj_full_golden = gen_radial_traj(golden=True, full_spoke=True)
show_trajectory(traj_full_golden, golden=True, figsize=(10,8))

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   | 1.999e+00 |
|   NMSE  | 2.000e+00 |
|   RMSE  | 1.414e+00 |
|  NRMSE  | 1.414e+00 |
|   PSNR  |   14.981  |
|   SSIM  | 0.0092604 |
+---------+-----------+

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 plt

Data 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)

show_image_grid(rss_rec(my_grid), figsize=(10,10))

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    |
+---------+-----------+

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

• Original CG-SENSE Paper by Pruessmann et. al., 2001

• CG-SENSE revisited: Results from the first ISMRM reproducibility challenge

• Dwight Nishimura, (Published Jan 10, 2010), “Principles of Magnetic Resonance Imaging,” Stanford University, Palo Alto, CA

• Prof. James V. Burke’s Lecture on “The Conjugate Gradient Algorithm”

• Magnetic Resonance Image Reconstruction: Theory, Methods, and Applications, (2022)

• sigPy github repo

Citing this work

DOI

Do, H. P. (2026). CG-SENSE Tutorial version 1.0.0 (1.0.0). Zenodo. https://doi.org/10.5281/zenodo.18526383