LEARNING THE DOMAIN SPECIFIC INVERSE NUFFT FOR ACCELERATED SPIRAL MRI

USING DIFFUSION MODELS

Trevor J. Chan1 , Chamith S. Rajapakse2

1. Department of Bioengineering, University of Pennsylvania, Philadelphia, United States

2. Department of Radiology, University of Pennsylvania, Philadelphia, United States
arXiv:2404.12361v1 [cs.AI] 18 Apr 2024

ABSTRACT ratio of undersampling, but comes at the cost of image qual-

ity. Sampling below the Nyquist limit introduces ambiguities
Deep learning methods for accelerated MRI achieve state-of- during reconstruction which manifest as artifacts in the final
the-art results but largely ignore additional speedups possible image. Reconstruction algorithms must therefore leverage ad-
with noncartesian sampling trajectories. To address this ditional information, including multicoil data in the case of
gap, we created a generative diffusion model-based recon- parallel imaging, and sparse priors, in the case of compressed
struction algorithm for multi-coil highly undersampled spiral sensing, in order to resolve these ambiguities [3, 4].
MRI. This model uses conditioning during training as well A third more recent approach to undersampled MRI re-
as frequency-based guidance to ensure consistency between construction lies in deep learning methods, which essentially
images and measurements. Evaluated on retrospective data, learn a set of image priors and use these to regularize so-
we show high quality (structural similarity >0.87) in re- lutions to the ill-posed reconstruction problem [5]. Within
constructed images with ultrafast scan times (0.02 seconds this category, diffusion models stand out for producing state-
for a 2D image). We use this algorithm to identify a set of of-the art results on image reconstruction tasks for faster
optimal variable-density spiral trajectories and show large scanning, motion correction, noise reduction, and others
improvements in image quality compared to conventional [6, 7, 8, 9]. Despite this, the vast majority of deep learn-
reconstruction using the non-uniform fast Fourier transform. ing approaches, and to our knowledge, all diffusion-based
By combining efficient spiral sampling trajectories, multicoil approaches to image reconstruction, focus on Cartesian-
imaging, and deep learning reconstruction, these methods sampled MRI, missing out on potential acceleration gains
could enable the extremely high acceleration factors needed attained by using more efficient non-Cartesian sampling tra-
for real-time 3D imaging. jectories. This work addresses this gap by introducing a
Index Terms— Accelerated MRI, Spiral MRI, Deep diffusion model-based method for trajectory-agnostic un-
Learning, Image Reconstruction dersampled MRI reconstruction of multicoil spiral MRI.

Contributions:
1. INTRODUCTION • Creation of a novel multi-conditioning strategy for solv-
ing the inverse non-uniform fast Fourier transform (nufft)
Despite the numerous advantages of MRI, inherent physical using a learned conditional score function and weak fre-
and hardware constraints cap acquisition speed and lead to quency guidance during sampling
long scan times. This creates myriad downstream obstacles, • Efficient hyperparameter search of the joint trajectory-
including low patient compliance, inefficient resource allo- reconstruction space and identification of optimal sam-
cation, and image motion artifacts, among others. For this pling trajectories
reason, methods to accelerate MRI acquisition have been and • Retrospective acquisition and reconstruction of a 2D,
continue to be an active area of research. Acceleration can be 256x256 pixel, 22x22 cm2 image with a readout duration
achieved through a combination of using more efficient scan- of 0.02 seconds.
ning sequences and reducing the number of measurements
made in the frequency space of the image, k-space. 2. BACKGROUND
Considering the former, successful techniques for faster
scanning include radial and spiral imaging methods, which Canonical score-based models consider the mapping between
exploit the unequal distribution of information across k-space a known distribution of independently and identically dis-
by sampling lower frequencies more densely [1, 2]. tributed samples of gaussian noise and an observed, but un-
Considering the latter, acquiring fewer measurments in k- known distribution of data p(x). These distributions book-
space rewards a decrease in scan time proportional the the end a Langevin diffusion process described by the following
 Fig. 2. Given measurements y0 , reconstruction follows a
