installation & examples

[mzaiss]

This looks very nice here.
Thanks for the PR MRsources/MRzero-Core#147
I wanted to test your examples but get weird errors after installation.

from mrpro.algorithms.reconstruction import DirectReconstruction
from mrpro.data import KData

throws the errors:

Cell In[14], line 11
      9 from cmap import Colormap
     10 from einops import rearrange
---> 11 from mrpro.algorithms.reconstruction import DirectReconstruction
     12 from mrpro.data import KData
     13 from mrpro.data import SpatialDimension

File ~\AppData\Local\Packages\PythonSoftwareFoundation.Python.3.11_qbz5n2kfra8p0\LocalCache\local-packages\Python311\site-packages\mrpro\__init__.py:2
      1 from mrpro._version import __version__
----> 2 from mrpro import algorithms, operators, data, phantoms, utils
      3 __all__ = [
      4     "__version__",
      5     "algorithms",
   (...)      9     "utils"
     10 ]

File ~\AppData\Local\Packages\PythonSoftwareFoundation.Python.3.11_qbz5n2kfra8p0\LocalCache\local-packages\Python311\site-packages\mrpro\algorithms\__init__.py:3
      1 """Algorithms for reconstructions, optimization, density and sensitivity map estimation, etc."""
----> 3 from mrpro.algorithms import csm, optimizers, reconstruction, dcf
      4 from mrpro.algorithms.prewhiten_kspace import prewhiten_kspace
      5 __all__ = ["csm", "dcf", "optimizers", "prewhiten_kspace", "reconstruction"]

File ~\AppData\Local\Packages\PythonSoftwareFoundation.Python.3.11_qbz5n2kfra8p0\LocalCache\local-packages\Python311\site-packages\mrpro\algorithms\optimizers\__init__.py:7
...
    211 def from_ismrmrd_acquisitions(
   (...)    215     additional_fields: Sequence[str],
    216     convert_time_stamp: _convert_time_stamp_type = convert_time_stamp_from_siemens,

AttributeError: module 'ismrmrd' has no attribute 'acquisition'

[ckolbPTB]

that is strange, could you maybe try to install it in a new environment with python 3.12? 3.12 is the version we are using to test the examples run correctly in GHA

[mzaiss]

I can confirm that it works with python 3.12.

Disclaimer: I work on windows with vscode, with mrseq:main and mr0: ckolbPTB/MRzero-Core/tree/extend_mrd, I am neither familiar with mrd nor mrpro.

Further ideas:

1. the tempfile copies the file somewhere in windows, as this is quite interesting intermediate data, I would rather copy this to the local folder, or into the output folder mr0 uses: tmp = tempfile.TemporaryDirectory(dir='./output'). Only with this I realized the next point.
2. Non-cartesian do have a header, Cartesian seqs not? is this standard in mrd? It feels a bit unharmonic. Especially if you think of eddy currents or other deviations that are not part of the sequence, but can be of a realistic simulation. Can I trigger a header for Cartesian?
3. The trajectory is extracted from pypulseq now, right?, mr0 uses an own extraction from the seq file. we saw sometimes deviations between pp and ours, especially for spirals with ours being more reliable. Must be checked again, as maybe resolved with all pp work. Yet it is probably cleaner to reconstruct with the exact trajectory the mr0 simulation assumes. Otherwise there could be artifacts solely due to this sim/reco mismatch which is then hard to find.
4. I tried to reconstruct simulated undersampled data (grappa undersampling R=2) with SENSE and compared it to grappa. The recon works, but is not as good as grappa. when comparing the estimated csm to the csm used in the phantom I see a size mismatch that i cannot explain.

here the reco comparisons:

5. the dimension flips as seen in above images. probably an mr0 bug, though. or different convention.

if you share ground thruth csm data of your choice i can create an example

[ckolbPTB]

Thank you for getting it to work. Regarding your points, here are some thoughts:

1. One reason to save it as a tempfile is to make sure it is always nicely cleaned up. ismrmrd as a bit tricky because if files already exist, new data is added to it rather than replacing it. Needs careful file handling, but open to suggestions.
2. ismrmrd files always have a header. When we export Cartesian pulseq data from the scanner, the ismrmrd header is correct and can be used for reconstruction. For non-Cartesian data this is not the case with pulseq and therefore we create this pseudo-ismrmrd-header with the correct information. Once we export the raw data from the scanner, we then replace some header information to "correct" it. We also do this in the examples to test that the pseudo-header we create contains the correct information.
3. In most cases we calculate the trajectory by hand (see e.g. here:

   mrseq/src/mrseq/scripts/radial_flash.py

   Line 263 in a1727e8

   radial_trajectory[:, 0] = k_radial_line * np.cos(rotation_angle_rad)

   ) Using this in the example allows us to test, that our manual calculation is correct.
4. I did try for a 2D case to use a SENSE reconstruction (this example is not online yet) and this worked well. Happy to have a look into your example case. I used coil maps from mrpro like this:

from mrpro.data import SpatialDimension
from mrpro.phantoms.coils import birdcage_2d
im_dims = SpatialDimension(z=1, y=image_matrix_size[1], x=image_matrix_size[0])
coil_maps = birdcage_2d(6, image_dimensions=im_dims, relative_radius=0.8)

phantom = mr0.util.load_phantom(image_matrix_size)
phantom.B1[:] = 1.0
phantom.coil_sens = rearrange(coil_maps[0, ...], 'coils z y x -> coils x y z') 

One thing to keep in mind of course is, that the coil maps will only be estimated where there is signal in the image data. Any extrapolation is not done. This might lead to differences in input and estimated coil maps.
5. The dimension flip and also the "1-voxel-shift" I saw might have just be a difference in fft-convention.

[mzaiss]

1. Got it. Still with tmp = tempfile.TemporaryDirectory(dir='./output') you also get the temp folders but close-by.
2. ok. so in the t2_multi_echo_se_single_line.ipynb this:

   kdata = KData.from_file(fname_mrd, trajectory=KTrajectoryCartesian())
   kdata.header.encoding_matrix = SpatialDimension(z=1, y=image_matrix_size[1], x=image_matrix_size[0])
   kdata.header.recon_matrix = SpatialDimension(z=1, y=image_matrix_size[1], x=image_matrix_size[0])
   

   is redundant or necessary? Can the traj here be extracted from pulseq?
3. ok.
4. ok, it seems the flip and the pixel shift actually causes this problem, which makes some sense. I manually adjusted csm and then the reco looks better. Regularization makes it worse again, but makes again sense as then an image with different convention is used.
5. see 4. and as I used a small 3D volume with 8 partitions, the "pixel shift" in partition dimension has a huge impact.

I clearly see a "larger" csm estimated than provided, but the difference is 0 within the mask, so it might be some blurring at the edges of the estimated, and no mismatch.
I could now get much better recons, and also could use the phantom csm for recon with best results.
but its still not perfect with my coils, yours behave more stable and lead to a fine recon. As I broadcasted yours in 3D, its probably the shift in z that still needs to be considered.

Seems adjusting our conventions could prevent some headaches.