Create numerical model data for multi-echo field maps

[mathieuboudreau]

Resolves #102

I think I've done enough here to warrant a sanity check to make sure I'm going in the right track. for this hackathon project. Likely not actually ready for a full PR review, just wanted to start a discussion here.

• I've created a class to for the simulations, NumericalModel
• If initialized with a string parameter, it will define the volume geometry associated to that key.
  • So far, I've only implemented a 2D Shepp-Logan phantom, with key Shepp-Logan
  • The phantom associates T2star and proton density values. I've defined some default ones for now.

• By default, the background delta B0 is set to 0, but you can define one with the method generate_deltaB0 after creating an instance of the class.
  • So far, I've only created a linear field in the X direction (mx+b)
  • Example usage: obj_name.generate_deltaB0('linear', [m, b]) where m and b are single values.

• After doing so, you can generate a 4D measurement volume (4th volume is different TEs)
  • These are treated as separate/independent measurements in the signal model.
  • Can TE can be an array
  • You can optionally define SNR to add noise during the method call.
  • There are some get functions to get the magnitude/phase/real/imaginary volumes of your choice.

Example call:

test_obj=NumericalModel('Shepp-Logan');
test_obj.generate_deltaB0('linear', [0.05 0.01]);
test_obj.simulate_measurement(15, [0.001 0.015], 10);
magn = test_obj.getMagnitude;
phase = test_obj.getPhase;
figure(),imagesc(magn(:,:,1,1))
figure(),imagesc(phase(:,:,1,1))

• You can export your simulated volume to NIFTI using the save method
  • Example call: test_obj.save('Magnitude', 'simulated_magnitude.nii');

Opened in FSLeyes:

Lastly, for all the code I wrote, I wrote tests for them. Since I didn't see many tests to base myself off of in your master branch, I used what I felt comfortable with, which is MATLAB's class-based unit test and tags workflow. You can run them by (after adding all the repo's folders to your MATLAB path) opening the tests directory and running result = runtests(pwd, 'Recursively', true).

>> result = runtests(pwd, 'Recursively', true)
Running NumericalModel_Test
.......... ....
Done NumericalModel_Test
__________

result = 

  1×14 TestResult array with properties:

    Name
    Passed
    Failed
    Incomplete
    Duration
    Details

Totals:
   14 Passed, 0 Failed, 0 Incomplete.
   0.42035 seconds testing time.

>> 

[po09i]

This looks really good to me! Nice format and I like the test suite!

Some ideas (not suggestions) I had when looking at you PR, not sure all are relevant

• We could input any grayscale image to get a resulting model
  • I see issues with the different MR values (T2*, PD, ...)
• Save the TE values as a property to refer to them later
• When saving the nifti, output a json sidecar to store relevant infos such as TEs
• Add possibility to input T1? (not sure if it would change much)
• Add a 5th dimension (acquisition number)
  • These could originally just change the noise from frame to frame but potentially get more realistic with slight movement of the "subject" , ...)

@rtopfer What do you think? any other ideas?

[mathieuboudreau]

Thanks for the insights @po09i !

Maybe as another sanity check to make sure I'm simulating correctly, I could add a test to make sure that a fitted B0 map from this data gives approximately the input B0 map. Is there a function/API already set up in the toolbox such that I could easily fit my data in a test? Dual-echo should be fine for this sanity check.

[jcohenadad]

Thanks for the insights @po09i !

Maybe as another sanity check to make sure I'm simulating correctly, I could add a test to make sure that a fitted B0 map from this data gives approximately the input B0 map. Is there a function/API already set up in the toolbox such that I could easily fit my data in a test? Dual-echo should be fine for this sanity check.

@mathieuboudreau we are in the middle of a major refactoring, and i am not sure pointing you to the currently-implemented fieldmap function is relevant (will break soon). So for now, i suggest you use a "hard-coded" fieldmap computation as we do in the up-to-date example file.

Note: this code relies on BIDS data (with json file that has TE info), but you can easily bypass that if you need. Ultimately the fieldmap API will be agnostic to BIDS and we will just pass the necessary argument (matrix, TE)

[mathieuboudreau]

