#### NOTE
New to DeepInverse? Get started with the basics with the
[5 minute quickstart tutorial](https://deepinv.org/auto_examples/basics/demo_quickstart.html.md#sphx-glr-auto-examples-basics-demo-quickstart-py)..
# Reconstruct accelerated non-Cartesian breast MRI acquisition data
This example reconstructs raw non-Cartesian multicoil kspace data from the FastMRI breast dataset [1](#footcite-solomonfastmri2025), for mammography.
The data was acquired using radial sampling, with 288 spokes. We compare image reconstruction
using all spokes vs. with 72 spokes (i.e. 4x acceleration).
We model the 2D non-uniform FFT physics with [`deepinv.physics.NonCartesianMRI`](https://deepinv.org/api/stubs/deepinv.physics.NonCartesianMRI.html.md#deepinv.physics.NonCartesianMRI), which uses [MRI-NUFFT](https://github.com/mind-inria/mri-nufft).
#### NOTE
This example requires the `mri-nufft` library to model the physics.
Install with `pip install mri-nufft[finufft]` (CPU or MPS) or `pip install mri-nufft[cufinufft]` (GPU).
You can choose between the various backends, see mri-nufft docs for more details. We suggest using `backend='cufinufft'` for cuda devices,
or `backend='finufft'` for CPU. For MPS, use `backend='mps'`, which uses `finufft` but bypasses a torch multithreading problem.
```Python
import torch
import deepinv as dinv
device = dinv.utils.get_device()
import importlib
if importlib.util.find_spec("mrinufft") is None:
raise ImportError(
"mri-nufft is required for NonCartesianMRI. Install with `pip install mri-nufft[finufft]` (CPU or MPS) or `pip install mri-nufft[cufinufft]` (GPU)."
)
if torch.device(device).type == "cuda":
backend = "cufinufft"
elif torch.device(device).type == "mps":
backend = "mps"
else:
backend = "finufft"
```
```none
Selected GPU 0 with 15237.25 MiB free memory
```
## Load the non-Cartesian data
The data is originally provided as `h5` files, where each file consists of raw kspace for one patient volume, of shape:
- `C` = channels = 2,
- `S` = number of spokes = 288,
- `Y` = number of samples per spoke = 640,
- `N` = num coils = 16,
- `P` = partitions = 83 (this will become part of slices = 192)
The data was acquired with a Cartesian fully-sampled partition (slice) readout, and is stored in the frequency domain.
Therefore, preprocessing steps are needed to zero-fill the partition axis, then iFFT the partition axis, then take a 2D slice.
For the demo, we perform these steps offline and provide on HuggingFace a sample slice available to download.
To reproduce this slice preprocessing, you can run the following code:
```python
import h5py
with h5py.File("/path/to/fastMRI_breast_001_1.h5", "r") as f:
y = torch.from_numpy(f["kspace"][:, :, :]).float() # C S Y N P
```
```Python
# y = torch.complex(y[0], y[1]).permute(3, 2, 0, 1) # P N S Y
# P, N, num_shots, num_samples = y.shape
# shift = 192 // 2 - 31
# W = dinv.utils.MRIMixin.ifft(torch.eye(102, dtype=y.dtype), dim=(0,)) # Z,Z centered ifft basis
# y = torch.einsum("p,pnsy->nsy", W[96, shift : shift + P], y) # Take middle slice -> N S Y
# y = dinv.utils.MRIMixin.to_torch_complex(y).permute(1, 2, 3, 0) # 2SYN
dinv.utils.download_example(
"fastMRI_breast_001_1_slice_96.pt", dinv.utils.get_data_home() / "fastMRI_breast"
)
y = torch.load(
dinv.utils.get_data_home() / "fastMRI_breast" / "fastMRI_breast_001_1_slice_96.pt"
).to(device)
print(
"Provided slice shape (2, S=num shots, Y=num samples per shot, N=num coils):",
y.shape,
)
```
```none
/local/jtachell/deepinv/deepinv/examples/external-libraries/demo_mrinufft_breast.py:77: DeprecationWarning: Function 'get_data_home' is deprecated and will be removed in a future version.
"fastMRI_breast_001_1_slice_96.pt", dinv.utils.get_data_home() / "fastMRI_breast"
/local/jtachell/deepinv/deepinv/examples/external-libraries/demo_mrinufft_breast.py:81: DeprecationWarning: Function 'get_data_home' is deprecated and will be removed in a future version.
dinv.utils.get_data_home() / "fastMRI_breast" / "fastMRI_breast_001_1_slice_96.pt"
Provided slice shape (2, S=num shots, Y=num samples per shot, N=num coils): torch.Size([2, 288, 640, 16])
```
The final kspace should be of shape `(1,2,N,S)` to be used with [`deepinv.physics.NonCartesianMRI`](https://deepinv.org/api/stubs/deepinv.physics.NonCartesianMRI.html.md#deepinv.physics.NonCartesianMRI).
```Python
y = y.reshape(2, y.shape[1] * y.shape[2], y.shape[3]).swapaxes(-2, -1).unsqueeze(0)
print("Ready slice shape (1, 2, N, S*Y):", y.shape)
```
```none
Ready slice shape (1, 2, N, S*Y): torch.Size([1, 2, 16, 184320])
```
## Fully-sampled reconstruction
We compute the root-sum-squares reconstruction with all 288 angles, bypassing the need to estimate coil maps.
This reconstruction consists of the adjoint of density compensated kspace, and can be seen as an approximate pseudo-inverse.
The FastMRI data was acquired with golden-angle radial sampling, with 640 samples per shot, and sampling from center outwards.
We use the standard reconstruction size of 320\*320.
```Python
physics_fs = dinv.physics.NonCartesianMRI(
img_size=(320, 320),
num_shots=288,
num_samples_per_shot=640,
coil_maps=y.shape[2], # = 16 to construct physics correctly. Not used for RSS
trajectory="radial",
tilt="golden",
in_out=True,
backend=backend,
device=device,
normalize=True,
)
with torch.no_grad():
x = physics_fs.A_dagger(y, density_compensate=True, rss=True) # 1, H, W
```
```none
/local/jtachell/deepinv/deepinv/.pixi/envs/docs/lib/python3.12/site-packages/mrinufft/_utils.py:67: UserWarning: Samples will be rescaled to [-pi, pi), assuming they were in [-0.5, 0.5)
warnings.warn(
Power iteration converged at iteration 9, ||A^T A||_2=117.25
```
## Reconstruct accelerated data
We undersample the radial data by taking the first few (golden angle) shots.
We construct the accelerated non-Cartesian MRI physics, and we empirically normalise the physics for
solvers that require this.
Then, we estimate the coil sensitivity maps using ESPIRiT (on a lower dimensional image).
```Python
undersampling_factor = 4
y = y[..., : y.shape[-1] // undersampling_factor]
physics = dinv.physics.NonCartesianMRI(
img_size=(320, 320),
num_shots=288 // undersampling_factor,
num_samples_per_shot=640,
coil_maps=y.shape[2], # dummy. Update with real maps below
trajectory="radial",
tilt="golden",
in_out=True,
backend=backend,
device=device,
normalize=True,
noise_model=dinv.physics.GaussianNoise(0.001),
)
coil_maps = physics.estimate_coil_maps(y, method="espirit", decim=4)
physics.update(coil_maps=coil_maps)
```
```none
Power iteration converged at iteration 12, ||A^T A||_2=29.31
0%| | 0/10 [00:00
## Reconstruct with conjugate-gradient
We reconstruct the data with the conjugate-gradient (CG) algorithm, which gives a least-squares solution.
Notice that streak artifacts are present, which are expected for CG on undersampled data.
We compare also to the adjoint (which should be overwhelmingly low-freq)
and the density-compensated adjoint (which gives a fast approximation to the least-squares solution by pre-filtering `y`).
```Python
with torch.no_grad():
x_cg = physics.A_dagger(y)
x_adj = physics.A_adjoint(y)
x_dc = physics.A_dagger(y, density_compensate=True)
# Magnitude images for computing metrics
x_cg = dinv.utils.complex_abs(x_cg)
x_adj = dinv.utils.complex_abs(x_adj)
x_dc = dinv.utils.complex_abs(x_dc)
metric = dinv.metric.PSNR(max_pixel=None, min_pixel=None)
dinv.utils.plot(
[x, x_cg, x_adj, x_dc],
titles=[
"Fully-sampled RSS",
"Conjugate-gradient 4x acc",
"Adjoint 4x acc",
"Density-comp adj 4x acc",
],
subtitles=[
"",
f"{metric(x_cg, x).item():.2f}dB",
f"{metric(x_adj, x).item():.2f}dB",
f"{metric(x_dc, x).item():.2f}dB",
],
plot_inset=True,
extract_loc=(0.2, 0.5),
inset_loc=(0.6, 0),
figsize=(12, 4),
)
```
## What’s next?
We didn’t show any model-based on deep learning reconstruction methods in this example to keep the example lightweight.
You can try out other types of reconstruction algorithms listed in the [user guide](https://deepinv.org/user_guide/reconstruction/introduction.html.md#reconstructors).
Interested in non-Cartesian MRI? You can dive deeper into `mri-nufft`, which has [extensive features](https://mind-inria.github.io/mri-nufft/index.html)
such as more advanced trajectories, trajectory estimation, various coil map estimation algorithms or off-resonance correction.
* **References:**
* Eddy Solomon, Patricia Johnson, Zhengguo Tan, Radhika Tibrewala, Yvonne Lui, Florian Knoll, Linda Moy, Sungheon Gene Kim, and Laura Heacock. Fastmri breast: a publicly available radial k-space dataset of breast dynamic contrast-enhanced mri. *Radiology Artificial Intelligence*, 2025. [doi:10.1148/ryai.240345](https://doi.org/10.1148/ryai.240345).
**Total running time of the script:** (0 minutes 8.183 seconds)
🏷 Tags: MRI
[`Download Jupyter notebook: demo_mrinufft_breast.ipynb`](https://deepinv.org/auto_examples/external-libraries/_downloads/3456e8709b3bcf09ac8086ea5cfcb925/demo_mrinufft_breast.ipynb)
[`Download Python source code: demo_mrinufft_breast.py`](https://deepinv.org/auto_examples/external-libraries/_downloads/2a8350da703287b878b7fd0eabae0ac6/demo_mrinufft_breast.py)
[`Download zipped: demo_mrinufft_breast.zip`](https://deepinv.org/auto_examples/external-libraries/_downloads/3325618e91c2fab9bee2a325c09fa669/demo_mrinufft_breast.zip)
[Gallery generated by Sphinx-Gallery](https://sphinx-gallery.github.io)