Skip to content

Improve NUFFT Gram for large datasets#989

Open
fzimmermann89 wants to merge 8 commits into
mainfrom
nufftgram
Open

Improve NUFFT Gram for large datasets#989
fzimmermann89 wants to merge 8 commits into
mainfrom
nufftgram

Conversation

@fzimmermann89
Copy link
Copy Markdown
Member

@fzimmermann89 fzimmermann89 commented May 1, 2026

This makes it possible to do iterative reconstruction on large non cartesian 3D datasets by

  1. memory efficient in-place operations
  2. avoiding fftshifts
  3. avoiding unnesesarry intermedate allocations

Additionally, the optional DCF weighting is now exposed in NonUniformFastFourierOpGramOp. It was already implemented in the backend, but there was no way to use it. Fourier.gram will still not use it, only if created manually with an additional argument weight

For subpsace compressed reconstruction, one needs PCA^H Fourier^H Fourier PCA
So for, for example, 500 different trajectories in the Fourier operator and 5 subspace coefficients, we would have to do 500 FFTs and save 500 Toeplitz kernels. This would be to slow for iterative recon inside a neural network.
Instead, we can pull the subspace compressin into the kernel. The kernel is not longer a real valued 2D image, but includes off-diagonal mixing; in this example 5x5x(2recon_y)x(2recon_x) and we only require 5 FFTs

Finally, it adds a .toeplitz function on the nufft operator as a .gram with more control about density compensation and subsapce compression. .gram uses it with the default parameters. Gram is a property and cannot take arguments, gram gets chained through, whereas topelitz needs to be called with the options on the nufft operator for more granular control

@github-actions
Copy link
Copy Markdown
Contributor

github-actions Bot commented May 1, 2026

Coverage

Coverage Report
FileStmtsMissCoverMissing
src/mrpro
   _version.py6267%7–8
src/mrpro/algorithms/csm
   inati.py25196%43
src/mrpro/algorithms/dcf
   dcf_voronoi.py55493%15, 55–56, 89
src/mrpro/algorithms/optimizers
   adam.py30680%108, 125–129
   cg.py52198%139
   pdhg.py81396%180–181, 187
   pgd.py53492%108, 153–156
src/mrpro/algorithms/reconstruction
   DirectReconstruction.py30293%64, 80
   Reconstruction.py56886%57–59, 123–126, 128, 136
   RegularizedIterativeSENSEReconstruction.py52296%125, 153
   TotalVariationRegularizedReconstruction.py54493%106, 112, 129, 132
src/mrpro/data
   AcqInfo.py165796%49, 56, 134–135, 137, 243, 367
   Dataclass.py3112692%58, 319, 335, 401, 465–467, 480, 575, 595–596, 598, 613–614, 616, 663–664, 669–670, 855–856, 881, 888, 893–894, 896
   DcfData.py33197%62
   EncodingLimits.py97397%37, 127, 130
   IData.py1781293%138, 155–156, 217–225, 258, 304, 322
   IHeader.py174895%49, 98–101, 284, 288, 292, 296
   KData.py2302490%122–123, 138, 145, 156–167, 178, 186, 197, 236, 258–260, 307–308, 380, 546, 548, 620
   KHeader.py1761393%115–121, 148, 196, 203–204, 231–238
   KNoise.py22195%44
   KTrajectory.py95397%163, 165, 185
   Rotation.py7374294%104, 202, 339, 437, 481, 499, 586, 588, 597, 634, 636, 699, 776, 781, 784, 799, 816, 821, 897, 1085, 1090, 1093, 1117, 1121, 1145, 1265, 1267, 1275–1276, 1340, 1422, 1626, 1633–1635, 1694, 1790, 1942, 1977, 1981, 2156, 2177
   SpatialDimension.py1501987%34, 103, 146, 158, 278–280, 293–295, 329, 347, 360, 373, 386, 399, 408–409, 437
src/mrpro/data/traj_calculators
   KTrajectoryCalculator.py26292%84, 95
   KTrajectoryCartesian.py31487%129–132, 136
   KTrajectoryIsmrmrd.py19195%57
   KTrajectorySpiral2D.py571377%63–66, 69, 71, 73, 75, 77, 105, 107, 134–136
src/mrpro/operators
   AveragingOp.py38295%73, 115
   B0InformedFourierOp.py75297%151–154
   CartesianSamplingOp.py112496%152, 191, 266, 387
   ConjugateGradientOp.py89792%62, 64, 100, 106, 228, 230, 233
   ConstraintsOp.py85495%78, 80, 250, 255
   EndomorphOperator.py28293%71, 77
   FiniteDifferenceOp.py29293%40, 126
   FourierOp.py105991%85–87, 192–193, 212, 257, 321, 326
   Functional.py70297%116, 118
   GridSamplingOp.py1651591%72–73, 82–83, 90–91, 94, 96, 98, 282, 290–291, 306, 312–313
   LinearOperator.py202697%217, 255, 296, 305, 313, 330
   LinearOperatorMatrix.py1741989%99, 136, 169, 178, 183, 192–195, 208–211, 219–220, 225–226, 238, 347, 377, 404
   MultiIdentityOp.py16288%58, 63
   NonUniformFastFourierOp.py3093190%75, 102, 226, 228, 266, 268, 370, 429–430, 461–464, 481, 534, 575–576, 589, 591, 595, 600, 612, 616, 638, 644–647, 662, 665, 720
   Operator.py88397%82, 115, 125
   PatchOp.py49394%93, 129, 144
   ProximableFunctionalSeparableSum.py44393%118, 213, 224
   SliceProjectionOp.py1781094%45, 62, 64, 70, 154, 180, 216, 253, 290, 330
   WaveletOp.py119397%184, 228, 254
   ZeroPadOp.py18194%30
