B0 mapping based on LLSf

[evaalonsoortiz]

I think this PR is ready for review/feedback.

Both field maps shown below are generated using high-resolution multi-echo (4 echoes) gradient echo data obtained at 7T, which can be found here: https://mfr.ca-1.osf.io/render?url=https://osf.io/s8qhk/?direct%26mode=render%26action=download%26mode=render

B0 map obtained from a linear fit of the phase evolution vs. time:

B0 map obtained using a the 4-quadrant arctangent of the first 2 echoes:

The windowing is the same for both field maps. No spatial unwrapping was used. Temporal unwrapping was used in the first case, but the smallest deltaTE is 2.7 ms so multiple phase wraps per echo spacing is expected in high deltaB0 regions.

[kousu]

I'm having trouble reproducing your example. So far I've put together what I know (and what you told me in the chatroom) into this script:

$ cat example/llsf.m 
% this should be run from inside the examples/ folder
addpath(genpath('..'));

% Relative path for the data
data = '7T_fieldmaps/';

% --------------------------------------
%% Download data when not already present
if ~isfolder( data )
    url = 'https://osf.io/s8qhk/download?version=1' ;
    fprintf( ['\n Downloading test data...\n URL=' url '\n'] ) ;
    unzip(url, '.') ;
end

% convert to nifti
% TODO: use imutils.dicom_to_nifti
!mkdir 7T_fieldmaps/ASPIRE_P_A_GRE_JASON_ASPIRE_100SLI_2MMISO_0006/nifti
!dcm2niix -o 7T_fieldmaps/ASPIRE_P_A_GRE_JASON_ASPIRE_100SLI_2MMISO_0006/nifti 7T_fieldmaps/ASPIRE_P_A_GRE_JASON_ASPIRE_100SLI_2MMISO_0006/
!mkdir 7T_fieldmaps/ASPIRE_M_A_GRE_JASON_ASPIRE_100SLI_2MMISO_0007/nifti
!dcm2niix -o 7T_fieldmaps/ASPIRE_M_A_GRE_JASON_ASPIRE_100SLI_2MMISO_0007/nifti 7T_fieldmaps/ASPIRE_M_A_GRE_JASON_ASPIRE_100SLI_2MMISO_0007/

[mag_data, mag_info, mag_json] = imutils.load_niftis('7T_fieldmaps/ASPIRE_M_A_GRE_JASON_ASPIRE_100SLI_2MMISO_0007/nifti');
[ph_data, ph_info, ph_json] = imutils.load_niftis('7T_fieldmaps/ASPIRE_P_A_GRE_JASON_ASPIRE_100SLI_2MMISO_0006/nifti');

mag_data = double(squeeze(mag_data));
ph_data = double(squeeze(ph_data));

ph_data = rescalePhaseImage(ph_data);

compl_vol(:,:,:,:) = mag_data(:,:,:,:).*exp(1i*ph_data(:,:,:,:));

[delf, offset, STDX, MSE] = B0_multiecho_linfit(compl_vol, ph_json)

which I run like this (on a unix machine):

$ (cd example; matlab -r "run('llsf.m');")

This runs, and I'm dropped to a matlab shell with

>> who

Your variables are:

MSE        STDX       compl_vol  data       delf       mag_data   mag_info   mag_json   offset     ph_data    ph_info    ph_json    
>> size(delf) 

ans =

   112   112   100

>> size(offset)

ans =

   112   112   100

>> size(STDX)

ans =

   112   112   100     2

>> size(MSE)

ans =

   112   112   100

I don't know what to do next though. How did you make those images? Did you have to do further processing?

This script could become a unit test, but before we do that I want to figure out how to treat the sample data as a fixture; I know @agahkarakuzu has advice on testing in https://github.com/agahkarakuzu/eda_organized, so maybe there's advice in there.

[evaalonsoortiz]

@kousu You could have a look at delf (the frequency map) using imagesc in matlab or save it as a nifti using niftiwrite(delf,'delf_map.nii') and then view it using fsleyes.

I will add an example script. One thing I wanted to consult with everyone about is where to keep the code. For now it is in /misc, but I'm not sure if that's where it should stay.