Metadata-Version: 2.1
Name: cgsense2023
Version: 0.0.2
Summary: Conjugate Gradient SENSE (CG-SENSE) MRI Reconstruction
Home-page: https://github.com/hdocmsu/cgsense2023/
Author: Hung Do, PhD
Author-email: clinicalcollaborations@gmail.com
License: MIT License
Keywords: nbdev
Classifier: Development Status :: 3 - Alpha
Classifier: Intended Audience :: Developers
Classifier: Natural Language :: English
Classifier: Programming Language :: Python :: 3.7
Classifier: Programming Language :: Python :: 3.8
Classifier: Programming Language :: Python :: 3.9
Classifier: Programming Language :: Python :: 3.10
Classifier: License :: OSI Approved :: MIT License
Requires-Python: >=3.7
Description-Content-Type: text/markdown
Provides-Extra: dev
License-File: LICENSE

# CG-SENSE HOME

<!-- WARNING: THIS FILE WAS AUTOGENERATED! DO NOT EDIT! -->

This tutorial is based on the
[RRSG_Challenge_01](https://github.com/ISMRM/rrsg_challenge_01) github
repository.

## pip install

The [cgsense2023](https://pypi.org/project/cgsense2023/) package was
uploaded to [PyPI](https://pypi.org/) and can be easily installed using
the below command.

`pip install cgsense2023`

## How to use

### Plot Module

``` python
# !pip install ../.
```

``` python
from cgsense2023.plot import *
import numpy as np
```

``` python
# Sequential - full - spoke radial trajectory
traj_full = gen_radial_traj(golden=False, full_spoke=True)
show_trajectory(traj_full, golden=True)
```

![](index_files/figure-commonmark/cell-4-output-1.png)

``` python
# golden - full - spoke radial trajectory
traj_full_golden = gen_radial_traj(golden=True, full_spoke=True)
show_trajectory(traj_full_golden, golden=True)
```

![](index_files/figure-commonmark/cell-5-output-1.png)

``` python
Nim, Nx, Ny = 12, 128, 256
x1 = np.random.randn(Nim, Nx, Ny)
x2 = np.random.randn(1, Nx, Ny)
```

``` python
show_image_grid(x1)
```

    Warning: number of images (12) is larger than number of panels (2x5)!

![](index_files/figure-commonmark/cell-7-output-2.png)

``` python
show_image_grid(x2)
```

![](index_files/figure-commonmark/cell-8-output-1.png)

### Metrics Module

``` python
from cgsense2023.metrics import *
import numpy as np
```

``` python
# 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

``` python
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

``` python
# 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'
```

``` python
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']

``` python
ref_cg.shape, ref_grid.shape
```

    ((300, 300), (12, 300, 300))

#### Setup Parameters

``` python
# one stop shop for all Parameters and Data
params = setup_params(fpath, R=1)
```

    ['Coils', 'InScale', 'rawdata', 'trajectory']

#### Setup MRI Operator

``` python
# 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'])
```

``` python
my_grid.shape, ref_grid.shape
```

    ((12, 300, 300), (12, 300, 300))

``` python
# test gridding recon results
np.allclose(ref_grid, my_grid)
```

    True

``` python
# test gridding recon results
np.array_equal(ref_grid, my_grid)
```

    False

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

``` python
show_image_grid(my_grid, figsize=(10,10), rows=3, cols=4)
```

![](index_files/figure-commonmark/cell-21-output-1.png)

### Gradient Descent

``` python
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

``` python
show_compared_images(np.abs(ref_cg), np.abs(SD_result), diff_fac=10, 
                     labels=['Reference', 'Steepest Descent', 'diff'])
```

![](index_files/figure-commonmark/cell-23-output-1.png)

``` python
np.allclose(ref_cg, SD_result)
```

    False

``` python
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

``` python
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

``` python
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();
```

![](index_files/figure-commonmark/cell-27-output-1.png)

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

``` python
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

``` python
show_compared_images(np.abs(ref_cg), np.abs(CG_result_vs_REF), diff_fac=200000, 
                     labels=['Reference', 'Conjugate Gradient', 'diff'])
```

![](index_files/figure-commonmark/cell-29-output-1.png)

``` python
np.allclose(ref_cg, CG_result_vs_REF)
```

    True

``` python
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](https://github.com/fastai/nbdev),
[fastcore](https://github.com/fastai/fastcore), and
[execnb](https://github.com/fastai/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](https://github.com/hdocmsu/cgsense2023/issues/new). I will
promptly correct it.

- [RRSG_Challenge_01](https://github.com/ISMRM/rrsg_challenge_01) github
  repository

- [fastMRI](https://github.com/facebookresearch/fastMRI) github
  repository

- [Original CG-SENSE Paper by Pruessmann et. al.,
  2001](https://onlinelibrary.wiley.com/doi/10.1002/mrm.1241)

- [CG-SENSE revisited: Results from the first ISMRM reproducibility
  challenge](https://onlinelibrary.wiley.com/doi/10.1002/mrm.28569)

- [Dwight Nishimura, (Published Jan 10, 2010), “Principles of Magnetic
  Resonance Imaging,” Stanford University, Palo Alto,
  CA](https://www.lulu.com/shop/dwight-nishimura/principles-of-magnetic-resonance-imaging/hardcover/product-1gvnv7yn.html?page=1&pageSize=4)

- [Prof. James V. Burke’s Lecture on “The Conjugate Gradient
  Algorithm”](https://sites.math.washington.edu/~burke/crs/408/lectures/L14-CG.pdf)

- [Magnetic Resonance Image Reconstruction: Theory, Methods, and
  Applications,
  (2022)](https://shop.elsevier.com/books/magnetic-resonance-image-reconstruction/akcakaya/978-0-12-822726-8)

- [sigPy github repo](https://github.com/mikgroup/sigpy)

# Invitation to Contributions

If you would like to contribute to improve this repository please feel
free [to propose
here](https://github.com/hdocmsu/cgsense2023/issues/new).

# License

MIT License