src/mrpro/operators/functionals
   SSIM.py73790%60–80, 82, 86, 114, 147
src/mrpro/operators/models
   EPG.py2561993%31–32, 283, 288, 304–306, 326–327, 332, 356, 361, 386, 391, 545, 599, 738, 755, 782
src/mrpro/phantoms
   EllipsePhantom.py43295%66, 131
   b0map.py262119%37–63
   brainweb.py2974087%276, 290–294, 349–359, 398, 454–457, 479–480, 485–486, 488–489, 493, 501, 508–509, 550, 621, 624–625, 644, 653–656, 667, 669, 700–701, 715, 723
   fastmri.py1061091%50–51, 59, 65, 162, 169–171, 174–175, 189
   m4raw.py74495%58–59, 74, 76
   mdcnn.py71790%58, 62–63, 70, 82, 88, 135
src/mrpro/utils
   RandomGenerator.py1551590%23–24, 36, 38, 188, 212, 428, 446, 528, 799, 829–832, 895, 898
   ema.py42490%49, 74, 78, 82
   filters.py61198%46
   indexing.py177199%321
   pad_or_crop.py39685%40, 43, 46, 49, 66, 73
   reshape.py1531093%176, 370, 482–484, 505, 508–509, 515, 530
   slice_profiles.py49688%21, 37, 119–122, 155
   split_idx.py10280%43, 47
   summarize.py57689%40–41, 70–73, 81
   unit_conversion.py721579%34, 44, 51, 53, 62, 69, 71, 78, 80, 89, 100, 121, 123, 144, 146
TOTAL822752294% 

Tests Skipped Failures Errors Time
3309 0 💤 0 ❌ 0 🔥 2m 32s ⏱️

@fzimmermann89 fzimmermann89 requested a review from ckolbPTB May 1, 2026 17:03
@github-actions
Copy link
Copy Markdown
Contributor

github-actions Bot commented May 1, 2026

📚 Documentation

📁 Download as zip
🔍 View online

Comment thread src/mrpro/operators/NonUniformFastFourierOp.py Outdated
Comment thread src/mrpro/operators/NonUniformFastFourierOp.py
@ckolbPTB
Copy link
Copy Markdown
Collaborator

ckolbPTB commented May 1, 2026

Can we think of a better way to select the loop case? Doing the loop for all 3D reconstructions could make many of these reconstructions unneccessarily slow.

@fzimmermann89
Copy link
Copy Markdown
Member Author

I tried it on cpu and a6000 gpu. for 3d non cartesian reconstruction with 128x128x128 voxels, the loop in the forward saturates the gpu and it is not slower, i.e. batching more compute at once dpes not help.
are there realistically trajectories that are non cartesian on all 3 directions and use less voxels?

The issue is more on the other side, that for some 2D reconstructions it could also make sense to loop (I have not tried that)

@fzimmermann89
Copy link
Copy Markdown
Member Author

fzimmermann89 commented May 3, 2026

mnhm... i reconsidered the loop case. I think i will remove this and instead, one will have to loop outside the operator.
This might be achieved by unpacking and using the LinearOperatorMatrix, which also loops over the operator.

@fzimmermann89 fzimmermann89 requested a review from ckolbPTB May 3, 2026 16:57
@fzimmermann89
Copy link
Copy Markdown
Member Author

For subpsace compressed reconstruction, one needs PCA^H Fourier^H Fourier PCA
So for, for example, 500 different trajectories in the Fourier operator and 5 subspace coefficients, we would have to do 500 FFTs and save 500 Toeplitz kernels. This would be to slow for iterative recon inside a neural network.
Instead, we can pull the subspace compressin into the kernel. The kernel is not longer a real valued 2D image, but includes off-diagonal mixing; in this example 5x5x(2recon_y)x(2recon_x) and we only require 5 FFTs

Finally, it adds a .toeplitz function on the nufft operator as a .gram with more control about density compensation and subsapce compression. .gram uses it with the default parameters. Gram is a property and cannot take arguments, gram gets chained through, whereas topelitz needs to be called with the options on the nufft operator for more granular control

For subpsace compressed reconstruction, one needs PCA^H Fourier^H Fourier PCA
So for, for example, 500 different trajectories in the Fourier operator and 5 subspace coefficients, we would have to do 500 FFTs and save 500 Toeplitz kernels. This would be to slow for iterative recon inside a neural network.

Instead, we can pull the subspace compressin into the kernel. The kernel is not longer a real valued 2D image, but includes off-diagonal mixing; in this example 5x5x(2*recon_y)x(2*recon_x) and we only require 5 FFTs

Additionaly, it adds a .toeplitz function on the nufft operator as a .gram with more control about density compensation and subsapce compression. .gram uses it with the default parameters. Gram is a property and cannot take arguments, gram gets chained through, whereas topelitz needs to be called with the options on the nufft operator for more granular control
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants