PyQMRI: An accelerated Python based Quantitative MRI toolbox

Various medical examinations are seeing a shift to a more patient centric and personalized view, based on quantitative instead of qualitative observations and comparisons. This trend has also affected medical imaging, and particularly quantitative MRI (qMRI) gained importance in recent years. qMRI aims to identify the underlying biophysical and tissue parameters that determine contrast in an MR imaging experiment. In addition to contrast information, qMRI permits insights into diseases by providing biophysical, microstructural, and functional information in absolute quantitative values. For quantification, biophysical models are used, which describe the relationship between image intensity and physical properties of the tissue for certain scanning sequences and sequence parameters. By performing several measurements with different sequence parameters (e.g. flip angle, repetition time, echo time) the related inverse problem of identifying the tissue parameters sought can be solved.


Summary
Various medical examinations are seeing a shift to a more patient centric and personalized view, based on quantitative instead of qualitative observations and comparisons. This trend has also affected medical imaging, and particularly quantitative MRI (qMRI) gained importance in recent years. qMRI aims to identify the underlying biophysical and tissue parameters that determine contrast in an MR imaging experiment. In addition to contrast information, qMRI permits insights into diseases by providing biophysical, microstructural, and functional information in absolute quantitative values. For quantification, biophysical models are used, which describe the relationship between image intensity and physical properties of the tissue for certain scanning sequences and sequence parameters. By performing several measurements with different sequence parameters (e.g. flip angle, repetition time, echo time) the related inverse problem of identifying the tissue parameters sought can be solved.
Quantitative MR typically suffers from increased measurement time due to repeated imaging experiments. Therefore, methods to reduce scanning time by means of optimal scanning protocols and subsampled data acquisition have been extensively studied. However, these approaches are typically associated with a reduced SNR, and can suffer from subsampling artifacts. To address both aspects, it has been shown that the inclusion of a biophysical model in the reconstruction process leads to much faster data acquisition, while simultaneously improving image quality. The inverse problem associated with this special reconstruction approach requires dedicated numerical solution strategies (Block et al., 2009;Doneva et al., 2010;Donoho, 2006;Lustig et al., 2007;Maier, Schoormans, et al., 2019;Roeloffs et al., 2016;Sumpf, 2012), commonly known as model-based reconstruction. Model-based reconstruction is based on variational modeling, and combines parallel imaging and compressed sensing to achieve acceleration factors as high as 10x the speed of fully sampled acquisitions. The method directly solves for the unknown parameter maps from raw k-space data. The repeated transition from k-space to image-space, combined with the involved non-linear iterative reconstruction techniques to identify the unknown parameters, often leads to prolonged reconstruction times. This effect gets even more demanding if 3D image volumes are of interest.
Recently, the upsurge of computationally powerful GPUs has led to a variety of GPU based implementations to speed up computation time of highly parallelizeable operations (e.g., the Fourier transformation in MRI (Knoll et al., 2014)). As model-based approaches possibly deal with hundred Gigabytes of data (e.g. diffusion tensor imaging), available memory on current GPUs (e.g. 12 GB) can be a limiting factor. Thus, most reconstruction and fitting algorithms are applied in a slice-by-slice fashion to the volumetric data by taking a Fourier transformation along a fully sampled acquisition direction, effectively yielding a set of 2D problems. Hence, the additional information in form of the third dimension of volumetric data is neglected, leading to a loss in performance.
To utilize full 3D information in advanced reconstruction and fitting algorithms on memory limited GPUs, special solutions strategies are necessary to leverage the speed advantage, e.g., hiding memory latency of repeated transfers to/from the GPU to host memory. This can be achieved using asynchronous execution strategies. However, correct synchronization of critical operations can be error prone. To this end, we propose PyQMRI, a simple to use Python toolbox for quantitative MRI.

Statement of need
PyQMRI aims at reducing the required reconstruction time by means of a highly parallelized PyOpenCL (Klöckner et al., 2012) implementation of a state-of-the-art model-based reconstruction and fitting algorithm, while maintaining the easy-to-use properties of a Python package. In addition to processing small data (e.g. 2D slices) completely on the GPU, an efficient double-buffering based solution strategy is implemented. Double-buffering allows overlapping computation and memory transfer from/to the GPU, thus hiding the associated memory latency. By overlapping the transferred blocks it is possible to pass on 3D information utilizing finite differences based regularization strategies . Figure 1 shows a schematic of the employed double-buffering scheme. To make sure that this asynchronous execution strategy yields the expected results, unit-testing is employed. Currently, 3D acquisitions with at least one fully sampled dimension can be reconstructed on the GPU, including stack-of-X acquisitions or 3D Cartesian based imaging. Of course 2D data can be reconstructed as well. The combination of reconstruction and non-linear fitting is based on an iteratively regularized Gauss-Newton (IRGN) approach combined with a primal-dual inner loop. Regularization strategies include total variation (TV) (Rudin et al., 1992) and total generalized variation (TGV) (Bredies et al., 2010;Knoll et al., 2011) using finite differences gradient operations. In addition to the combined reconstruction and fitting algorithm from k-space data, PyQMRI can also be used to speed-up non-linear parameter fitting of complex or real valued image data. The main advantage of fitting the complex (k-space) data is that the assumed Gaussian noise characteristics for the commonly used L 2 data fidelity term are valid. This is especially important for problems suffering from poor SNR, e.g. Diffusion Tensor Imaging, where the wrong noise assumption can lead to significant errors in the quantification process (Jones & Basser, 2004).
PyQMRI comes with several pre-implemented quantiative models. In addition, new models can be introduced via a simple text file, utilizing the power of SymPy to generate numerical models as well as their partial derivatives in Python. Fitting can be initiated via a command line interface (CLI) or by importing the package into a Python script. Due to PyQMRI's OpenCL backend, no vendor specific hardware restrictions are present. However, current limitations of the gpyfft package, used to wrap the clfft, constrain the use to GPU devices only. A switch to other clfft wrappers might solve this limitation in future releases, but gpyfft is the only one that currently supports fast non-power-of-two transformations up to 13.

