The Documentation should gives more information on the mathematics of the NUFFT, e.g. the gridding/spreading of non uniform points, the Oversampled grid, the deconvolution etc ..

Wild NUFFT appears! Catch them all!

https://github.com/spinicist/riesling

It is CPU-based, With CLI similar to BART (but aims for better 3D support)

For now imports seem really slow, make it faster for faster debug cycles

With the common interface for the numerous backend library, it would be silly not to implement a nice benchmark.

The Benchmark should focuses on the speed performance under various setup:

• monocoil
• multicoil with SENSE
• multicoil without SENSE
• with/without Density Compensation
• Kernel interpolation size.
• Density of the non uniform points (uniform sampling density is not relevant in MRI)

The accuracy of the methods is already tested by the backend themselves, and to a lesser extend in our unit-tests.

The current python code in README is a bit off-putting, I wonder why. Documenting the issue here so that we fix it:

Here is an image, notice the underline and pretty much no color theme:

Another non-cartesian trajectory that has some nice properties

An example:

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8616809/

Wild NUFFT backend appears! Gotta Catch them all

https://gitlab.mpcdf.mpg.de/mtr/ducc

Similar to finufft, but aims to be faster and more accurate (Benchmark is needed)

The goal of this issue is to list and discuss trajectories from the literature to implement. I gathered a bunch of articles, and I would be curious to know

• Which ones should be prioritized ? (would be nice for your future works)
• Which ones should not be reproduced ? (out-of-scope, too demanding computationally or in term of work, etc)

Feel also free to raise new questions, suggest new papers for discussion etc. and I will edit this post.

Here is a list of tasks related to trajectories from the literature that are not already implemented in MRI-NUFFT, along with a checklist for stuff that needs further exploration:

Trajectories

• 2D TWIRL from

  • "Twisting radial lines with application to robust magnetic resonance imaging of irregular flow", John I. Jackson et al. (1992)

• 3D radial trajectory from

  • "A Strategy for Sampling on a Sphere Applied to 3D Selective RF Pulse Design", Sam T. S. Wong and Mark S. Roos (1994)
  • "Fast isotropic volumetric coronary MR angiography using free‐breathing 3D radial balanced FFE acquisition", Stehning, C. et al. (2004)
  • "A radial sampling strategy for uniform k-space coverage with retrospective respiratory gating in 3D ultrashort-echo-time lung imaging", Jinil Park et al. (2016)

• 3D TPI from

  • "Fast Three Dimensional Sodium Imaging", Fernando E. Boada (1997)

• 2D Fibonacci spiral (implemented in #105) from

  • "Uniform K-Space Sampling with an Interleaved Fibonacci Spiral Acquisition", Harvey E Cline and Thomas R Anthony (1999)

• 2D WHIRL from

  • "An Optimized Center-Out k-Space Trajectory for Multishot MRI: Comparison With Spiral and Projection Reconstruction", James G. Pipe (1999)

• 2D Teardrop from

  • 'Teardrop, a novel non-raster readout for True FISP", Christopher Kumar Anand et al. (2001)
  • "Optimizing Teardrop, an MRI Sampling Trajectory", Christopher Kumar Anand et al. (2008)

• 3D genetic trajectory from

  • "Three dimensional k-space trajectory design using genetic algorithms", Sebastian Sabat et al. (2003)

• 3D guided missile trajectory from

  • "Fast Three Dimensional k-space Trajectory Design Using Missile Guidance", Ideas R. Mir et al. (2004)

• Improve 2D waves based on Bunched Phase Encoding from

  • "Bunched Phase Encoding (BPE): A New Fast Data Acquisition Method in MRI", Hisamoto Moriguchi and Jeffrey L. Duerk (2006)

• 3D genetic trajectory from

  • "Random Volumetric MRI Trajectories via Genetic Algorithms", Andrew Thomas Curtis and Christopher Kumar Anand (2008)

• 3D DURGA from

  • "DURGA: A heuristically-optimized data collection strategy for volumetric Magnetic Resonance Imaging", Christopher Kumar Anand et al. (2008)

• 3D Golden Means from

  • "Temporal Stability of Adaptive 3D Radial MRI Using Multidimensional Golden Means", Rachel W. Chan (2009)
  • "Generalization of three‑dimensional golden‑angle radial acquisition to reduce eddy current artifacts in bSSFP CMR imaging", Alexander Fyrdahl et al. (2021)
  • Golden Means from "Golden-Angle Radial MRI: Basics, Advances, and Applications", Li Feng (2022)

• 3D "Spiral" Phyllotaxis from

  • "Spiral Phyllotaxis: The Natural Way to Construct a 3D Radial Trajectory in MRI" Davide Piccini et al. (2011)
  • "Golden-Angle Radial MRI: Basics, Advances, and Applications", Li Feng (2022)

• 3D RAZER from

  • "RAZER: A Pulse Sequence for Whole-Brain Bolus Tracking at High Frame Rates", Sumeeth V. Jonathan et al. (2014)

• Improve 2D and 3D concentric rings based on

  • "Density-Weighted Concentric Circle Trajectories for High Resolution Brain Magnetic Resonance Spectroscopic Imaging at 7T", Lukas Hingerl et al. (2018)
  • "Three-dimensional, 2.5-minute, 7T phosphorus magnetic resonance spectroscopic imaging of the human heart using concentric rings", William T. Clarke (2022)

• 3D T-Hex spirals from

  • "T-Hex: Tilted hexagonal grids for rapid 3D imaging", Maria Engel et al. (2020)
  • "Mono-planar T-Hex: Speed and flexibility for high-resolution 3D imaging", Maria Engel et al. (2021)

• Reproduce Looping Stars from

  • "Looping Star fMRI in Cognitive Tasks and Resting State", Beatriz Dionisio-Parra et al. (2020)

• AZTEK tools and 3D trajectories from

  • "AZTEK: Adaptive zero TE k-space trajectories", Tanguy Boucneau et al. (2020)

• Distinguish 3D Yarnball from currently implemented 3D Seiffert spirals based on

  • "Quantifying Lung Water Density with Ultrashort Echo Time Yarn-Ball MRI", William Quinn Meadus (2020)

• 3D REPI from (#49)

  • "Three Dimensional Radial Echo Planar Imaging for Functional MRI", Christoph A. Rettenmeier et al. (2022)

• 2D PILOT from

  • "PILOT: Physics-Informed Learned Optimized Trajectories for Accelerated MRI", Tomer Weiss et al. (2021)

• 2D trajectory from

  • "Learning Optimal K-space Acquisition and Reconstruction using Physics-Informed Neural Networks", Wei Peng et al. (2022)

• 2D trajectory from

  • "Bayesian Optimization of Sampling Densities in MRI", Alban Gossard (2023)

• 2D BJORK from

  • "B-spline Parameterized Joint Optimization of Reconstruction and K-space Trajectories (BJORK) for Accelerated 2D MRI", Guanhua Wang et al. (2023)

• 3D SNOPY from

  • "Stochastic optimization of three-dimensional non-Cartesian sampling trajectory", Guanhua Wang et al. (2023)

• 2D and 3D ECCENTRIC from

  • "ECCENTRIC: a fast and unrestrained approach for high-resolution in vivo metabolic imaging at ultra-high field MR", Antoine Klauser et al. (2023)

• 3D PETALUTE from

  • "An Accelerated PETALUTE MRI Sequence for In Vivo Quantification of Sodium Content in Human Articular Cartilage at 3T", Cameron X. Villarreal et al. (2024)
  • "Accerated Preclinical UHF Abdominal T1 Mapping using Novel Rosette Accelerated Preclinical UHF Abdominal T1 Mapping using Novel Rosette Ultrashort Echo Time (PETALUTE) Ultrashort Echo Time (PETALUTE)", Alexandra Lipka et al. (2024)

• 3D TURBINE (#50) from

  • "Multi-shot 3D diffusion MRI sequence for a fast and high-resolution imaging at 3T", Sajjad Feizollah et al. (2024)

• 3D Arc-ZTE from

  • "Arc-ZTE: Incoherent k-space sampling in time using continuously-slewed gradients for exible, dynamic, quiet Zero TE MRI", Shreya Ramachandran et al. (2024)

Check-list

• "Optimal Design of k-Space Trajectories Using a MultiObjective Genetic Algorithm", Brian M. Dale et al. (2004)
• "Single TrAjectory Radial (STAR) Imaging", Gordon E. Sarty (2004)
• "A 3D trajectory for undersampling k-space in MRSI applications", Sergio Uribe et al. (2007)
• "A Fast Method for Designing Time-Optimal Gradient Waveforms for Arbitrary k-Space Trajectories", Michael Lustig et al. (2008)
• "Spiral Trajectory Design: A Flexible Numerical Algorithm and Base Analytical Equations" James G. Pipe et al. (2014)
• "Fast and Robust Design of Time-Optimal k-Space Trajectories in MRI", Mathias Davids et al. (2015)
• "Squeezed Trajectory Design for Peak RF and Integrated RF Power Reduction in Parallel Transmission MRI", Qing Li et al. (2018)
• "Partial fourier shells trajectory for non-cartesian MRI", Tao et al. (2019)
• "High resolution free-breathing respiratory-resolved volumetric lung imaging at 0.55T using stack-of-spiral out-in bSSFP", Ziwei Zhao et al. (2024)
• "Silent Looping Star fMRI with enhanced Encoding and Reconstruction Performance", Ana Beatriz Solana et al. (2024)
• "Fast and Motion-Robust Non-Contrast MR Angiography using Centric k –k Trajectory.", Vadim Malis et al. (2024)
• "ABRICOTINE MRI: Enhancing Sparsity Across the Three Dimensions of the Fourier Domain in Cartesian Sampling", Antoine Klauser et al. (2024)
• "3D UTE With Twisted Trajectories", Michael Carl et al. (2024)

Following the discussion of getkeops/keops#328, adding a proper NDFT nufft (with GPU support! ) should be doable now.

GPU CI is broken (due to pynfft).

Along with the ISMRMRD format (#64) we could also refactor/expose the routine used in the BART interface.

Currently we also do FFT which is not ideal.

Some paradigm of acquisition uses Stack of non cartesian trajectory (e.g. Stack of Spiral, stack of radial, etc). even if we know that fully 3D non cartesian are better, having support for such trajectory for comparison would be nice. In particular the "stacking" can be leverage for faster reconstruction: by performing a sequence of (potentially batched) 2D-NUFFT followed by a 1D FFT (or DFT if only a few slice are acquired).

We need a baseline to compare the results against. For this, we need to have NDFT rightly implemented the same way, which is slow but then can help a lot in benchmarking.

PyNUFFT has support for both CPU and GPU:
https://jyhmiinlin.github.io/pynufft/index.html

CPU bindings is already done, but adding another GPU library bindings won't hurt.

Presented in:

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC9546489/

This is not equivalent to a stack of 2D-radial trajectory due to the interleaving of the blades (see figure)

Before any proper release on pypi, we should provide unit-test of the bindings.

Simply checking each interface vs the NUDFT on the small set of point should be enough and fairly easy. However the parametrization of each interface might be tricky.

A package like pytest-cases might help.

Currently the pipe density needs normalization. The best I feel is to take a template data, do FH ( D * F (x)) and normalize based on difference in means
This can help have consistent lipschitz constants and also a good consitent candidate lamda/mu values for normalization

torch-kbnufft is used in the DL community. Yet better nufft libraries are available (e.g. cufinufft, gpunufft). The idea would be add torch-kbnufft in order to the run the benchmark, and show what is the best nufft currently.

Note that torch-kbnufft has also some interesting capabilities, notably differentiation and support for toeplitz embedding (speeds things up in data consistency evaluation)

After cloning & installing the current repository, I get the following error when trying to import the package as a whole:

Traceback (most recent call last):
File "", line 1, in
File "/home/guillaume/miniconda3/envs/general/lib/python3.10/site-packages/mrinufft/init.py", line 9, in
from .operators import (
File "/home/guillaume/miniconda3/envs/general/lib/python3.10/site-packages/mrinufft/operators/init.py", line 2, in
from .interfaces import (
File "/home/guillaume/miniconda3/envs/general/lib/python3.10/site-packages/mrinufft/operators/interfaces/init.py", line 4, in
from .cufinufft import MRICufiNUFFT, CUFINUFFT_AVAILABLE
File "/home/guillaume/miniconda3/envs/general/lib/python3.10/site-packages/mrinufft/operators/interfaces/cufinufft.py", line 8, in
from .utils.gpu_utils import (
File "/home/guillaume/miniconda3/envs/general/lib/python3.10/site-packages/mrinufft/operators/interfaces/utils/init.py", line 7, in
from .gpu_utils import (
File "/home/guillaume/miniconda3/envs/general/lib/python3.10/site-packages/mrinufft/operators/interfaces/utils/gpu_utils.py", line 53, in
with open(str(Path(file).parent / "css_colors.txt")) as f:
FileNotFoundError: [Errno 2] No such file or directory: '/home/guillaume/miniconda3/envs/general/lib/python3.10/site-packages/mrinufft/operators/interfaces/utils/css_colors.txt'

I guess the css_colors.txt file was simply not pushed.

Hi all,

I have just started trying out MRI-NUFFT after ISMRM, previously having used just BART.

As far as I understand, BART expects non-uniform trajectories to be normalized to 1/FOV. In most cases for me this ends up such that distance between k-space points is just 1 (dk=1), which seems to work.
However, MRI-NUFFT by default normalizes trajectories to [-0.5,0.5] and then scales this based on the requested image size in traj2cfl (traj_ = traj * (np.array(shape) - 1)). This couples the image size to the reconstructed FOV, and generally complicates things. Also, it seems for my BART trajectory, the wrapping (new_traj = (new_traj + 0.5) % 1 - 0.5) step in proper_trajectory essentially corrupts the trajectory, as the later scaling does not undo the wrapping.

Maybe I am missing something, however if not an easy fix would be removing the call to proper_trajectory(samples, normalize="unit") from the BART operator, as the normalization does not seem to fit well with the BART formalism.

Curious to see what others think.

Adding Propeller trajectories would be a nice addition to the current collection. Also this might requires some additional reconstruction methods to benefit from the motion correction.

References:

• https://www.ajronline.org/doi/10.2214/AJR.07.3657
• https://mriquestions.com/propellerblade.html

There are a bunch of scripts to better view k-space sampling trajectories, particularly but not limited to 3D:

• Temporal colored plots: Helps to verify how good trajectory is with respect to off-resonance
• K-space sampling density plots: Helps identify k-space sampling holes
• Plotting only a bunch of shots within a range etc.
• Have a gallery of trajectories presented based on above details. It should be automatic generation from the trajectories library

Ideally, we need a script to handle this and this can be a good GUI to handle most of what we need (with or without the help of napari).
@paquiteau Opinions?

With all the recent development we should deploy a new release. But there is still some work to be done

• Merge #12 for initial tests
• Describe the (hopefully temporary) custom installation of cufinufft
• Merge #10 for IO support
• Add tests for trajectories (Just check that all trajectories generated lies in the [-0.5, 0.5] region, maybe some checks on the constraints).
• Add support for the new gpuNUFFT, patched by @chaithyagr
• [ ] Add support for the nufft in mrrt

Since you're using the finufft package, is it possible to wrap their type-3 nufft?

At ISMRM 2024, peoples asks if there were any Pulseq integration.

I am not familiar with Pulseq, but let's have this here for future reference and discussion.

The documentation is good shape but can be better.

Switch to a better like https://sphinx-book-theme.readthedocs.io/en/latest/index.html (see the config from benchopt for instance)

Clearly separate the contents following the convention of https://diataxis.fr/:

• References (aka the API)
• Explanation (aka the Math and design choices)
• Tutorials (empty for now)
• How-Tos (currently the minimal examples)

We could have possibly the following examples:

• A simple U-Net (pulled from fastMRI dataset possibly), trained? @Lenoush can you handle this in your free time?
• A network built with DeepInv library for some already preset reconstruction (PnP possibly?) @matthieutrs can you handle this this in your free time? I think you can also have something basic passed to @Lenoush and she can work on it forward.
• A simple example where we could learn the trajectory possibly? @alineyyy can you handle this in your free time?

Hi,

I have the following error when using "import snkf" on a windows computer:

import snkf
Traceback (most recent call last):
File "", line 1, in
File "C:\Users\benja\anaconda3\envs\fRMI\Lib\site-packages\snkf_init_.py", line 6, in
from .handlers import (
File "C:\Users\benja\anaconda3\envs\fRMI\Lib\site-packages\snkf\handlers_init_.py", line 25, in
importlib.import_module("." + name, name)
File "C:\Users\benja\anaconda3\envs\fRMI\Lib\importlib_init_.py", line 90, in import_module
return bootstrap.gcd_import(name[level:], package, level)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "C:\Users\benja\anaconda3\envs\fRMI\Lib\site-packages\snkf\handlers\acquisition_init.py", line 3, in
from .base import (
File "C:\Users\benja\anaconda3\envs\fRMI\Lib\site-packages\snkf\handlers\acquisition\base.py", line 28, in
from .workers import acq_cartesian, acq_noncartesian
File "C:\Users\benja\anaconda3\envs\fRMI\Lib\site-packages\snkf\handlers\acquisition\workers.py", line 11, in
from fmri.operators.fourier import FFT_Sense
File "C:\Users\benja\anaconda3\envs\fRMI\Lib\site-packages\fmri\operators\fourier.py", line 12, in
from mrinufft import get_operator
File "C:\Users\benja\anaconda3\envs\fRMI\Lib\site-packages\mrinufft_init.py", line 9, in
from .operators import (
File "C:\Users\benja\anaconda3\envs\fRMI\Lib\site-packages\mrinufft\operators_init_.py", line 21, in
importlib.import_module(".interfaces." + name, name)
File "C:\Users\benja\anaconda3\envs\fRMI\Lib\importlib_init_.py", line 90, in import_module
return _bootstrap._gcd_import(name[level:], package, level)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "C:\Users\benja\anaconda3\envs\fRMI\Lib\site-packages\mrinufft\operators\interfaces\bart.py", line 20, in
BART_AVAILABLE = not subp.call(
^^^^^^^^^^
File "C:\Users\benja\anaconda3\envs\fRMI\Lib\subprocess.py", line 389, in call
with Popen(*popenargs, **kwargs) as p:
^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "C:\Users\benja\anaconda3\envs\fRMI\Lib\subprocess.py", line 1026, in init
self._execute_child(args, executable, preexec_fn, close_fds,
File "C:\Users\benja\anaconda3\envs\fRMI\Lib\subprocess.py", line 1538, in _execute_child
hp, ht, pid, tid = _winapi.CreateProcess(executable, args,
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
FileNotFoundError: [WinError 2] Le fichier spécifié est introuvable

Hi, while looking for other stuff, we found an another NUFFT backend with @chaithyagr:

https://github.com/mritools/mrrt.nufft

Its a interesting one: fully written with a Python Stack (Cython / numba /Cupy) and is basically a port of Fessler's NUFFT from MATLAB. It has:

• both support for CPU and GPU
• single or double precision.
• Multi coil support (but without SENSE reconstruction)

MRI-NUFFT documentation is great, but it lack a document describing the common interfaces for the operators (basically describing the bare minimum API that makes a MRI-NUFFT operator), as well as the expected data format.

This is also the case for the trajectory generation (where we are always in a (Nc,Ns,d) setup).

This would help the community to include new backends/trajectories in MRI-NUFFT, and the other way around, helping the integration of MRI-NUFFT in existing and future tools.

Density compensation weights act as preconditionner on the adjoint operator, and leads to better image quality (and thus faster convergence in iterative schemes)

Currently we support Voronoi estimation and Pipe's iterative scheme. Yet those estimation are simple heuristic, and does optimize a criteria, thus there is no way of knowing if its "best" in one form or the other. Also the current Voronoi method implementation is not applicable in 3D (too much memory and computation requirement).
Jeff Fessler's work: https://web.eecs.umich.edu/~fessler/book/c-four.pdf provides a good overview of the available methods (see notably section 6.4.2)

Implementation of these methods would be a great asset to mri-nufft, as they can be made agnostic to the NUFFT backend. Most of them are iterative, so a benchmark using benchopt could also be done.

The density compensation vector can be estimated using several methods:

• Pipe et al. (Fixed point method)
• Voronoi Cell area
• ...

A good list of methods is provided by Jeff Fessler: https://web.eecs.umich.edu/~fessler/book/c-four.pdf , section 6.4.2

After creating 3D cones trajectories, I have observed that he min value of the trajectory if out f the [-0.5,0.5] range. Here I attached the code to reproduce this issue:

import mrinufft

# Trajectory parameters
Nc = int(208*384/10)  # Number of shots
Ns = 384 # Number of samples per shot
in_out = False  # Choose between in-out or center-out trajectories
tilt = "uniform"  # Angular distance between shots

trajectory = mrinufft.initialize_3D_cones(Nc, Ns, in_out=in_out)

trajectory = trajectory.reshape(Nc*Ns,3)

# here min value out of [-0.5,0.5] range
trajectory.min()
-0.5000150076299451

ISMRMRD is the standard for sharing k-space data acquired using an MRI machine ¹.
There exists Python bindings ² , as well as a toolbox for handling the data ³, the examples in it are quite instructive on how we could read (and write !) ISMRMRD files.

Ideally, MRI-NUFFT should be able to read the kspace-data and trajectories in ISMRMRD files, and process them to make them compatible with our standards. Ideally, it should also be able to write ISMRMRD, for interoperability with other toolboxes.

There is also vendor conversion tools, for instance for Siemens⁴

Here, estimate_density just estimates the density, not the density_compensators. What we need is a reciprocal of this values.

All we need to do is to change this line to:
tf.math.reciprocal_no_nan(tfmri.estimate_density(
samples, shape, method="pipe", max_iter=15
))

Also, some tests on this would be helpful.

The main idea is to implement the differentation of NUFFT forward and adjoint operator with respect to k-space trajectory locations.

Some earlier implementations:

1. TorchKBNUFFT: https://github.com/mmuckley/torchkbnufft/blob/main/torchkbnufft/_autograd/interp.py

2. TFKBNUFFT: https://github.com/zaccharieramzi/tfkbnufft/blob/master/tfkbnufft/kbnufft.py

3. tensorflow-nufft: https://github.com/mrphys/tensorflow-nufft/blob/7d2f502c6f5f4ba305e7a389ef8794b65a5b4f29/tensorflow_nufft/python/ops/nufft_ops.py#L126-L232

4. and 2) are very slow and inefficient. 3) while is much faster, it is based on Google_CUDA and reimplementation of cufinufft.

When implemented in mri-nufft, we expect a much better interface with all the backends, giving a much faster, efficient and flexible way to compute.

Underlying math: https://web.eecs.umich.edu/~fessler/papers/files/jour/23/wang-23-eao.pdf

Currently we are using pyNFFT as the reference backend to check the fidely of the other supported backend.
This result from a speed/precision compromise, because comparing all backends to the NDFT would be expensive

However this causes some troubles:

• as the method are not strictly NDFT, the triangle of comparison NDFT/NFFT/another_backend introduce more imprecision.
• PyNFFT does not seem to be actively maintained, and the installion is broken with new version of cython (see pyNFFT/pyNFFT#75)

An alternative would be to use finufft as the reference implementation, as they have prime support for python bindings.

Interface observed: gpuNUFFT, could potentially impact others also as it seems cupy issue

I recently saw nearly double to triple the usage of GPU memory than what was needed. On closely checking, I think the problem is when we try to convert a cupy array to complex64

This causes creation of new unwanted memory. We need a better converter, which should either delete the earlier memory or not even create one in the ideal scenario that we give the right datasizes. I feel we need a to complex64 as a specific separate utils function to handle this.