Add simulation of B1+-heterogeneity

[JanWP]

This is a first attempt at adding a B1+ map to the phantom structure and
adapt the simulations accordingly. I still need to add documentation and
a specific test case.
Closes #475

[JanWP]

I'll look into the test failures. I have trouble running the CUDA tests on my machine, so didn't run them before the pull request (simulations with CUDA seem to run fine).

[codecov]

Codecov Report

Attention: Patch coverage is 63.15789% with 7 lines in your changes missing coverage. Please review.

Project coverage is 85.76%. Comparing base (927de27) to head (4ea5c51).

Files with missing lines	Patch %	Lines
KomaMRIPlots/src/ui/DisplayFunctions.jl	50.00%	4 Missing ⚠️
...RICore/src/simulation/SimMethods/Bloch/BlochGPU.jl	0.00%	3 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #482      +/-   ##
==========================================
- Coverage   89.86%   85.76%   -4.11%     
==========================================
  Files          54       54              
  Lines        3050     3056       +6     
==========================================
- Hits         2741     2621     -120     
- Misses        309      435     +126     

Flag	Coverage Δ
base	86.24% <ø> (ø)
core	72.46% <70.00%> (-20.61%)	⬇️
files	91.69% <ø> (ø)
komamri	93.98% <100.00%> (ø)
plots	87.64% <50.00%> (-0.79%)	⬇️

Flags with carried forward coverage won't be shown. Click here to find out more.

Files with missing lines	Coverage Δ
KomaMRIBase/src/datatypes/Phantom.jl	84.29% <ø> (ø)
KomaMRICore/src/datatypes/Spinor.jl	88.23% <ø> (ø)
...RICore/src/simulation/SimMethods/Bloch/BlochCPU.jl	100.00% <100.00%> (ø)
...c/simulation/SimMethods/BlochSimple/BlochSimple.jl	100.00% <100.00%> (ø)
src/KomaUI.jl	99.27% <100.00%> (ø)
src/ui/ExportMATFunctions.jl	96.66% <ø> (ø)
...RICore/src/simulation/SimMethods/Bloch/BlochGPU.jl	0.00% <0.00%> (-100.00%)	⬇️
KomaMRIPlots/src/ui/DisplayFunctions.jl	88.26% <50.00%> (-0.86%)	⬇️

... and 8 files with indirect coverage changes

[cncastillo]

Wow, congratulations! this is a very good first PR. 🥳 How was the experience of generating a PR? where the instructions clear? any other comment that is not in #382?

EDIT: Forgot to ask, but did you try the new Bloch() simulation method? it should be a lot faster!

I have a few concerns, but it is better to support it like it is now at the start, and change it later, as more user could find it useful.

The concerns are: (1) if this makes it considerably less efficient if the user are not interested in B1+ (we can run the benchmarks for that), (2) problems with motion, for example for flowing spins, by putting B1+ in the Phantom it would make the B1+ "follow the spin", where I think it makes more sense to keep "B1 static" as it is more of a property of the scanner coils. Also, (3) that to keep adding B1+ channels to the Phantom may not be the best for parallel transmit (pTx). If you could put an issue with these points so I remember, that would be great :)

Do you think that plotting B1 (or any other complex property) as colors for the phase, and intensity of the color for magnitude (like optical flow) would be useful?

Here the colors represent displacement direction, but we could encode the phase.

Or maybe add the ability to see both magnitude and phase.

[JanWP]

Thanks for your quick reply! Regarding your concerns:

1. I was worried about efficiency as well, but since the existing code aggregates RF and local Bz into a spinor, there is already one spinor per spin and per time point, so this doesn't change much. I have run the benchmarks here locally and did not see any difference (<1%).
2. Yes, this was an issue that was a bit puzzling to me at first, also for the other phantom properties. When is a phantom property determined by the spin or the parts of the sample that move with the spin, and when is it an extraneous property? I don't think there is a general answer. Personally, I am interested in B1 maps to simulate the effects of RF currents in implants. In that case, if the patient moves, the implant moves as well, so it makes sense to have B1 follow the spins. But I realize that this is a pretty special case, and generally B1 is a composite effect of RF coils and dielectric properties of the sample. If the sample moves in that case, all bets are off as to what B1 should look like. I think the same is true for B0 heterogeneity from sources other than the parts of the sample that move with the spin. In that case d\omega shouldn't be a spin property.
3. Fully agree about pTx. I think there is room for both, an effect of the sample, and an effect from the coils.

As for plotting complex values with intensity + color, I debated trying that. My past attempts for this have given results that were rather underwhelming. It becomes really tricky to get the scale for saturation or intensity ranges correct if the output is going to be really informative, so I didn't bother for this first attempt and was more thinking on how to add both magnitude and phase separately (like for reconstructed images). The images you are showing are very nice though, so it might be useful to try again.

Not sure why the benchmarks fail, they worked on my machine, at least for CPU and CUDA.

[JanWP]

With the new Bloch() method, you mean BlochCPU.jl and BlochGPU.jl? Yes, I added B1+ to both of them. Unless there was something I missed in the merge from new-motion yesterday? But I didn't compare speed since the merge. Are both GPU and CPU accelerated? It was fast already!

[cncastillo]

Yes! The new methods should be at least 4x faster, and I have seen a speed-up of 100-200x for some sequences (the older one is now called BlochSimple()).

Thanks for the insight. It is tricky to say what moves with the object and what doesn't; maybe as you say, we can have one B1 component that belongs to the phantom and another related to the coils, and similarly for B0.

Don't worry about the benchmarks failing. That is just a permission issue as you are not yet an "official collaborator" and are doing a PR from a fork (and that is what is causing it to fail). I will run the benchmark locally when I have time, and merge it. Congratulations again on the nice work!

[cncastillo]

Hi, I will look into this later this week; as this is technically a breaking change, we can incorporate it and have a quick v1.10 version.

Answering some questions in the comments:

s_α_i[i_l] = -(Bz[i_g, t] / (B[i_g, t] + eps(T))) * sin_φ # why not zero?

That check is just an artifact of how the RF waveforms are considered,

  x-x-x
x/     \x

a square gradient would be sampled like the "figure above," so at the start and end (or even at the middle depending on the waveform), there could be one sample with B1 = 0. If more than one sample is 0, that goes to a precession block.

With #458, that is no longer needed as we do a trapezoidal integration B1_trap = (B1_n + B1_n+1)/2 , meaning we would never obtain B1_trap = 0. So the check will be removed in future releases.

... # TODO should this be more explicit on which symbols to include rather than using indices?

Good catch. This is super hardcoded. We now have a variable KomaMRIBase.VECTOR_PHANTOM_FIELDS that should do the trick.

[JanWP]

Thanks for the explanations! As far as I am concerned, there is no rush to merge this, in case you prefer to wait for other features to be ready before tagging the next release.