Thanks for the insights @po09i ! Maybe as another sanity check to make sure I'm simulating correctly, I could add a test to make sure that a fitted B0 map from this data gives approximately the input B0 map. Is there a function/API already set up in the toolbox such that I could easily fit my data in a test? Dual-echo should be fine for this sanity check.

@mathieuboudreau we are in the middle of a major refactoring, and i am not sure pointing you to the currently-implemented fieldmap function is relevant (will break soon). So for now, i suggest you use a "hard-coded" fieldmap computation as we do in the up-to-date example file.

Note: this code relies on BIDS data (with json file that has TE info), but you can easily bypass that if you need. Ultimately the fieldmap API will be agnostic to BIDS and we will just pass the necessary argument (matrix, TE)

Ok thanks @jcohenadad ! I just threw the initial question out there in case there was already something in the software - I didn't want to duplicate code. I've been testing it on my side using the dual-angle B0 map equation I'm familiar with too, which is that one.

[rtopfer]

@mathieuboudreau what i just posted on slack may be of interest : links here for posterity : https://www.mathworks.com/matlabcentral/fileexchange/30853-field-mapping-toolbox https://www.mathworks.com/matlabcentral/fileexchange/22508-siemens-dicom-sort-and-convert-to-nifti https://github.com/shimming-toolbox/shimming-toolbox/blob/f320602c07c97216176daaee3ebf8ca0542368da/Img/FieldEval.m#L744

[mathieuboudreau]

Thanks for the insights everyone!

My "sanity check" just raised an issue, possibly related to my sanity, but likely to the physics of this problem.

When I calculated a B0 from simulated measurements using the equations shown above in this file, the resulting B0 map is off by a sign (+/-).

>> a=NumericalModel('Shepp-Logan');
>> B0_hz = 13;
>> a.generate_deltaB0('linear', [0.0 B0_hz]); % constant B0 of 10 Hz. deltaB0 is in Tesla.
>> TR = 0.025;
>> TE = [0.004 0.008];
>> a.simulate_measurement(TR, TE);
>> phaseMeas = a.getPhase();
>> phaseTE1 = squeeze(phaseMeas(:,:,1,1));
>> phaseTE2 = squeeze(phaseMeas(:,:,1,2));
>> B0_meas = (phaseTE2(64, 64) - phaseTE1(64, 64))/(TE(2) - TE(1))

B0_meas =

  -81.6814

>> B0_meas_hz = B0_meas/(2*pi)

B0_meas_hz =

  -13.0000

>> 

Note that this this is my simulation of the signal (taken from this paper, referenced by @rtopfer in this comment:

https://github.com/shimming-toolbox/shimming-toolbox/blob/c2bf18f649bb47d43c91c6ed59fb3e6691ff0798/numerical_model/NumericalModel.m#L159

Now here's where I need a sanity check. Although most textbooks show the Larmor equation in this form (taken from Nishimura):

which would lead some to believe that (since omega = dPhase/dt), you should calculate a dual angle B0 map from B0 = (Phase2 - Phase1)/(TE2-TE1) (divided by 2Pi to get linear Hz instead of rad*Hz). However, omega and B0 are actually vectors, and in vector form, the have opposite directions. From Haacke's green bible:

So the direction of change in phase is actually in the opposite direction from B0. So should the equation to calculate B0 actually be B0 = -(Phase2 - Phase1)/(TE2-TE1) (or in a more confusion notation, B0 = (Phase1 - Phase2)/(TE2-TE1) )?

@evaalonsoortiz @jcohenadad @po09i ?

[rtopfer]

So the direction of change in phase is actually in the opposite direction from B0. So should the equation to calculate B0 actually be B0 = -(Phase2 - Phase1)/(TE2-TE1) (or in a more confusion notation, B0 = (Phase1 - Phase2)/(TE2-TE1) )?

Not sure but the answer might also be mentioned in passing somewhere in those textbooks, in any case, it depends on the "handedness" convention the MRI manufacturer adopts for phase rotation. We've only used Siemens images in the past ("left-handed") but I believe the sign would be flipped for GE systems, eg see

[mathieuboudreau]

So the direction of change in phase is actually in the opposite direction from B0. So should the equation to calculate B0 actually be B0 = -(Phase2 - Phase1)/(TE2-TE1) (or in a more confusion notation, B0 = (Phase1 - Phase2)/(TE2-TE1) )?

Not sure but the answer might also be mentioned in passing somewhere in those textbooks, in any case, it depends on the "handedness" convention the MRI manufacturer adopts for phase rotation. We've only used Siemens images in the past ("left-handed") but I believe the sign would be flipped for GE systems, eg see

Ok that make sense - then maybe I should be adding a parameter in the simulation to determine the handedness (done by changing the sign in the complex exponential), and set the default one to whichever your B0 mapping algorithms are defined for.

[mathieuboudreau]

I've defined a class attribute called handedness with a default of 'left' (for Siemens), which flips the sign in the complex exponential of the simulation. This should make the simulations compatible with your B0 fitting algorithms. It might be a good idea to provide that same interface option to your B0 mapping algorithms, so that people can use them correctly between vendors.

>> a=NumericalModel('Shepp-Logan');
>> B0_hz = 13;
>> a.generate_deltaB0('linear', [0.0 B0_hz]); % constant B0 of 10 Hz. deltaB0 is in Tesla.
>> TR = 0.025;
>> TE = [0.004 0.008];
>> a.simulate_measurement(TR, TE);
>> phaseMeas = a.getPhase();
>> phaseTE1 = squeeze(phaseMeas(:,:,1,1));
>> phaseTE2 = squeeze(phaseMeas(:,:,1,2));
>> B0_meas = (phaseTE2(64, 64) - phaseTE1(64, 64))/(TE(2) - TE(1))

B0_meas =

   81.6814

>> B0_meas_hz = B0_meas/(2*pi)

B0_meas_hz =

   13.0000

And now switching the handedness will switch the sign of the B0 with the same dual-echo B0 calculation:

>> a.handedness = 'right'

a = 

  NumericalModel with properties:

              gamma: 2.6752e+08
      fieldStrength: 3
         handedness: 'right'
             numVox: 128
    starting_volume: [128×128 double]
             volume: [1×1 struct]
        measurement: [128×128×1×2 double]
             T2star: [1×1 struct]
      protonDensity: [1×1 struct]
            deltaB0: [128×128 double]

>> a.simulate_measurement(TR, TE);
>> phaseMeas = a.getPhase();
>> phaseTE1 = squeeze(phaseMeas(:,:,1,1));
>> phaseTE2 = squeeze(phaseMeas(:,:,1,2));
>> B0_meas = (phaseTE2(64, 64) - phaseTE1(64, 64))/(TE(2) - TE(1))

B0_meas =

  -81.6814

>> B0_meas_hz = B0_meas/(2*pi)

B0_meas_hz =

  -13.0000

[evaalonsoortiz]

Indeed, Siemens uses the left-handed coordinate system . Left vs. right coordinate systems flip the direction of the z-axis. So in the left-handed coordinate system a clockwise rotation about z has a negative rate of change of theta (as shown in Eq. 2.28 from Haacke). The handedness attribute is perfect. That would change the sign in front of deltaB0 here, correct?

[mathieuboudreau]

Indeed, Siemens uses the left-handed coordinate system . Left vs. right coordinate systems flip the direction of the z-axis. So in the left-handed coordinate system a clockwise rotation about z has a negative rate of change of theta (as shown in Eq. 2.28 from Haacke). The handedness attribute is perfect. That would change the sign in front of deltaB0 here, correct?

Great! Yup, see here:

https://github.com/shimming-toolbox/shimming-toolbox/blob/e9af472c529b7f756783a7a9b77657e56968daed/numerical_model/NumericalModel.m#L169

[mathieuboudreau]

Now that this simulation class seems to be working as (I) intended, I played around with it a bit and generated a couple of gifs using it:

Linear B0 field, high-ish SNR:

Script:

```MATLAB clear all clc close all a=NumericalModel('Shepp-Logan'); B0_hz = 13; a.generate_deltaB0('linear', [0.1 B0_hz]); % constant B0 of 10 Hz. deltaB0 is in Tesla. TR = 0.025; TE = [0.001:0.001:0.050]; SNR = 100; h=figure() figTitle = sgtitle(['']); filename = 'output2.gif' for ii = 1:length(TE) figTitle.String = ['Echo time = ', num2str(TE(ii)*1000), ' ms']; a.simulate_measurement(TR, TE(ii), SNR); magMeas = a.getMagnitude(); phaseMeas = a.getPhase(); subplot(2,2,1); imagesc(magMeas(:,:,1)), colorbar, caxis([0 0.03]), axis image title('Magnitude') axis off subplot(2,2,2); imagesc(phaseMeas(:,:,1)), colorbar, caxis([-pi pi]), axis image title('Phase') axis off subplot(2,2,3); imagesc(a.deltaB0*(a.gamma/(2*pi))), colorbar, caxis([0 30]), axis image title('Actual deltaB0 (Hz)') axis off phaseMeas = a.getPhase(); if ii == 1 refPhase = phaseMeas; refTE = TE(ii); else subplot(2,2,4); currentPhase = phaseMeas; currentTE = TE(ii); B0_meas = (currentPhase- refPhase)./(currentTE - refTE); B0_meas_hz = B0_meas/(2*pi); imagesc(B0_meas_hz), colorbar, caxis([0 30]), axis image title('B0 map (current TE relative to initial TE) (Hz)') axis off end drawnow % Capture the plot as an image frame = getframe(h); im = frame2im(frame); [imind,cm] = rgb2ind(im,256); % Write to the GIF File if ii == 1 imwrite(imind,cm,filename,'gif', 'Loopcount',inf); else imwrite(imind,cm,filename,'gif','DelayTime',0.2,'WriteMode','append'); end end ```

Constant B0 field, low SNR:

Script:

```MATLAB clear all clc close all a=NumericalModel('Shepp-Logan'); B0_hz = 13; a.generate_deltaB0('linear', [0 B0_hz]); % constant B0 of 10 Hz. deltaB0 is in Tesla. TR = 0.025; TE = [0.001:0.001:0.050]; SNR = 20; h=figure() figTitle = sgtitle(['']); filename = 'output.gif' for ii = 1:length(TE) figTitle.String = ['Echo time = ', num2str(TE(ii)*1000), ' ms']; a.simulate_measurement(TR, TE(ii), SNR); magMeas = a.getMagnitude(); phaseMeas = a.getPhase(); subplot(2,2,1); imagesc(magMeas(:,:,1)), colorbar, caxis([0 0.03]), axis image title('Magnitude') axis off subplot(2,2,2); imagesc(phaseMeas(:,:,1)), colorbar, caxis([-pi pi]), axis image title('Phase') axis off subplot(2,2,3); imagesc(a.deltaB0*(a.gamma/(2*pi))), colorbar, caxis([0 30]), axis image title('Actual deltaB0 (Hz)') axis off phaseMeas = a.getPhase(); if ii == 1 refPhase = phaseMeas; refTE = TE(ii); else subplot(2,2,4); currentPhase = phaseMeas; currentTE = TE(ii); B0_meas = (currentPhase- refPhase)./(currentTE - refTE); B0_meas_hz = B0_meas/(2*pi); imagesc(B0_meas_hz), colorbar, caxis([0 30]), axis image title('B0 map (current TE relative to initial TE) (Hz)') axis off end drawnow % Capture the plot as an image frame = getframe(h); im = frame2im(frame); [imind,cm] = rgb2ind(im,256); % Write to the GIF File if ii == 1 imwrite(imind,cm,filename,'gif', 'Loopcount',inf); else imwrite(imind,cm,filename,'gif','DelayTime',0.2,'WriteMode','append'); end end ```

I didn't do any unwrapping, since I know that another part of the work currently going on, so you can see at which point the B0 map breaks down if no unwrapping is done.

[evaalonsoortiz]

Very nice!

[rtopfer]

@mathieuboudreau that's so cool! going to make a point of copying your .gif plot routine whenever/wherever possible 🙏

[mathieuboudreau]

At your convenience, I think this PR is ready for review as it resolves all the points listed in the issue #102 by @po09i and myself. No rush, I know you folks are busy with the other parts of the hackathon. Cleanup is likely required to adhere to your coding conventions. The test fails because I have MATLAB-based tests, but I could convert these to Octave-compatible if you'd prefer.

[evaalonsoortiz]

I am currently trying to test your code by running using my B0 mapping code on data simulated with your code. One thing I need to be able to do is load nifti header or json sidecar information that contains the echo times.

I am trying to figure out how to save the echo times used in the simulation in the nifti header for the files I save. For now the (messy) solution I've come up with is

• saving the data in a nifti file using test_obj.save('Magnitude', 'simulated_magnitude.nii');

• then loading the nifti header using matlab's info = niftiinfo('simulated_magnitude.nii');

• adding the echo times to the Description field

• then saving the data and updated header with niftiwrite(magn,'simulated_magnitude.nii',info);

There's got to be a cleaner way, ie, generating a nifti header structure when initially creating the nii file. Does anyone know of a way to do this?

[jcohenadad]

There's got to be a cleaner way, ie, generating a nifti header structure when initially creating the nii file. Does anyone know of a way to do this?

IIUC, the problem is that @mathieuboudreau's output from the simulation does not have a json sidecar, right? if that's the case, then i would fix the problem at the source, i.e. write a json using matlab's jsonencode

[mathieuboudreau]

I am currently trying to test your code by running using my B0 mapping code on data simulated with your code. One thing I need to be able to do is load nifti header or json sidecar information that contains the echo times.

I am trying to figure out how to save the echo times used in the simulation in the nifti header for the files I save. For now the (messy) solution I've come up with is

• saving the data in a nifti file using test_obj.save('Magnitude', 'simulated_magnitude.nii');
• then loading the nifti header using matlab's info = niftiinfo('simulated_magnitude.nii');
• adding the echo times to the Description field
• then saving the data and updated header with niftiwrite(magn,'simulated_magnitude.nii',info);

There's got to be a cleaner way, ie, generating a nifti header structure when initially creating the nii file. Does anyone know of a way to do this?

We could probably save the TEs to a class attribute, and then set them after making the nifti file variable from the empty template here:

https://github.com/shimming-toolbox/shimming-toolbox/pull/127/files#diff-2ce54f71af5871173fada60703cb0ecbR129-R131

i.e. edit the Description field like you suggest, butin the nii_vol variable after creating it in line 130 above.

We could export a JSON sidecar in the same save routine as well, as Julien mentioned.

[mathieuboudreau]

@evaalonsoortiz JSON sidecars are now created with the TE and FA values (EchoTime and FlipAngle). I didn't add them to the NIfTI header, since I looked it up and it appears these are never/rarely saved in NIfTI headers (I was getting mixed up above because I'm used to MINC, which do save them). But if you want them in NIfTI's too, just let me know where in the header and what name you'd like them to be and I can save them too.

