espirit csm

[fzimmermann89]

A work-in-progress for espirit algorithm for cartesian images.
This is blocked until we have

• a Cartesian FFT Operator (the code currently does an ifft on kdata.data) and
• a notion of the ACS region ('calib' in the code).

Other TODOs/Notes:

• image_shape should removed by using the shape of the result of the adjoint operator
• the code is completely untested and will most likely contain logic and syntax errors.
• not sure if the unfolds are correct.

The logic is derived from sigpy
we would have to check the license

Why implement it and not wrap sigpy?
Sigpy only works with Cartesian data and is sometimes annoying to install, as it uses cupy for the GPU implementation.

Closes #34

[fzimmermann89]

https://onlinelibrary.wiley.com/doi/full/10.1002/mrm.28386

[github-actions]

Coverage Report

File	Stmts	Miss	Cover	Missing
src/mrpro/algorithms/csm
inati.py	24	1	96%	44
walsh.py	16	1	94%	34
src/mrpro/algorithms/dcf
dcf_voronoi.py	53	4	92%	15, 48–49, 76
src/mrpro/algorithms/optimizers
adam.py	20	1	95%	69
src/mrpro/algorithms/reconstruction
DirectReconstruction.py	28	16	43%	51–71, 85
IterativeSENSEReconstruction.py	13	1	92%	76
Reconstruction.py	50	22	56%	42, 54–56, 80–87, 104–113
RegularizedIterativeSENSEReconstruction.py	41	17	59%	96–100, 114–139
src/mrpro/data
AcqInfo.py	128	3	98%	26, 169, 207
CsmData.py	44	14	68%	17, 104–121, 125–127
DcfData.py	45	8	82%	18, 66, 78–83
IData.py	67	9	87%	119, 125, 129, 159–167
IHeader.py	75	7	91%	75, 109, 127–131
KHeader.py	153	12	92%	25, 119–123, 150, 199, 210, 217–218, 221, 228
KNoise.py	31	15	52%	39–52, 56–61
KTrajectory.py	69	5	93%	178–182
MoveDataMixin.py	137	18	87%	15, 113, 129, 143–145, 207, 305–307, 320, 399, 419–420, 422, 437–438, 440
QData.py	39	7	82%	42, 65–73
Rotation.py	674	35	95%	100, 198, 335, 433, 477, 495, 581, 583, 592, 626, 628, 691, 768, 773, 776, 791, 808, 813, 889, 1077, 1082, 1085, 1109, 1113, 1240, 1242, 1250–1251, 1315, 1397, 1690, 1846, 1881, 1885, 1996
SpatialDimension.py	230	20	91%	33, 103, 135, 141, 261–263, 276–278, 312, 330, 343, 356, 369, 382, 391–392, 407, 416
acq_filters.py	12	1	92%	47
src/mrpro/data/_kdata
KData.py	134	18	87%	109–110, 125, 132, 142, 150, 204–205, 243, 248–249, 268–279
KDataRemoveOsMixin.py	29	2	93%	44, 46
KDataSelectMixin.py	19	2	89%	48, 63
KDataSplitMixin.py	48	3	94%	53, 84, 93
src/mrpro/data/traj_calculators
KTrajectoryCalculator.py	25	2	92%	23, 45
KTrajectoryIsmrmrd.py	13	2	85%	41, 50
KTrajectoryPulseq.py	29	1	97%	54
src/mrpro/operators
CartesianSamplingOp.py	89	3	97%	118, 157, 280
ConstraintsOp.py	60	2	97%	46, 48
EndomorphOperator.py	65	2	97%	228, 234
FiniteDifferenceOp.py	27	2	93%	40, 105
FourierOp.py	157	3	98%	263, 381, 386
Functional.py	71	5	93%	20–22, 117, 119
GridSamplingOp.py	136	9	93%	72–73, 82–83, 90–91, 94, 96, 98
LinearOperator.py	171	12	93%	55, 91, 190, 220, 261, 270, 278, 287, 295, 320, 418, 423
LinearOperatorMatrix.py	158	16	90%	82, 119, 152, 161, 166, 175–178, 191–194, 203, 215, 304, 331, 359
MultiIdentityOp.py	13	2	85%	43, 48
Operator.py	78	2	97%	25, 74
ProximableFunctionalSeparableSum.py	39	3	92%	50, 103, 110
SliceProjectionOp.py	173	8	95%	44, 61, 63, 69, 206, 227, 260, 300
WaveletOp.py	120	5	96%	152, 170, 205, 210, 233
ZeroPadOp.py	16	1	94%	30
src/mrpro/utils
filters.py	62	2	97%	44, 49
slice_profiles.py	46	6	87%	20, 36, 113–116, 149
sliding_window.py	34	1	97%	34
split_idx.py	10	2	80%	43, 47
summarize_tensorvalues.py	11	2	82%	23, 25
typing.py	18	11	39%	8–23
zero_pad_or_crop.py	31	6	81%	26, 30, 54, 57, 60, 63
TOTAL	4908	352	93%

Tests	Skipped	Failures	Errors	Time
2010	0 💤	1 ❌	0 🔥	2m 3s ⏱️

[SherineBrahma]

Background:

ESPIRiT paper

In eq.(9), we need to calculate V_|| that is associated with the row space of the Calibration matrix A. This is later used to construct the forward matrix W. In principle, V_|| is obtained from V by trimming off the row vectors corresponding to zero eigen-values. However, here, they have done it based on this mask: "S > thresh * S.max()". This makes sense as eigen-values may not go to 0 because of the presence of noise.

Problem:
This trimming of VH[S > thresh * S.max()] based on the mask (S > thresh * S.max()) is not supported inside "vmap" operations. I have a feeling that operations like this change the size of the matrix, and that poses challenges to perform batch operations, given that all elements of the batches should have the same shape.

Possible Solution:
If we want to adopt a parallelized approach like before, maybe it makes sense to perform vmap only till the SVD is calculated, since SVD is usually a computationally demanding algorithm for large matrices. Thereafter, we can use loops to construct V_||'s, where different V_||'s can have different shapes based on the batched SVD results. For this, I can use the smap function defined here. However, based on the layout of this file (define _iterative_walsh_csm and then estimate_walsh, define estimate_espirit and then _espirit_csm) , I am not sure if opening up this code to perform first part in one place (SVD calculation using vmap), while performing the other part in another place (forward operation using loops), will break the stye of the code.

If we want to retain the original coding style, I can use the smap function over the whole _espirit_csm function where the SVD is included, but then we would not be parallelizing the SVD.

@ckolbPTB

[fzimmermann89]

I would just switch to fully serial calculation of the csms for different others using smap.

[github-actions]

📚 Documentation

📁 Download as zip
🔍 View online