Fig. 1. Example trajectories (A) and the corresponding read- modified diffusion sampling process. At each timestep, a
out gradients in kx and ky (B). All trajectories shown cover noisy latent xt is concatenated with a prior p0 and passed to
the frequency space of a 256x256 image and have a readout the denoising model to obtain x̃t−1 . To enforce consistency
duration of 10.0 ms. with y0 , we compute a frequency gradient ∇yt−1 and solve
for the image gradient using a modified iterative inverse nufft
stochastic differential equation representing the trajectory of (section 3.3). A weighted sum of xt−1 and ∇xt−1 yields the
a sample from our data distribution into a sample from our corrected image xt−1 . This is repeated until t = 0.
noise distribution:
we simulated spiral acquisition by retrospectively interpolat-
dx = f (x, t)dt + g(t)dw. (1) ing in k-space to attain complex-valued measurements along
generated spiral trajectories.
Here, functions f (·, t) and g(·) are the drift and diffusion co-
efficiets of x(t) respectively, and w is a standard Wiener pro- 3.2. Generating spiral trajectories
cess, or Brownian motion.
In order to generate a novel sample from our data distri- Following Kim et al. [13], we consider spiral trajectories of
bution, we can generate a random noise vector and attempt to the form
Z τ
solve the reverse-time SDE, but this is generally intractable. 1
k(τ ) = dϕejωτ ≈ λτ α ejωτ . (3)
However, we can approximate this process by estimating the 0 ρ(ϕ)
noise-conditioned score function, ∇x log pt (x), which com-
Here, ρ denotes sampling density, τ is a function of time, ϕ is
putes the likelihood of a sample x existing between the noise
angular position, ω = 2πn is frequency, with n the number of
and image distributions. With this, the reverse-time SDE be-
turns in k-space, λ a scaling factor equal to matrix size/(2 ∗
comes:
F OV ), and α is a bias term for oversampling the center of k-
dx = f (x, t) − g(t)2 ∇x log pt (x) dt + g(t)dw̄ (2)
  space relative to the edges. Solving this parametric equation
under the constraints of capped gradient slew rate and capped
The score function can be trained using a score matching with gradient amplitude yields gradients (gx (t) and gy (t)) as well
Langevin dynamics algorithm [10, 11]. as a spiral trajectory in the kx ,ky -plane (figure 1). In doing so,
we can tune sampling parameters to control for factors such as
read out duration and dwell time, while varying the number of
3. METHODS interleaves and ratio of low-to-high frequency oversampling.

This research study was conducted retrospectively using hu-

3.3. Image reconstruction is inverse problem solving
man subject data made available in open access by [12]. Eth-
ical approval was not required. MRI undersampled acquisition amounts to measuring an un-
known signal x through some imperfect sampling function A:
3.1. Data y = Ax + ϵ. Here, y is the measured multicoil k-space data,
and A is the non-uniform fourier transform. ϵ is measurement
We use the NYU FastMRI dataset [12] consisting of 6970 noise and exists in the same domain as the y; in MRI, noise
fully sampled 2D brain scans on hardware ranging from 4 is gaussian-distributed across the real and imaginary compo-
to 24 coils. For training and testing, we consider axial T2 nents of y for each coil.
weighted turbo spin echo sequences characterized by the Reconstruction is an ill-posed inverse problem of recover-
following sequence parameters: scan time=140 s, TR=6 ing an image signal x from a set of incomplete k-space mea-
s, TE=113 ms, slices=30, slice thickness=5 mm, field of surements y. As x and y exist in different domains, x is hid-
view=22 cm, matrix size=320x320. Effective scan time for den behind a sampling operator A. Solving this problem ne-
a 2D slice at 2562 resolution is 140s/320 ∗ 256/30 ≈ 3.7s. cessitates prior knowledge. In our case, we learn an underly-
As this data is initially acquired using Cartesian sequences, ing conditional distribution of images and seek to reconstruct
Fig. 3. Representative reconstruction results for a single 2D
16 coil image. Retrospective k-space data was sampled with
an optimized 23 interleave sequence with a total readout dura-
tion of 0.02 s. Rows 1 and 2 show the RSS-reconstructed im-
ages and log-scaled k-space magnitudes for the ground truth,
inverse nufft, and proposed model reconstructions. Below
are the individual coil magnitude and phase images for the Fig. 4. We performed a grid hyperparameter search over a
fully sampled image, the inverse nufft reconstructions, and 2D trajectory space. We fixed readout duration at 0.02 sec-
the model predictions. onds and varied the number of interleaves from 1 to 125
and alpha from 1 to 4. Based on structural similarity of the
samples from this distribution consistent with the measure- model-reconstructed images, we found multiple trajectories
ments. Information is supplied in two forms: first, we learn a that yield improved image quality. In comparison, the naive
conditional score function ∇x log pt (xt |Ã−1 y0 ), where y0 is Archimedean spiral, corresponding to 1 interleave and α = 1,
the measurement in frequency space and Ã−1 is an approxi- performs very poorly.
mate inverse of A, in our case the inverse nufft solved itera-
tively using conjugate gradients. We find that adding this su- In practice, due to the non-invertibility of the nufft, imper-
pervision during training helps to constrain the model when fections in the approximate inverse nufft bleed into the final
faced with a large number of input image channels and the image reconstruction, introducing artifacts and reducing qual-
periodic ambiguity inherent when operating on complex num- ity. To avoid this, we anneal the guidance signal following an
bers. empirically chosen linear schedule γ(t) = β(1 − t), ensuring
Second, we use frequency space gradients to weakly strong guidance at the outset of sampling and minimal arti-
guide the sampling process. At each time step during sam- facts at the end of sampling. A consequence of this choice is
pling, we compute the forward nufft of an uncorrected x̃t−1 that we do not strongly enforce that Ax0 = y0 at time 0.
and take a difference between that and the measurements
y0 . To minimize this difference, the approximate gradient in
4. RESULTS
image space is calculated by solving a modified approximate
−1
inverse nufft Ãt , which corrects for low frequency biases Model reconstruction performance was evaluated on a held-
and applies time step-dependent noising determined by the out test dataset. Test trajectories have a fixed readout duration
noise schedule σ(t), which is necessary for sampling with of 0.02 seconds, in which time the measurements needed to
langevin diffusion. reconstruct a 256x256 pixel, 22x22 cm2 2D image are ac-
quired. Reconstructed image quality was scored using struc-
−1 Ã−1 xt
Ãt xt = + N (0, σ(t)2 ) (4) tural similarity (SSIM) (figure 3).
c1 e−c2 r2
To investigate the effect choosing different scanning tra-
Following [11], we choose a linear noise schedule and ob- jectories has on the quality of reconstructed images, we
serve that the underlying ordinary differential equation de- also perform a grid hyperparameter search of spiral trajec-
scribing transit from latent to image is locally linear, so sum- tories with a fixed readout duration of 0.02 seconds (figure
mation of x̃t−1 and ∇xt−1 to obtain a frequency-corrected 4) and varying α and interleaves. Surprisingly, the com-
image xt−1 is akin to gradient descent (figure 2). mon ’naive’ trajectory, a single interleave Archimedean spi-
Fig. 5. (A) Effect of sampling trajectory optimization, model reconstruction without frequency guidance, and model recon-
struction with frequency guidance. For the non-optimized trajectory, we used a single interleave Archimedean spiral with a
readout duration of 0.02 s. The optimized trajectory uses a 23 interleave, α = 1.23 sequence with an identical readout duration.
(B) Snapshots of the image latent xt and the gradient signal ∇xt taken during a diffusion sampling process.

ral, corresponding to α = 1, performs very poorly when to acquire prospective data with custom sequences and use
sampled below the Nyquist limit. Trajectories which per- it to validate image reconstruction. Currently, a concrete di-
form better tended to lie along two logarithmic curves rect comparison is between the proposed sequences and their
roughly characterized by α = 1.33 log(0.39 interleaves) Nyquist-sampled counterparts, which run roughly 3x longer.
and α = 0.87 log(0.54 interleaves). Apart from the short-term task of matching sequence pa-
Finally, we conduct an ablation study to disentangle the rameters between spiral and Cartesian sequences, the choice
effects of optimal sampling trajectory without model recon- of spiral sequence leaves a considerable amount of flexibil-
struction, model reconstruction without frequency guidance, ity even within the space of optimized interleave and density
and model reconstruction with frequency guidance (figure 5). pairs identified above. Variation in number of interleaves, and
We find that all three contribute to noticeable increases in im- by extension the duration of a single interleave, allows for tai-
age quality, both visual, and quantitative based on SSIM. The loring of sequence contrast, signal, and speed to task-specific
combination of choosing an optimal trajectory, performing requirements. An area in which these sequences could pro-
model reconstruction with conditioning, and using annealed vide additional benefit, even outside of sheer acceleration,
frequency guidance results in large improvements in image would be in imaging tissues with a very short T2* , as acceler-
quality, up to and exceeding a 0.15 boost in SSIM. ation within an interleave allows for proportionally more data
acquisition to occur before signal has decayed.

5. DISCUSSION
6. CONCLUSION
While initial results are promising, the main limitation of this
project is the reliance on retrospective, Cartesian-sampled Here we introduce a new method and show preliminary re-
data. Implementing the sequences outlined in this work will sults for reconstructing spiral MRI using a diffusion model.
likely require customizing spiral sequences so as to match Combining multicoil imaging, spiral scanning, and under-
the contrast and signal of the original Cartesian sequences, sampling enables dramatically faster imaging speeds. Ap-
which will constrain the space of realizable trajectories. Until plications of this work are widespread; in addition to the
a dataset of raw non-Cartesian MRI data becomes available, numerous typical benefits associated with faster scanning,
this will continue to be an obstacle. For a similar reason, including better patient compliance and fewer motion arti-
it is difficult to make head-to-head comparisons between facts, these methods have the potential to reach the extremely
the original sequence and the proposed sequences without high acceleration factors necessary to achieve high resolution
prospective validation. For this reason, the immediate task is real-time 3D imaging.
 7. ACKNOWLEDGMENTS [10] Yang Song, Jascha Sohl-Dickstein, Diederik P
Kingma, Abhishek Kumar, Stefano Ermon, and Ben
No funding was received for conducting this study. The au- Poole, “Score-based generative modeling through
thors have no relevant financial or non-financial interests to stochastic differential equations,” arXiv preprint
disclose. arXiv:2011.13456, 2020.
[11] Tero Karras, Miika Aittala, Timo Aila, and Samuli
8. REFERENCES Laine, “Elucidating the design space of diffusion-based
generative models,” Advances in Neural Information
[1] Mark J Blum, Michael Braun, and Dov Rosenfeld, “Fast Processing Systems, vol. 35, pp. 26565–26577, 2022.
magnetic resonance imaging using spiral trajectories,”
in Medical Imaging. SPIE, 1987, vol. 767, pp. 40–46. [12] Florian Knoll, Jure Zbontar, Anuroop Sriram, Matthew J
Muckley, Mary Bruno, Aaron Defazio, Marc Parente,
[2] Stefanie Winkelmann, Tobias Schaeffter, Thomas Krzysztof J Geras, Joe Katsnelson, Hersh Chandarana,
Koehler, Holger Eggers, and Olaf Doessel, “An optimal et al., “fastmri: A publicly available raw k-space and
radial profile order based on the golden ratio for time- dicom dataset of knee images for accelerated mr image
resolved mri,” IEEE transactions on medical imaging, reconstruction using machine learning,” Radiology: Ar-
vol. 26, no. 1, pp. 68–76, 2006. tificial Intelligence, vol. 2, no. 1, pp. e190007, 2020.
[3] Daniel K Sodickson and Warren J Manning, “Simul- [13] Dong-hyun Kim, Elfar Adalsteinsson, and Daniel M
taneous acquisition of spatial harmonics (smash): fast Spielman, “Simple analytic variable density spiral de-
imaging with radiofrequency coil arrays,” Magnetic res- sign,” Magnetic Resonance in Medicine: An Official
onance in medicine, vol. 38, no. 4, pp. 591–603, 1997. Journal of the International Society for Magnetic Reso-
nance in Medicine, vol. 50, no. 1, pp. 214–219, 2003.
[4] Michael Lustig, David Donoho, and John M Pauly,
“Sparse mri: The application of compressed sensing for
rapid mr imaging,” Magnetic Resonance in Medicine:
An Official Journal of the International Society for Mag-
netic Resonance in Medicine, vol. 58, no. 6, pp. 1182–
1195, 2007.

[5] Julio A Oscanoa, Matthew J Middione, Cagan Alkan,

Mahmut Yurt, Michael Loecher, Shreyas S Vasanawala,
and Daniel B Ennis, “Deep learning-based reconstruc-
tion for cardiac mri: A review,” Bioengineering, vol. 10,
no. 3, pp. 334, 2023.

[6] Yang Song, Liyue Shen, Lei Xing, and Stefano Er-
mon, “Solving inverse problems in medical imaging
with score-based generative models,” arXiv preprint
arXiv:2111.08005, 2021.

[7] Patricia M Johnson, Michael P Recht, and Florian Knoll,

“Improving the speed of mri with artificial intelligence,”
in Seminars in musculoskeletal radiology. Thieme Med-
ical Publishers, 2020, vol. 24, pp. 012–020.

[8] Zhuo-Xu Cui, Chentao Cao, Jing Cheng, Sen Jia,

Hairong Zheng, Dong Liang, and Yanjie Zhu, “Spirit-
diffusion: Self-consistency driven diffusion model for
accelerated mri,” arXiv preprint arXiv:2304.05060,
2023.

[9] Asad Aali, Marius Arvinte, Sidharth Kumar, and

Jonathan I Tamir, “Solving inverse problems with score-
based generative priors learned from noisy data,” arXiv
preprint arXiv:2305.01166, 2023.