Like most fields of science, in magnetic resonance imaging (MRI) our appetite for data is never satiated – spatiotemporal resolution and signal-to-noise ratio can never be too high, and MRI scan times can never be too low. But, with big data comes big compute. MRI image reconstruction, for example, often involves parallel acquisition of multiple (3+1)D MRI images as measured by each of the 32+ scanner readout coils which are then combined through a complex series of iterative optimization problems. These types of algorithms are typically prototyped in high-level programming languages like Matlab or Python, and subsequently translated to C/C++ for deployment. The Julia community is familiar with this old story, though – in fact the fantastic package MRIReco.jl (https://github.com/MagneticResonanceImaging/MRIReco.jl) by Knopp et al. provides Julia implementations of MRI reconstruction algorithms which are competitive with C/C++. Post-processing these reconstructed MRI images faces similar computational challenges, and in this talk, we will describe our experience of implementing a parameter inference algorithm from the MRI subfield of myelin water imaging (MWI) in Julia.

In MWI, one analyses (3+1)D time series of image volumes acquired on cartesian spatiotemporal grids with dimensions 250x250x250x64 or more. The MRI time signals acquired in each voxel, which exhibit multi-exponential decay, are decomposed into a spectrum of decay rates. These decay rates are used to compute, among other useful metrics, the myelin water fraction (MWF) which is known to correlate with local myelin content in the brain. This inverse Laplace transform-like computation involves fitting an MRI signal model – the extended phase graph (EPG) model – to each time signal by solving an L2 regularized nonnegative least squares (NNLS) optimization problem. Note that, with images typically consisting of 10^7 voxels or more, this computation requires solving upwards of 10^7 optimization problems.

Prior to this work, the NNLS procedure used for MWI was implemented in Matlab – a closed-source high-level programming language – as is common in MRI research. This computation is particularly poorly suited for Matlab, however. First, similar to the Python library numpy, Matlab encourages computations on vectors or matrices, as opposed to explicit for-loops, in order to call out to BLAS or LAPACK libraries and ameliorate the overhead of the Matlab interpreter. For this reason, the EPG algorithm – which is most efficiently expressed in terms of nested for-loops – was previously implemented in terms of sparse matrix-vector products in order to avoid Matlab’s slow loops. Second, while the solving of independent optimization problems from each voxel is embarrassingly parallel, Matlab provides little control over multiprocessing optimizations such as reusing thread-local memory buffers and task scheduling. Lastly, Matlab does not provide a statically sized array type, which would be beneficial for micro-optimizing 3x3 matrix-vector products present in the EPG algorithm.

Julia excels in these types of computations. In the DEcomposition and Component Analysis of Exponential Signals (DECAES.jl) package (https://github.com/jondeuce/DECAES.jl, https://doi.org/10.1016/j.zemedi.2020.04.001), we provide optimized procedures for computing MWI which address the aforementioned limitations of Matlab, and additionally include command line and Matlab interfaces for ease of interoperability. In all, DECAES.jl reduced computation times approximately 60X from 1.5-2.5 hours down to 1.5-2.5 mins. This large speedup demonstrates that it is possible to perform this analysis directly on the MRI scanner, removing the need for researchers to (manually) process the acquired data.

Among the many additional benefits from the Julia translation is the synergy with other Julia packages: we experimented with explicit SIMD vectorization in the EPG algorithm using the SIMD.jl package; we make liberal use of statically sized vectors and matrices using the StaticArrays.jl package; we use a pure-Julia implementation of NNLS using the NNLS.jl package. Furthermore, the EPG algorithm is independently useful outside of MWI, and can e.g. be efficiently differentiated trivially using the automatic differentiation packages ForwardDiff.jl or Zygote.jl.

In conclusion, we have found that the combination of high-performance and high-expressibility present in Julia is well suited to MRI research, particularly in comparison to existing Matlab-based workflows, and we believe that our experience will resonate with the scientific computing community more broadly. We look forward to the opportunity to share our experience.