[jcohenadad]

No need to encode that in the nifti, that’s what the json is for. Thank you Mathieu!!

[evaalonsoortiz]

This is great Mathieu! I've tested it out with my multi-echo fit. Here are my results:

Code

``` test_obj=NumericalModel('Shepp-Logan'); test_obj.generate_deltaB0('linear', [10 0]); test_obj.simulate_measurement(15, [0.0015 0.002 0.0025 0.003 0.0035 0.004 0.0045 0.005 0.0055 0.006 0.0065 0.007], 100); magn = test_obj.getMagnitude; phase = test_obj.getPhase; test_obj.save('Phase', 'simulated_phase.nii'); test_obj.save('Magnitude', 'simulated_magnitude.nii'); compl_vol(:,:,:,:) = magn(:,:,:,:).*exp(1i*phase(:,:,:,:)); ph_json = jsondecode(fileread('simulated_phase.nii.json')); [delf] = B0_multiecho_linfit(compl_vol, ph_json);

[po09i]

This is really nice @mathieuboudreau

[evaalonsoortiz]

@mathieuboudreau Just one thing I noticed: the .json file is saved with a double extension, ex: filename.nii.json

Can this be fixed?

[mathieuboudreau]

This is great Mathieu! I've tested it out with my multi-echo fit. Here are my results:

Thanks! Awesome!

@mathieuboudreau Just one thing I noticed: the .json file is saved with a double extension, ex: filename.nii.json

Can this be fixed?

Sorry about that @evaalonsoortiz ! I just pushed changes + more tests to my branch which should resolve this issue / your use-case.