Related Work
The increased importance of qMRI is reflected by a multitude of open-source toolboxes, each focusing on a subset or combination of qMRI applications. Most tools show a strong focus on neurological applications (mostly brain) but are usually not limited to this application area. hMRI (Tabelow et al., 2019) is Matlab based and builds upon SPM (Friston, 2007). It extends the spatial registration and statistical inference capabilities of SPM by relaxometry and quantification of the magnetisation transfer effect. mrQ (Mezer et al., 2013) offers relaxometry combined with the ability to quantify the macromolecular tissue volume, apparent volume of interacting water protons, and the water-surface interaction rate and is completely written in Matlab. Another Matlab based project is qMRlab (Karakuzu et al., 2020) which offers a multitude of quantification algorithms including relaxometry, diffusion imaging, quantitative susceptibility mapping, field mapping, and quantitative magnetization transfer. It further offers routines for visualization, simulation, and protocol optimization of quantitative MRI examinations. Another Matlab based software is qmap (Hurley et al., 2011) which offers a collection of tools for quantitative MRI.
qMRI is also present in the Python community with PyMRT (Metere & Möller, 2017) offering tools for image analysis and relaxometry. Another powerful Pyhton package, with a focus on neuroimaing, is DIPY (Garyfallidis et al., 2014), offering a multitude of ways to evaluated diffusion and perfusion MRI data. Other software packages focus on fast execution and fitting, like the QUIT (Wood, 2018) toolbox, which is entirely written in C++ to speed up the computations. All of the above mentioned qMRI toolboxes have in common that they usually require image data for the fitting process and, thus, are not suitable for accelerated acquired data or require dedicated reconstruction algorithms prior to fitting.
A recent extensions to BART (Uecker et al., 2015) allows for T 1 quantification from 2D radially acquired inversion recovery Look-Locker data. The approach utilizes a model-based reconstruction algorithm to estimate T 1 directly from k-space (Wang et al., 2018(Wang et al., , 2019. Even though this approach can handle undersampled data and incorporates the whole MRI acquisition pipeline, it is currently limited to this single quantification model.
To the best of the authors knowledge PyQMRI is the only available Python toolbox that offers real 3D regularization in an iterative solver for model-based qMRI problems and for arbitrary large volumetric data, while simultaneously utilizing the computation power of recent GPUs. Further, the ability to use symbolic equations to generate new models seems to be unique as other tools require modifications of the code to include new quantification models.

Algorithms
PyQMRI deals with the following general problem structure: which includes a non-linear forward operator (A), mapping the parameters u to (complex) data space d, and a non-smooth regularization functional due to the L 1 -norms of the T(G)V functional (Bredies et al., 2010;Knoll et al., 2011). Setting α 1 = 0 and v = 0 the problem becomes simple TV regularization (Rudin et al., 1992). The gradient ∇ and symmetrized gradient E operators are implemented using finite differences. To further improve the quality of the reconstructed parameter maps PyQMRI uses a Frobenius norm to join spatial information from all maps in the T(G)V functionals (Bredies, 2014;Knoll et al., 2017;Maier, Schoormans, et al., 2019). Box constraints, limiting each unknown parameter in u to a physiological meaningful range, can be set in conjunction with real or complex value constraints.
Following the Gauss-Newton approach a sequence k of linearized sub-problems of the form needs to be solved to find a solution of the overall problem. The matrix DA| u=u k = ∂A ∂u (u k ) resembles the Jacobian of the system. The subproblems can be recast into a saddle-point structure by application of the Fenchel duality min u max y ⟨Ku, y⟩ + G(u) − F * (y), and solved utilizing a well established primal-dual algorithm (Chambolle & Pock, 2011) combined with a line-search (Malitsky & Pock, 2018) to speed-up convergence. Constant terms, stemming from the linearization, are precomputed and fused with the data d, yieldingd k . The inclusion of the additional L 2 -norm penalty improves convexity of the subproblem and resembles a Levenberg-Marquardt update for M k = diag(DA| T u=u k DA| u=u k ). A graphical representation of the involved steps is given in Figure 2. The regularization weights, regularization type (TV/TGV), and the number of outer and inner iterations can be changed using a plain text configuration file. It was shown by (Salzo & Villa, 2012) that the GN approach converges with linear rate to a critical point for non-convex problems with non-differential penalty functions if the initialization is sufficiently close. Thus a meaningful initial guess based on physiological knowledge on the parameters u should be used to initialize the fitting, e.g. mean T 1 value of the tissue of interest. Figure 2: Graphical representation of the employed regularized non-linear fitting procedure shown for an exemplary T1 quantification problem. Ci describes complex coil sensitivity information, F amounts to the sampling process including the Fourier transformation, and Sp equals the non-linear relationship between image intensity and the unknown physical quantities (T1 and Proton Density (PD)).