Synopsis

We present SigPy, a Python package designed for high performance iterative reconstruction. Its main features include:

- A unified CPU and GPU Python interface to signal processing functions, including convolution, FFT, NUFFT, wavelet transform, and thresholding functions.

- Convenient classes (Linop, Prox, Alg, App) to build more complicated iterative reconstruction algorithms.

- Commonly used MRI reconstruction methods as Apps, including SENSE, L1-wavelet regularized reconstruction, total-variation regularized reconstruction, and JSENSE.

- MRI-specific functions, including poisson-disc sampling, ESPIRiT calibration, and non-Cartesian preconditioners.

- Simple installation via pip and conda.

Introduction

To develop high performance and clinically feasible iterative reconstruction methods, MRI researchers usually rely on self written C or C++ libraries or existing ones such as BART [1] and Gadgetron [2]. While C/C++ implementations are faster than implementations in higher level languages, they are more difficult to use and extend for most developers. This prolongs the turnaround time to translate research methods to clinical practice.

On the other hand, all major machine learning packages now use Python as their main interface, and have shown competitive runtimes. The user programs her/his high-level algorithm directly in Python, and the package automatically translates the code to fast CPU/GPU runtime implementations. Inspired by this, we are interested in developing a Python package for high performance iterative reconstruction that shares this design philosophy.

We present SigPy, a Python package designed for high performance iterative reconstruction. Its main features include:

- A unified CPU and GPU Python interface to signal processing functions, including convolution, FFT, NUFFT, wavelet transform, and thresholding functions.

- Convenient classes (Linop, Prox, Alg, App) to build more complicated iterative reconstruction algorithms.

- Commonly used MRI reconstruction methods as Apps, including SENSE, L1-wavelet regularized reconstruction, total-variation regularized reconstruction, and JSENSE.

- MRI-specific functions, including poisson-disc sampling, ESPIRiT calibration, and non-Cartesian preconditioners.

- Simple installation via pip and conda.

Resources

SigPy is available through pip and conda, and can be installed with the following commands:

$ pip install sigpy

$ conda install -c frankong sigpy

The source code is hosted on Github with BSD 3-Clause License:

https://github.com/mikgroup/sigpy

And the documentation is hosted on:

http://sigpy.readthedocs.io

Basic Usage

SigPy is designed to have as little learning curve as possible. Since almost all Python users already use NumPy, SigPy operates on NumPy [3] arrays directly on CPU, and avoids defining any redundant functions. When NumPy implementation is slow, SigPy uses Numba [4] instead to translate Python functions to optimized machine code at runtime. For example, gridding functions in SigPy are implemented using Numba. For GPU, SigPy operates on CuPy arrays [5] , which have the same interface as NumPy but are implemented in CUDA.

Computing devices can easily be switched on the fly in SigPy. In particular, it uses the same context switching mechanism used in TensorFlow, PyTorch, and CuPy as shown in Figure 1.

Building iterative reconstruction methods

SigPy provides four abstraction classes (Linop, Prox, Alg, and App) to help users building iterative methods. Such abstraction is inspired by similar structure in BART [1].

The Linop class abstracts a linear operator, and supports adjoint, addition, composing, and stacking. For example, given a Linop A, A.H * A defines the normal operator. Prepackaged Linops include FFT, NUFFT, and wavelet, and common array manipulation functions.

The Prox class abstracts a proximal operator, and can do stacking and conjugation. Prepackaged Proxs include L1/L2 regularization and projection functions.

The Alg class abstracts iterative algorithms. Prepackaged Algs include conjugate gradient, (accelerated/proximal) gradient method, and primal dual hybrid gradient.

Finally, the App class wraps the above three classes into a final deliverable application. Users can run an App without knowing the internal implementation. Prepackaged MRI Apps include SENSE reconstruction, l1 wavelet reconstruction, total variation reconstruction, and JSENSE reconstruction. The relationship between the classes is shown in Figure 2.

These classes together allow users to quickly build iterative methods. Figure 3 shows an example of building and running an l1-wavelet reconstruction App using 12 lines of Python code.

Performance comparison with BART

Here we briefly compare the reconstruction speed of SigPy versus BART. Being a Python library, SigPy should not be expected to be faster than BART, a C library. However, for reconstructions that are dominated by FFTs, SigPy can often closely match BART in terms of speed because they all rely on the same FFT libraries.

In particular, we ran a GPU Cartesian SENSE reconstruction of image size 256x256x256 with 8-channel coil array, solved with 30 iterations of the conjugate gradient method. BART took 15.94s, and SigPy took 16.50s.

Conclusion

SigPy provides a simple Python interface to build iterative reconstruction methods on CPU and GPU. Figure 4 shows an example GPU 3D cones reconstruction result implemented only using SigPy. Because of SigPy’s pure Python interface, complicated reconstruction methods can be easily implemented.

Similar to Tensorflow and Pytorch, the ability to use GPUs allows SigPy to achieve fast reconstruction times despite being implemented in a high-level language.

Finally, the Python implementation allows users to use SigPy along with other Python packages, such as Tensorflow and Pytorch, without switching between different programming languages.

Acknowledgements

References

[1] Martin Uecker, Frank Ong, Jonathan I Tamir, Dara Bahri, Patrick Virtue, Joseph Y Cheng, Tao Zhang, and Michael Lustig, Berkeley Advanced Reconstruction Toolbox, Annual Meeting ISMRM, Toronto 2015, In Proc. Intl. Soc. Mag. Reson. Med. 23:2486

[2] Hansen, M.S. and Sørensen, T.S., 2013. Gadgetron: an open source framework for medical image reconstruction. Magnetic resonance in medicine, 69(6), pp.1768-1776.

[3] Walt, S.V.D., Colbert, S.C. and Varoquaux, G., 2011. The NumPy array: a structure for efficient numerical computation. Computing in Science & Engineering, 13(2), pp.22-30.

[4] Lam, S.K., Pitrou, A. and Seibert, S., 2015, November. Numba: A LLVM-based python JIT compiler. In Proceedings of the Second Workshop on the LLVM Compiler Infrastructure in HPC (p. 7). ACM.

[5] Ryosuke Okuta, Yuya Unno, Daisuke Nishino, Shohei Hido and Crissman Loomis. CuPy: A NumPy-Compatible Library for NVIDIA GPU Calculations. Proceedings of Workshop on Machine Learning Systems (LearningSys) in The Thirty-first Annual Conference on Neural Information Processing Systems (NIPS), (2017)