Allow separate B1 map for MTsat B1 correction

CHANGELOG.md

@@ -7,29 +7,29 @@ Most recent version numbers *should* follow the [Semantic Versioning](https://se
 ## [unreleased]
 ### Added
 - option to choose different models and parameters for B1-correction of MTsat
-- set default WM percent value in hmri_defaults
+- set default WM percent value in `hmri_defaults`
 - spatial processing: add explicit mask creation and fix implicit mask (0 to NaN in float images)
 - update FIL seste seq parameters in get_metadata_val_classic
 - denoising module-first part: Java-Matlab interface for LCPCA denoising
 - read EffectiveEchoTime in new TerraX Dicom format
 - save LCPCA-denoising supplementary files as nifti instead of .mat
 - parameter error maps
 - robust combination of two runs using error maps
-
+- option to use a separate B1 map for B1 correction of MTsat; useful if pTx used for excitation pulses and CP mode for the MT pulse
 
 ### Fixed
 - replace `datestr(now)` with `datetime('now')` in line with [MATLAB recommendation](https://mathworks.com/help/matlab/matlab_prog/replace-discouraged-instances-of-serial-date-numbers-and-date-strings.html)
 - fix crash if input images have different matrix sizes, and warn
 - make B1-map creation using 3DEPI SE/STE and AFI methods fall back to defaults without sidecar files, rather than crash
-- Modify the filenames as files are copied to RFsensCalc to prevent overwriting in further processing
+- modify the filenames as files are copied to RFsensCalc to prevent overwriting in further processing
 - batch interface now enforces the number of B1 input images correctly for B1 mapping methods which only need two images
 - more informative error if optimization toolbox not present during NLLS R2* calculation
 - fix 3D-EPI B1 mapping not using b1defaults for Triotim scanner
-- use cell- instead of char- array to accommodate filenames of unequal length in RFsens
+- use cell- instead of char-array to accommodate filenames of unequal length in RFsens
 
 ## [v0.6.1]
 ### Fixed
-- The local config files have been converted to scripts for compatibility with compiled version
+- local config files have been converted to scripts for compatibility with compiled version
 - function-evaluate SPM-struct (preproc8.val) for SPM development version compatibility
 - copy acquisition metadata to TE=0 volumes in Results/Supplementary folder after map creation so they can be used as input to the toolbox if needed
 

hmri_create_MTProt.m

@@ -99,6 +99,13 @@
 % P_trans(1,:) = magnitude image (anatomical reference for coregistration)
 % P_trans(2,:) = B1 map (p.u.)
 P_trans = jobsubj.b1_trans_input;
+if isfield(jobsubj,'b1_MT_trans_input')
+    P_trans_MT = jobsubj.b1_MT_trans_input;
+    separate_MT_B1map = true;
+else
+    P_trans_MT = [];
+    separate_MT_B1map = false;
+end
 
 % index number for each contrast - zero index means no images available
 PDwidx = mpm_params.PDwidx;
@@ -279,6 +286,17 @@
     V_trans = spm_vol(P_trans);
 end
 
+V_trans_MT = [];
+if ~isempty(P_trans_MT)
+    % Load B1 mapping data if available and coregister to PDw
+    % P_trans(1,:) = magnitude image (anatomical reference for coregistration) 
+    % P_trans(2,:) = B1 map (p.u.)
+    if mpm_params.coreg
+        hmri_coreg(Pavg{PDwidx}, P_trans_MT, mpm_params.coreg_bias_flags);
+    end
+    V_trans_MT = spm_vol(P_trans_MT);
+end
+
 % parameters saved for quality assessment
 if mpm_params.QA.enable
     if exist(mpm_params.QA.fnam,'file')
@@ -826,13 +844,27 @@
                 mpm_params.small_angle_approx);
             MT = hmri_calc_MTsat(struct('data',MTw,'fa',fa_mtw_rad,'TR',TR_mtw,'B1',B1_mtw), A_forMT, R1_forMT);
 
+            % Use MT-pulse B1 map if available
+            if ~isempty(V_trans_MT) % separate B1 mapping data provided for MT pulse
+                f_T_MT = hmri_read_vols(V_trans_MT(2,:),Ni,p,mpm_params.interp)/100; % divide by 100, since p.u. maps
+                switch mpm_params.MTsatB1CorrectionModel
+                    % need to undo implicit B1^2 scaling from MPM B1 map
+                    case 'helms'
+                        MT = MT.*(f_T./f_T_MT).^2;
+                end
+            elseif separate_MT_B1map % separate B1 mapping data for MT pulse selected, but with "no B1 correction" selected
+                f_T_MT = [];
+            else % same B1 mapping data for excitation and MT pulse selected
+                f_T_MT = f_T;
+            end
+
             % f_T correction is applied either if:
             % - f_T has been provided as separate B1 mapping measurement (not
             % UNICORT!) or
             % - f_T has been calculated using UNICORT *AND* the UNICORT.MT flag
             % is enabled (advanced user only! method not validated yet!)
-            if (~isempty(f_T))&&(~mpm_params.UNICORT.R1 || mpm_params.UNICORT.MT)
-                [MT,dMTsatcorr] = hmri_correct_MTsat(MT,f_T,mpm_params.MTsatB1CorrectionModel,mpm_params.MTsatB1CorrectionC);
+            if (~isempty(f_T_MT))&&(~mpm_params.UNICORT.R1 || mpm_params.UNICORT.MT)
+                [MT,dMTsatcorr] = hmri_correct_MTsat(MT,f_T_MT,mpm_params.MTsatB1CorrectionModel,mpm_params.MTsatB1CorrectionC);
             end
 
             % truncate MT maps
@@ -844,7 +876,7 @@
 
                 dMT = hmri_calc_dMT(PDw,T1w,MTw,Edata.PDw,Edata.T1w,Edata.MTw,fa_pdw_rad,fa_t1w_rad,fa_mtw_rad,TR_pdw,TR_t1w,TR_mtw,mpm_params.small_angle_approx) * 100;
 
-                if (~isempty(f_T))&&(~mpm_params.UNICORT.R1 || mpm_params.UNICORT.MT)
+                if (~isempty(f_T_MT))&&(~mpm_params.UNICORT.R1 || mpm_params.UNICORT.MT)
                     % use chain rule to include MTsat B1 correction in error maps
                     dMT = dMT.*abs(dMTsatcorr);
                 end
@@ -1291,7 +1323,12 @@ function PDcalculation(fA, fTPM, mpm_params)
 end
 
 % UNICORT settings:
-mpm_params.UNICORT.R1 = isfield(jobsubj.b1_type,'UNICORT'); % uses UNICORT to estimate B1 transmit bias
+if isfield(jobsubj.b1_type,'b1_MT')
+    b1_type = jobsubj.b1_type.b1_MT.b1_type;
+else
+    b1_type = jobsubj.b1_type;
+end
+mpm_params.UNICORT.R1 = isfield(b1_type,'UNICORT'); % uses UNICORT to estimate B1 transmit bias
 tmp = hmri_get_defaults('UNICORT');
 mpm_params.UNICORT.PD = tmp.PD; % uses B1map estimated as biasfield for R1 to correct for B1 transmit bias in PD
 mpm_params.UNICORT.MT = tmp.MT; % uses B1map estimated as biasfield for R1 to correct for B1 transmit bias in MT
@@ -1642,7 +1679,7 @@ function PDcalculation(fA, fTPM, mpm_params)
 
 % Prepare output for R1, PD, MT and R2* maps
 RFsenscorr = fieldnames(mpm_params.proc.RFsenscorr);
-B1transcorr = fieldnames(jobsubj.b1_type);
+B1transcorr = fieldnames(b1_type);
 coutput = 0;
 
 % R1 output: requires PDw and T1w input and optional B1 transmit (or

hmri_create_b1map.m

@@ -97,9 +97,14 @@
 %   otherwise copyfile does not find the files!!
 
 if ~isempty(P_trans)
+    if ~isfield(jobsubj,'b1_suffix')
+        b1_suffix = '';
+    else
+        b1_suffix = jobsubj.b1_suffix;
+    end
     P_trans = spm_file(P_trans,'number','');
-    P_trans_copy{1} = fullfile(jobsubj.path.b1respath, spm_file(P_trans(1,:), 'filename'));
-    P_trans_copy{2} = fullfile(jobsubj.path.b1respath, spm_file(P_trans(2,:), 'filename'));
+    P_trans_copy{1} = fullfile(jobsubj.path.b1respath, spm_file(spm_file(P_trans(1,:), 'filename'),'suffix',b1_suffix));
+    P_trans_copy{2} = fullfile(jobsubj.path.b1respath, spm_file(spm_file(P_trans(2,:), 'filename'),'suffix',b1_suffix));
     copyfile(deblank(P_trans(1,:)), P_trans_copy{1});
     try copyfile([spm_str_manip(P_trans(1,:),'r') '.json'],[spm_str_manip(P_trans_copy{1},'r') '.json']); end %#ok<*TRYNC>
     copyfile(deblank(P_trans(2,:)), P_trans_copy{2});

hmri_run_create.m

@@ -88,13 +88,16 @@
 % define other (temporary) paths for processing data
 b1path = fullfile(outpath, 'B1mapCalc');
 if ~exist(b1path,'dir'); mkdir(b1path); end
+b1_MT_path = fullfile(outpath, 'B1mapCalc_MT');
+if ~exist(b1_MT_path,'dir'); mkdir(b1_MT_path); end
 rfsenspath = fullfile(outpath, 'RFsensCalc');
 if ~exist(rfsenspath,'dir'); mkdir(rfsenspath); end
 mpmpath = fullfile(outpath, 'MPMCalc');
 if ~exist(mpmpath,'dir'); mkdir(mpmpath); end
 
 % save all these paths in the job.subj structure
 job.subj.path.b1path = b1path;
+job.subj.path.b1_MT_path = b1_MT_path;
 job.subj.path.b1respath = supplpath; % copy B1 maps to supplementary folder
 job.subj.path.rfsenspath = rfsenspath;
 job.subj.path.mpmpath = mpmpath;
@@ -125,7 +128,25 @@
 spm_jsonwrite(fullfile(supplpath,'hMRI_map_creation_job_create_maps.json'),job,struct('indent','\t'));
 
 % run B1 map calculation for B1 bias correction
-P_trans = hmri_create_b1map(job.subj);
+if isfield(job.subj.b1_type,'b1_MT')
+    % excitation pulse B1 map
+    tmpsub = job.subj;
+    tmpsub.b1_suffix = '';
+    tmpsub.b1_type = job.subj.b1_type.b1_MT.b1_type;
+    job.subj.b1_trans_input = hmri_create_b1map(tmpsub);
+
+    % MT pulse B1 map
+    tmpsub = job.subj;
+    tmpsub.b1_suffix = '_MT';
+    tmpsub.path.b1path = b1_MT_path;
+    tmpsub.b1_type = job.subj.b1_type.b1_MT.b1_type_MT;
+    job.subj.b1_MT_trans_input = hmri_create_b1map(tmpsub);
+else
+    tmpsub = job.subj;
+    tmpsub.b1_suffix = '';
+    tmpsub.b1_type = job.subj.b1_type;
+    job.subj.b1_trans_input = hmri_create_b1map(tmpsub);
+end
 
 % check, if RF sensitivity profile was acquired and do the recalculation
 % accordingly
@@ -134,7 +155,6 @@
 end
 
 % run hmri_create_MTProt to evaluate the parameter maps
-job.subj.b1_trans_input = P_trans;
 [fR1, fR2s, fMT, fA, PPDw, PT1w, PMTw, Perror]  = hmri_create_MTProt(job.subj);
 
 % collect outputs:
@@ -166,6 +186,7 @@
 % clean after if required
 if hmri_get_defaults('cleanup')
     rmdir(job.subj.path.b1path,'s');
+    rmdir(job.subj.path.b1_MT_path,'s');
     rmdir(job.subj.path.rfsenspath,'s');
     rmdir(job.subj.path.mpmpath,'s');
 end
@@ -175,4 +196,4 @@
 
 hmri_log(sprintf('\t============ CREATE hMRI MAPS MODULE: completed (%s) ============', datetime('now')),flags);
 
-end
+end

tbx_scfg_hmri_create.m

@@ -78,7 +78,7 @@
 % ---------------------------------------------------------------------
 % menu b1_type
 % ---------------------------------------------------------------------
-[b1_type,b1parameters] = tbx_scfg_hmri_B1_menu;
+[b1_type_orig,b1parameters] = tbx_scfg_hmri_B1_menu;
 
 % ---------------------------------------------------------------------
 % New B1 mapping methods should be added to tbx_scfg_hmri_B1_menu.m
@@ -113,6 +113,7 @@
 % Add extra B1 mapping methods which cannot be run independently of the
 % MPM map creation to the menu
 % ---------------------------------------------------------------------
+b1_type = b1_type_orig;
 b1_type.values  = [b1_type.values {b1_input_UNICORT, b1_input_noB1}];
 b1_type.help=[b1_type.help; {...
     [' - UNICORT: Use this option when B1 maps not available. ' ...
@@ -123,6 +124,34 @@
     'Consider using UNICORT instead.']
     }];
 
+% ---------------------------------------------------------------------
+% Separate B1 map for MT pulse
+% ---------------------------------------------------------------------
+% potentially separate B1 map type for MT pulse (no UNICORT option)
+b1_type_MT = b1_type_orig;
+b1_type_MT.values  = [b1_type_MT.values {b1_input_noB1}];
+b1_type_MT.help=[b1_type_MT.help; {...
+    ' - no B1 correction: This option should only be used when there is no matching B1 mapping data. '}];
+b1_type_MT.name = 'MT pulse B1 mapping method';
+b1_type_MT.tag = 'b1_type_MT';
+
+b1_MT      = cfg_branch;
+b1_MT.tag  = 'b1_MT';
+b1_MT.name = 'Separate excitation and MT pulse B1 maps';
+b1_MT.help = {['One set of B1 mapping data is acquired for the excitation pulse ' ...
+    'and another for the MT pulse.']};
+b1_MT.val  = {setfield(b1_type,'name','Excitation pulse B1 mapping method') b1_type_MT}; %#ok<SFLD> 
+
+% ---------------------------------------------------------------------
+% Add separate B1 mapping for MT pulse to the menu
+% ---------------------------------------------------------------------
+b1_type.values  = [b1_type.values {b1_MT}];
+b1_type.help=[b1_type.help; {...
+    [' - separate excitation and MT pulse B1 maps: This option should be used ' ...
+     'when the excitation and MT pulses have different B1 inhomogeneities, e.g. ' ...
+     'due to the use of pTx ']
+    }];
+
 % ---------------------------------------------------------------------
 % Input images for RF sensitivity - RF sensitivity maps for MTw images
 % ---------------------------------------------------------------------