diff --git a/EPG_for_hmri_tbx/Unit_Test_Protocol.m b/EPG_for_hmri_tbx/Unit_Test_Protocol.m new file mode 100755 index 00000000..94cd91f8 --- /dev/null +++ b/EPG_for_hmri_tbx/Unit_Test_Protocol.m @@ -0,0 +1,8 @@ +global_hrmi_def +hmri_def.MPMacq_set.name = 'Unit_Test_Protocol'; +hmri_def.MPMacq_set.tag = 'Unit_Test_Protocol'; +hmri_def.MPMacq_set.val = [25 25 6 21]; +hmri_def.imperfectSpoilCorr.Unit_Test_Protocol.tag = 'Unit_Test_Protocol'; +hmri_def.imperfectSpoilCorr.Unit_Test_Protocol.P2_a = [41.5044 -56.2888 28.6704]; +hmri_def.imperfectSpoilCorr.Unit_Test_Protocol.P2_b = [-0.0972 0.0973 0.9649]; +hmri_def.imperfectSpoilCorr.Unit_Test_Protocol.enabled = hmri_def.imperfectSpoilCorr.enabled; diff --git a/config/hmri_b1_standard_defaults.m b/config/hmri_b1_standard_defaults.m index 3e66d692..46b5199d 100644 --- a/config/hmri_b1_standard_defaults.m +++ b/config/hmri_b1_standard_defaults.m @@ -47,6 +47,12 @@ hmri_def.b1map.i3D_AFI.b1acq.TR2TR1ratio = 5; hmri_def.b1map.i3D_AFI.b1acq.alphanom = 60; +% 'SDAM' +hmri_def.b1map.SDAM.b1type = 'SDAM'; +hmri_def.b1map.SDAM.b1avail = true; +hmri_def.b1map.SDAM.procreq = true; +hmri_def.b1map.SDAM.b1acq.alphanom = 60; + % 'pre_processed_B1' hmri_def.b1map.pre_processed_B1.b1type = 'pre_processed_B1'; hmri_def.b1map.pre_processed_B1.b1avail = true; @@ -100,4 +106,4 @@ hmri_def.b1map.rf_map.b1avail = true; hmri_def.b1map.rf_map.procreq = true; -end \ No newline at end of file +end diff --git a/config/hmri_defaults.m b/config/hmri_defaults.m index 67702158..d1147189 100644 --- a/config/hmri_defaults.m +++ b/config/hmri_defaults.m @@ -59,7 +59,7 @@ % text. The following settings are recommended. No modification currently % foreseen as useful... hmri_def.json = struct('extended',false,'separate',true,'anonym','none',... - 'overwrite',true, 'indent','\t'); + 'overwrite',true, 'indent','\t'); % recommended TPM for segmentation and spatial processing. The hMRI toolbox % provides a series of tissue probability maps. These TPMs could be % replaced by other TPMs, to better match the population studied. @@ -280,116 +280,50 @@ hmri_def.MPMacq.fa_pdw = 6; % [deg] hmri_def.MPMacq.tag = 'v2k'; -% IMPREFECT RF SPOILING CORRECTION PARAMETERS -% (Preibisch and Deichmann, MRM 61:125-135 (2009)) -% No run-time calculation of the correction parametes is currently -% performed in the toolbox. They've been calculated and stored here for a -% series of standard protocols, and are selected at run-time according to -% effective TR and FA values (more precisely: [TR_pdw TR_t1w fa_pdw -% fa_t1w]). If the used TR and FA values don't match any of the predefined -% types, no correction is applied. Additional protocol types and correction -% factors can be computed off-line (see more details below) and added to -% the local defaults file to apply proper RF correction. -% ADVANCED USER ONLY. - -% A) Defining the MPMacq parameters distinguishing the different protocols -%-------------------------------------------------------------------------- -% Using the following parameter order: [TR_pdw TR_t1w fa_pdw fa_t1w] -% NOTE: all tags MUST -% - start with a letter, and -% - include only letters, numbers or underscore, i.e. NO space. -% as these names are used to define a structure fieldname with the protocol -% parameters. -% -% 1) classic FIL protocol (Weiskopf et al., Neuroimage 2011): -% PD-weighted: TR=23.7ms; a=6deg; T1-weighted: TR=18.7ms; a=20deg -hmri_def.MPMacq_set.names{1} = 'Classic FIL protocol'; -hmri_def.MPMacq_set.tags{1} = 'ClassicFIL'; -hmri_def.MPMacq_set.vals{1} = [23.7 18.7 6 20]; -% 2) new FIL/Helms protocol -% PD-weighted: TR=24.5ms; a=5deg; T1-weighted: TR=24.5ms; a=29deg -hmri_def.MPMacq_set.names{2} = 'New FIL/Helms protocol'; -hmri_def.MPMacq_set.tags{2} = 'NewFILHelms'; -hmri_def.MPMacq_set.vals{2} = [24.5 24.5 5 29]; -% 3) Siemens product sequence protocol used in Lausanne (G Krueger) -% PD-weighted: TR=24ms; a=6deg; T1-weighted: TR=19ms; a=20deg -hmri_def.MPMacq_set.names{3} = 'Siemens product Lausanne (GK) protocol'; -hmri_def.MPMacq_set.tags{3} = 'SiemPrLausGK'; -hmri_def.MPMacq_set.vals{3} = [24.0 19.0 6 20]; -% 4) High-res (0.8mm) FIL protocol: -% PD-weighted: TR=23.7ms; a=6deg; T1-weighted: TR=23.7ms; a=28deg -hmri_def.MPMacq_set.names{4} = 'High-res FIL protocol'; -hmri_def.MPMacq_set.tags{4} = 'HResFIL'; -hmri_def.MPMacq_set.vals{4} = [23.7 23.7 6 28]; -% 5)NEW High-res (0.8mm) FIL protocol: -% PD-weighted: TR=25.25ms; a=5deg; T1-weighted: TR=TR=25.25ms; a=29deg -hmri_def.MPMacq_set.names{5} = 'New High-res FIL protocol'; -hmri_def.MPMacq_set.tags{5} = 'NHResFIL'; -hmri_def.MPMacq_set.vals{5} = [25.25 25.25 5 29]; -% 6)NEW 1mm protocol - seq version v2k: -% PD-weighted: TR=24.5ms; a=6deg; T1-weighted: TR=24.5ms; a=21deg -hmri_def.MPMacq_set.names{6} = 'v2k protocol'; -hmri_def.MPMacq_set.tags{6} = 'v2k'; -hmri_def.MPMacq_set.vals{6} = [24.5 24.5 6 21]; -% 7) 800um protocol - seq version v3* released used by MEG group: -% TR = 25ms for all volumes; flipAngles = [6, 21 deg] for PDw and T1w -% Correction parameters below were determined via Bloch-Torrey -% simulations but end result agrees well with EPG-derived correction -% for this RF spoiling increment of 137 degrees. -% See: Callaghan et al. ISMRM, 2015, #1694 -hmri_def.MPMacq_set.names{7} = 'v3star protocol'; -hmri_def.MPMacq_set.tags{7} = 'v3star'; -hmri_def.MPMacq_set.vals{7} = [25 25 6 21]; - -% B) Defining the imperfect spoiling correction parameters for the -% different protocols -%-------------------------------------------------------------------------- -% Antoine Lutti 15/01/09 -% The following correction parameters are based on -% [Preibisch and Deichmann, Magn Reson Med 61:125-135 (2009)]. The values -% for P2_a and P2_b below were obtained using the code supplied by R. -% Deichmann with the experimental parameters used for the standard MPM -% protocol and assuming T2 = 64 ms at 3T. - -% Given the possible confusion and resulting mistake (imperfect spoiling -% correction applied to the wrong sequence), when TR and FA values match -% one of the listed cases below, the option is disabled by default. +% IMPREFECT RF SPOILING CORRECTION +% The RF- and gradient-spoiled gradient echo sequences used to acquire the +% multiparametric mapping (MPM) data apply RF and gradient spoiling to +% destroy unwanted transverse magnetisation. Imperfect spoiling can leave +% residual bias in the apparent R1 map if no further correction is applied +% [Preibisch and Deichmann, MRM 61:125-135, 2009]. +% Correction coefficients are sequence-specific and can be determined by +% simulation [Preibisch and Deichmann 2009] to account for deviation from +% the Ernst equation. +% Since no run-time calculation of the correction parametes is currently +% performed in the toolbox, they've been calculated and stored for a +% series of standard protocols in a separate folder: +% config/local/ISC_parameters +% Additional protocol types and correction factors can be computed on-line +% (see more details below) and added to the local defaults file to apply +% proper RF correction, respectively. +% For the most standard MPM protocols using customised FLASH sequences on Siemens +% scanners, the above mentioned files folder contains files with spoiling correction coefficients +% The correction values P2_a and P2_b were obtained using the code supplied by R. Deichmann with the experimental +% parameters used for the standard MPM protocol and assuming T2 = 64 ms at 3T. % When enabling the imperfect spoiling correction, make sure the -% coefficients retrieved in the list below are definitely calculated for +% coefficients used by choosing the appropriate file are definitely calculated for % the protocol used! +% WARNING: although the TR and FA values are key parameters in simulating the +% imperfect spoiling, correction coefficients cannot be chosen based on TR +% and FA only. Amplitude and duration of gradient spoilers, phase increment +% for RF spoiling and diffusion effects must be taken into account. +% Therefore, two different sequences using identical TR anf FA are unlikely to +% use identical correction coefficients. +% Instead, the coefficients can be calculated seperately for other protocols +% by the additional hMRI module 'Imperfect Spoiling Calc.' that efficiently +% calculates the protocol-specific correction parameters required to account for +% imperfect spoiling. The produced output file can then be set in the local +% defaults file to be incorporated in the processing. +% If not all acquisition parameters are known and vendor sequences at 3 T MRI are +% used, an alternative correction according to Baudrexel et al., 79:3082–3092, +% 2018 can be applied by specifying the RF spoiling phase increment, which +% differs per vendor. Please consult the paper for assumptions and full set +% of correction parameters. +% By default the spoiling correction is disabled. For enabling please set +% up the method of choice in the local defaults file. +% ADVANCED USER ONLY. hmri_def.imperfectSpoilCorr.enabled = false; -% 1) classic FIL protocol (Weiskopf et al., Neuroimage 2011): -hmri_def.imperfectSpoilCorr.ClassicFIL.tag = 'Classic FIL protocol'; -hmri_def.imperfectSpoilCorr.ClassicFIL.P2_a = [78.9228195006542,-101.113338489192,47.8783287525126]; -hmri_def.imperfectSpoilCorr.ClassicFIL.P2_b = [-0.147476233142129,0.126487385091045,0.956824374979504]; -% 2) new FIL/Helms protocol -hmri_def.imperfectSpoilCorr.NewFILHelms.tag = 'New FIL/Helms protocol'; -hmri_def.imperfectSpoilCorr.NewFILHelms.P2_a = [93.455034845930480,-120.5752858196904,55.911077913369060]; -hmri_def.imperfectSpoilCorr.NewFILHelms.P2_b = [-0.167301931434861,0.113507432776106,0.961765216743606]; -% 3) Siemens product sequence protocol used in Lausanne (G Krueger) -hmri_def.imperfectSpoilCorr.SiemPrLausGK.tag = 'Siemens product Lausanne (GK) protocol'; -hmri_def.imperfectSpoilCorr.SiemPrLausGK.P2_a = [67.023102027100880,-86.834117103841540,43.815818592349870]; -hmri_def.imperfectSpoilCorr.SiemPrLausGK.P2_b = [-0.130876849571103,0.117721807209409,0.959180058389875]; -% 4) High-res (0.8mm) FIL protocol: -hmri_def.imperfectSpoilCorr.HResFIL.tag = 'High-res FIL protocol'; -hmri_def.imperfectSpoilCorr.HResFIL.P2_a = [1.317257319014170e+02,-1.699833074433892e+02,73.372595677371650]; -hmri_def.imperfectSpoilCorr.HResFIL.P2_b = [-0.218804328507184,0.178745853134922,0.939514554747592]; -% 5)NEW High-res (0.8mm) FIL protocol: -hmri_def.imperfectSpoilCorr.NHResFIL.tag = 'New High-res FIL protocol'; -hmri_def.imperfectSpoilCorr.NHResFIL.P2_a = [88.8623036106612,-114.526218941363,53.8168602253166]; -hmri_def.imperfectSpoilCorr.NHResFIL.P2_b = [-0.132904017579521,0.113959390779008,0.960799295622202]; -% 6)NEW 1mm protocol - seq version v2k: -hmri_def.imperfectSpoilCorr.v2k.tag = 'v2k protocol'; -hmri_def.imperfectSpoilCorr.v2k.P2_a = [71.2817617982844,-92.2992876164017,45.8278193851731]; -hmri_def.imperfectSpoilCorr.v2k.P2_b = [-0.137859046784839,0.122423212397157,0.957642744668469]; -% 7) 800um protocol - seq version v3* released used by MEG group: -hmri_def.imperfectSpoilCorr.v3star.tag = 'v3star protocol'; -hmri_def.imperfectSpoilCorr.v3star.P2_a = [57.427573706259864,-79.300742898810441,39.218584751863879]; -hmri_def.imperfectSpoilCorr.v3star.P2_b = [-0.121114060111119,0.121684347499374,0.955987357483519]; -% Unknwon protocol -hmri_def.imperfectSpoilCorr.Unknown.tag = 'Unknown protocol. No spoiling correction defined.'; - %-------------------------------------------------------------------------- % B1 mapping processing parameters %-------------------------------------------------------------------------- diff --git a/config/local/ISC_parameters/README.txt b/config/local/ISC_parameters/README.txt new file mode 100644 index 00000000..9de8f8bb --- /dev/null +++ b/config/local/ISC_parameters/README.txt @@ -0,0 +1,14 @@ +IMPREFECT RF SPOILING (ISC) CORRECTION PARAMETERS +(Preibisch and Deichmann, MRM 61:125-135 (2009)) +Correction parameters formerly calculated by Antoine Lutti (15/01/09) +for a series of standard protocols can be selected manually. +Effective TR and FA values (more precisely: [TR_pdw TR_t1w fa_pdw +fa_t1w]) are stored within each ISC file. +Additional protocol types and correction factors can be computed online +and directly chosen as dependency input file to apply proper RF correction. +The stored values for P2_a and P2_b were obtained using the code supplied +by R. Deichmann with the experimental parameters used for the standard +MPM protocol and assuming T2 = 64 ms at 3T. +When enabling the imperfect spoiling correction, make sure the +calculated coefficients are definitely calculated for the protocol used! +ADVANCED USER ONLY. diff --git a/config/local/ISC_parameters/protocol_01_classic_FIL.m b/config/local/ISC_parameters/protocol_01_classic_FIL.m new file mode 100644 index 00000000..e6a0b2cd --- /dev/null +++ b/config/local/ISC_parameters/protocol_01_classic_FIL.m @@ -0,0 +1,9 @@ +% 1) classic FIL protocol (Weiskopf et al., Neuroimage 2011): +% PD-weighted: TR=23.7ms; a=6deg; T1-weighted: TR=18.7ms; a=20deg +global hmri_def +hmri_def.MPMacq_set.name = 'Classic FIL protocol'; +hmri_def.MPMacq_set.tag = 'ClassicFIL'; +hmri_def.MPMacq_set.val = [23.7 18.7 6 20]; +hmri_def.imperfectSpoilCorr.ClassicFIL.tag = 'Classic FIL protocol'; +hmri_def.imperfectSpoilCorr.ClassicFIL.P2_a = [78.9228195006542,-101.113338489192,47.8783287525126]; +hmri_def.imperfectSpoilCorr.ClassicFIL.P2_b = [-0.147476233142129,0.126487385091045,0.956824374979504]; diff --git a/config/local/ISC_parameters/protocol_02_new_FIL_Helms.m b/config/local/ISC_parameters/protocol_02_new_FIL_Helms.m new file mode 100644 index 00000000..4965176d --- /dev/null +++ b/config/local/ISC_parameters/protocol_02_new_FIL_Helms.m @@ -0,0 +1,9 @@ +% 2) new FIL/Helms protocol +% PD-weighted: TR=24.5ms; a=5deg; T1-weighted: TR=24.5ms; a=29deg +global hmri_def +hmri_def.MPMacq_set.name = 'New FIL/Helms protocol'; +hmri_def.MPMacq_set.tag = 'NewFILHelms'; +hmri_def.MPMacq_set.val = [24.5 24.5 5 29]; +hmri_def.imperfectSpoilCorr.NewFILHelms.tag = 'New FIL/Helms protocol'; +hmri_def.imperfectSpoilCorr.NewFILHelms.P2_a = [93.455034845930480,-120.5752858196904,55.911077913369060]; +hmri_def.imperfectSpoilCorr.NewFILHelms.P2_b = [-0.167301931434861,0.113507432776106,0.961765216743606]; diff --git a/config/local/ISC_parameters/protocol_03_Siemens_product_sequence_Lausanne.m b/config/local/ISC_parameters/protocol_03_Siemens_product_sequence_Lausanne.m new file mode 100644 index 00000000..b8d7f1fb --- /dev/null +++ b/config/local/ISC_parameters/protocol_03_Siemens_product_sequence_Lausanne.m @@ -0,0 +1,9 @@ +% 3) Siemens product sequence protocol used in Lausanne (G Krueger) +% PD-weighted: TR=24ms; a=6deg; T1-weighted: TR=19ms; a=20deg +global hmri_def +hmri_def.MPMacq_set.name = 'Siemens product Lausanne (GK) protocol'; +hmri_def.MPMacq_set.tag = 'SiemPrLausGK'; +hmri_def.MPMacq_set.val = [24.0 19.0 6 20]; +hmri_def.imperfectSpoilCorr.SiemPrLausGK.tag = 'Siemens product Lausanne (GK) protocol'; +hmri_def.imperfectSpoilCorr.SiemPrLausGK.P2_a = [67.023102027100880,-86.834117103841540,43.815818592349870]; +hmri_def.imperfectSpoilCorr.SiemPrLausGK.P2_b = [-0.130876849571103,0.117721807209409,0.959180058389875]; diff --git a/config/local/ISC_parameters/protocol_04_high_res_0.8mm_FIL.m b/config/local/ISC_parameters/protocol_04_high_res_0.8mm_FIL.m new file mode 100644 index 00000000..ebb45aef --- /dev/null +++ b/config/local/ISC_parameters/protocol_04_high_res_0.8mm_FIL.m @@ -0,0 +1,9 @@ +% 4) High-res (0.8mm) FIL protocol: +% PD-weighted: TR=23.7ms; a=6deg; T1-weighted: TR=23.7ms; a=28deg +global hmri_def +hmri_def.MPMacq_set.name = 'High-res FIL protocol'; +hmri_def.MPMacq_set.tag = 'HResFIL'; +hmri_def.MPMacq_set.val = [23.7 23.7 6 28]; +hmri_def.imperfectSpoilCorr.HResFIL.tag = 'High-res FIL protocol'; +hmri_def.imperfectSpoilCorr.HResFIL.P2_a = [1.317257319014170e+02,-1.699833074433892e+02,73.372595677371650]; +hmri_def.imperfectSpoilCorr.HResFIL.P2_b = [-0.218804328507184,0.178745853134922,0.939514554747592]; diff --git a/config/local/ISC_parameters/protocol_05_NEW_High_res_0.8mm_FIL.m b/config/local/ISC_parameters/protocol_05_NEW_High_res_0.8mm_FIL.m new file mode 100644 index 00000000..cbf9a936 --- /dev/null +++ b/config/local/ISC_parameters/protocol_05_NEW_High_res_0.8mm_FIL.m @@ -0,0 +1,9 @@ +% 5)NEW High-res (0.8mm) FIL protocol: +% PD-weighted: TR=25.25ms; a=5deg; T1-weighted: TR=TR=25.25ms; a=29deg +global hmri_def +hmri_def.MPMacq_set.name = 'New High-res FIL protocol'; +hmri_def.MPMacq_set.tag = 'NHResFIL'; +hmri_def.MPMacq_set.val = [25.25 25.25 5 29]; +hmri_def.imperfectSpoilCorr.NHResFIL.tag = 'New High-res FIL protocol'; +hmri_def.imperfectSpoilCorr.NHResFIL.P2_a = [88.8623036106612,-114.526218941363,53.8168602253166]; +hmri_def.imperfectSpoilCorr.NHResFIL.P2_b = [-0.132904017579521,0.113959390779008,0.960799295622202]; diff --git a/config/local/ISC_parameters/protocol_06_NEW_1mm_v2k.m b/config/local/ISC_parameters/protocol_06_NEW_1mm_v2k.m new file mode 100644 index 00000000..85629883 --- /dev/null +++ b/config/local/ISC_parameters/protocol_06_NEW_1mm_v2k.m @@ -0,0 +1,10 @@ +% 6)NEW 1mm protocol - seq version v2k: +% PD-weighted: TR=24.5ms; a=6deg; T1-weighted: TR=24.5ms; a=21deg +global hmri_def +hmri_def.MPMacq_set.name = 'v2k protocol'; +hmri_def.MPMacq_set.tag = 'v2k'; +hmri_def.MPMacq_set.val = [24.5 24.5 6 21]; +hmri_def.imperfectSpoilCorr.v2k.tag = 'v2k protocol'; +hmri_def.imperfectSpoilCorr.v2k.P2_a = [71.2817617982844,-92.2992876164017,45.8278193851731]; +hmri_def.imperfectSpoilCorr.v2k.P2_b = [-0.137859046784839,0.122423212397157,0.957642744668469]; + diff --git a/config/local/ISC_parameters/protocol_07_800um_v3star.m b/config/local/ISC_parameters/protocol_07_800um_v3star.m new file mode 100644 index 00000000..73a604b5 --- /dev/null +++ b/config/local/ISC_parameters/protocol_07_800um_v3star.m @@ -0,0 +1,13 @@ +% 7) 800um protocol - seq version v3* released used by MEG group: +% TR = 25ms for all volumes; flipAngles = [6, 21 deg] for PDw and T1w +% Correction parameters below were determined via Bloch-Torrey +% simulations but end result agrees well with EPG-derived correction +% for this RF spoiling increment of 137 degrees. +% See: Callaghan et al. ISMRM, 2015, #1694 +global hmri_def +hmri_def.MPMacq_set.names{7} = 'v3star protocol'; +hmri_def.MPMacq_set.tags{7} = 'v3star'; +hmri_def.MPMacq_set.vals{7} = [25 25 6 21]; +hmri_def.imperfectSpoilCorr.v3star.tag = 'v3star protocol'; +hmri_def.imperfectSpoilCorr.v3star.P2_a = [57.427573706259864,-79.300742898810441,39.218584751863879]; +hmri_def.imperfectSpoilCorr.v3star.P2_b = [-0.121114060111119,0.121684347499374,0.955987357483519]; diff --git a/config/local/hmri_local_defaults.m b/config/local/hmri_local_defaults.m index 89df2151..472f3cfe 100644 --- a/config/local/hmri_local_defaults.m +++ b/config/local/hmri_local_defaults.m @@ -282,123 +282,59 @@ hmri_def.MPMacq.fa_pdw = 6; % [deg] hmri_def.MPMacq.tag = 'v2k'; -% IMPREFECT RF SPOILING CORRECTION PARAMETERS -% (Preibisch and Deichmann, MRM 61:125-135 (2009)) -% No run-time calculation of the correction parametes is currently -% performed in the toolbox. They've been calculated and stored here for a -% series of standard protocols, and are selected at run-time according to -% effective TR and FA values (more precisely: [TR_pdw TR_t1w fa_pdw -% fa_t1w]). If the used TR and FA values don't match any of the predefined -% types, no correction is applied. Additional protocol types and correction -% factors can be computed off-line (see more details below) and added to -% the local defaults file to apply proper RF correction. -% ADVANCED USER ONLY. - -% A) Defining the MPMacq parameters distinguishing the different protocols -%-------------------------------------------------------------------------- -% Using the following parameter order: [TR_pdw TR_t1w fa_pdw fa_t1w] -% NOTE: all tags MUST -% - start with a letter, and -% - include only letters, numbers or underscore, i.e. NO space. -% as these names are used to define a structure fieldname with the protocol -% parameters. -% -% 1) classic FIL protocol (Weiskopf et al., Neuroimage 2011): -% PD-weighted: TR=23.7ms; a=6deg; T1-weighted: TR=18.7ms; a=20deg -hmri_def.MPMacq_set.names{1} = 'Classic FIL protocol'; -hmri_def.MPMacq_set.tags{1} = 'ClassicFIL'; -hmri_def.MPMacq_set.vals{1} = [23.7 18.7 6 20]; -% 2) new FIL/Helms protocol -% PD-weighted: TR=24.5ms; a=5deg; T1-weighted: TR=24.5ms; a=29deg -hmri_def.MPMacq_set.names{2} = 'New FIL/Helms protocol'; -hmri_def.MPMacq_set.tags{2} = 'NewFILHelms'; -hmri_def.MPMacq_set.vals{2} = [24.5 24.5 5 29]; -% 3) Siemens product sequence protocol used in Lausanne (G Krueger) -% PD-weighted: TR=24ms; a=6deg; T1-weighted: TR=19ms; a=20deg -hmri_def.MPMacq_set.names{3} = 'Siemens product Lausanne (GK) protocol'; -hmri_def.MPMacq_set.tags{3} = 'SiemPrLausGK'; -hmri_def.MPMacq_set.vals{3} = [24.0 19.0 6 20]; -% 4) High-res (0.8mm) FIL protocol: -% PD-weighted: TR=23.7ms; a=6deg; T1-weighted: TR=23.7ms; a=28deg -hmri_def.MPMacq_set.names{4} = 'High-res FIL protocol'; -hmri_def.MPMacq_set.tags{4} = 'HResFIL'; -hmri_def.MPMacq_set.vals{4} = [23.7 23.7 6 28]; -% 5)NEW High-res (0.8mm) FIL protocol: -% PD-weighted: TR=25.25ms; a=5deg; T1-weighted: TR=TR=25.25ms; a=29deg -hmri_def.MPMacq_set.names{5} = 'New High-res FIL protocol'; -hmri_def.MPMacq_set.tags{5} = 'NHResFIL'; -hmri_def.MPMacq_set.vals{5} = [25.25 25.25 5 29]; -% 6)NEW 1mm protocol - seq version v2k: -% PD-weighted: TR=24.5ms; a=6deg; T1-weighted: TR=24.5ms; a=21deg -hmri_def.MPMacq_set.names{6} = 'v2k protocol'; -hmri_def.MPMacq_set.tags{6} = 'v2k'; -hmri_def.MPMacq_set.vals{6} = [24.5 24.5 6 21]; -% 7) 800um protocol - seq version v3* released used by MEG group: -% TR = 25ms for all volumes; flipAngles = [6, 21 deg] for PDw and T1w -% Correction parameters below were determined via Bloch-Torrey -% simulations but end result agrees well with EPG-derived correction -% for this RF spoiling increment of 137 degrees. -% See: Callaghan et al. ISMRM, 2015, #1694 -hmri_def.MPMacq_set.names{7} = 'v3star protocol'; -hmri_def.MPMacq_set.tags{7} = 'v3star'; -hmri_def.MPMacq_set.vals{7} = [25 25 6 21]; - -% B) Defining the imperfect spoiling correction parameters for the -% different protocols -%-------------------------------------------------------------------------- -% Antoine Lutti 15/01/09 -% The following correction parameters are based on -% [Preibisch and Deichmann, Magn Reson Med 61:125-135 (2009)]. The values -% for P2_a and P2_b below were obtained using the code supplied by R. -% Deichmann with the experimental parameters used for the standard MPM -% protocol and assuming T2 = 64 ms at 3T. - -% Given the possible confusion and resulting mistake (imperfect spoiling -% correction applied to the wrong sequence), when TR and FA values match -% one of the listed cases below, the option is disabled by default. +% IMPREFECT RF SPOILING CORRECTION +% The RF- and gradient-spoiled gradient echo sequences used to acquire the +% multiparametric mapping (MPM) data apply RF and gradient spoiling to +% destroy unwanted transverse magnetisation. Imperfect spoiling can leave +% residual bias in the apparent R1 map if no further correction is applied +% [Preibisch and Deichmann, MRM 61:125-135, 2009]. +% Correction coefficients are sequence-specific and can be determined by +% simulation [Preibisch and Deichmann 2009] to account for deviation from +% the Ernst equation. +% Since no run-time calculation of the correction parametes is currently +% performed in the toolbox, they've been calculated and stored for a +% series of standard protocols in a separate folder: +% config/local/ISC_parameters +% Additional protocol types and correction factors can be computed on-line +% (see more details below) and added to the local defaults file to apply +% proper RF correction, respectively. +% For the most standard MPM protocols using customised FLASH sequences on Siemens +% scanners, the above mentioned files folder contains files with spoiling correction coefficients +% The correction values P2_a and P2_b were obtained using the code supplied by R. Deichmann with the experimental +% parameters used for the standard MPM protocol and assuming T2 = 64 ms at 3T. % When enabling the imperfect spoiling correction, make sure the -% coefficients retrieved in the list below are definitely calculated for +% coefficients used by choosing the appropriate file are definitely calculated for % the protocol used! -hmri_def.imperfectSpoilCorr.enabled = false; - -% 1) classic FIL protocol (Weiskopf et al., Neuroimage 2011): -hmri_def.imperfectSpoilCorr.ClassicFIL.tag = 'Classic FIL protocol'; -hmri_def.imperfectSpoilCorr.ClassicFIL.P2_a = [78.9228195006542,-101.113338489192,47.8783287525126]; -hmri_def.imperfectSpoilCorr.ClassicFIL.P2_b = [-0.147476233142129,0.126487385091045,0.956824374979504]; -hmri_def.imperfectSpoilCorr.ClassicFIL.enabled = hmri_def.imperfectSpoilCorr.enabled; -% 2) new FIL/Helms protocol -hmri_def.imperfectSpoilCorr.NewFILHelms.tag = 'New FIL/Helms protocol'; -hmri_def.imperfectSpoilCorr.NewFILHelms.P2_a = [93.455034845930480,-120.5752858196904,55.911077913369060]; -hmri_def.imperfectSpoilCorr.NewFILHelms.P2_b = [-0.167301931434861,0.113507432776106,0.961765216743606]; -hmri_def.imperfectSpoilCorr.NewFILHelms.enabled = hmri_def.imperfectSpoilCorr.enabled; -% 3) Siemens product sequence protocol used in Lausanne (G Krueger) -hmri_def.imperfectSpoilCorr.SiemPrLausGK.tag = 'Siemens product Lausanne (GK) protocol'; -hmri_def.imperfectSpoilCorr.SiemPrLausGK.P2_a = [67.023102027100880,-86.834117103841540,43.815818592349870]; -hmri_def.imperfectSpoilCorr.SiemPrLausGK.P2_b = [-0.130876849571103,0.117721807209409,0.959180058389875]; -hmri_def.imperfectSpoilCorr.SiemPrLausGK.enabled = hmri_def.imperfectSpoilCorr.enabled; -% 4) High-res (0.8mm) FIL protocol: -hmri_def.imperfectSpoilCorr.HResFIL.tag = 'High-res FIL protocol'; -hmri_def.imperfectSpoilCorr.HResFIL.P2_a = [1.317257319014170e+02,-1.699833074433892e+02,73.372595677371650]; -hmri_def.imperfectSpoilCorr.HResFIL.P2_b = [-0.218804328507184,0.178745853134922,0.939514554747592]; -hmri_def.imperfectSpoilCorr.HResFIL.enabled = hmri_def.imperfectSpoilCorr.enabled; -% 5)NEW High-res (0.8mm) FIL protocol: -hmri_def.imperfectSpoilCorr.NHResFIL.tag = 'New High-res FIL protocol'; -hmri_def.imperfectSpoilCorr.NHResFIL.P2_a = [88.8623036106612,-114.526218941363,53.8168602253166]; -hmri_def.imperfectSpoilCorr.NHResFIL.P2_b = [-0.132904017579521,0.113959390779008,0.960799295622202]; -hmri_def.imperfectSpoilCorr.NHResFIL.enabled = hmri_def.imperfectSpoilCorr.enabled; -% 6)NEW 1mm protocol - seq version v2k: -hmri_def.imperfectSpoilCorr.v2k.tag = 'v2k protocol'; -hmri_def.imperfectSpoilCorr.v2k.P2_a = [71.2817617982844,-92.2992876164017,45.8278193851731]; -hmri_def.imperfectSpoilCorr.v2k.P2_b = [-0.137859046784839,0.122423212397157,0.957642744668469]; -hmri_def.imperfectSpoilCorr.v2k.enabled = hmri_def.imperfectSpoilCorr.enabled; -% 7) 800um protocol - seq version v3* released used by MEG group: -hmri_def.imperfectSpoilCorr.v3star.tag = 'v3star protocol'; -hmri_def.imperfectSpoilCorr.v3star.P2_a = [57.427573706259864,-79.300742898810441,39.218584751863879]; -hmri_def.imperfectSpoilCorr.v3star.P2_b = [-0.121114060111119,0.121684347499374,0.955987357483519]; -hmri_def.imperfectSpoilCorr.v3star.enabled = hmri_def.imperfectSpoilCorr.enabled; -% Unknwon protocol -hmri_def.imperfectSpoilCorr.Unknown.tag = 'Unknown protocol. No spoiling correction defined.'; -hmri_def.imperfectSpoilCorr.Unknown.enabled = false; +% WARNING: although the TR and FA values are key parameters in simulating the +% imperfect spoiling, correction coefficients cannot be chosen based on TR +% and FA only. Amplitude and duration of gradient spoilers, phase increment +% for RF spoiling and diffusion effects must be taken into account. +% Therefore, two different sequences using identical TR anf FA are unlikely to +% use identical correction coefficients. + +% Example setting for choosing the classic FIL protocol (uncomment following two lines to use): +% hmri_def.imperfectSpoilCorr.enabled = true; +% hmri_def.imperfectSpoilCorr.iscfile = fullfile(fileparts(fileparts(mfilename('fullpath'))),'ISC_parameters','protocol_01_classic_FIL.m'); + +% Instead, the coefficients can be calculated seperately for other protocols +% by the additional hMRI module 'Imperfect Spoiling Calc.' that efficiently +% calculates the protocol-specific correction parameters required to account for +% imperfect spoiling. The produced output file can then be set in the local +% defaults file to be incorporated in the processing. +% If not all acquisition parameters are known and vendor sequences at 3 T MRI are +% used, an alternative correction according to Baudrexel et al., 79:3082–3092, +% 2018 can be applied by specifying the RF spoiling phase increment, which +% differs per vendor. Please consult the paper for assumptions and full set +% of correction parameters. Following parameter values (in degree) are possible: +% 50 (Siemens), 117 (General Electrics), and 150 (Philips). + +% Example setting for choosing the phase increment at Philips scanners (uncomment following two lines to use): +% hmri_def.imperfectSpoilCorr.enabled = true; +% hmri_def.imperfectSpoilCorr.iscphase = 150; + +% ADVANCED USER ONLY. + + %========================================================================== % Maps processing parameters @@ -426,4 +362,4 @@ % Number of Gaussians per tissue class hmri_def.proc.nGauss = [2 2 2 3 4 2]; % originally in SPM [1 1 2 3 4 2] -end \ No newline at end of file +end diff --git a/config/local/hmri_local_defaults_iscfile.m b/config/local/hmri_local_defaults_iscfile.m new file mode 100644 index 00000000..d2e9c1c1 --- /dev/null +++ b/config/local/hmri_local_defaults_iscfile.m @@ -0,0 +1,366 @@ +function hmri_local_defaults +% PURPOSE +% To set user-defined (site- or protocol-specific) defaults parameters +% which are used by the hMRI toolbox. Customized processing parameters can +% be defined, overwriting defaults from hmri_defaults. Acquisition +% protocols can be specified here as a fallback solution when no metadata +% are available. Note that the use of metadata is strongly recommended. +% +% RECOMMENDATIONS +% Parameters defined in this file are identical, initially, to the ones +% defined in hmri_defaults.m. It is recommended, when modifying this file, +% to remove all unchanged entries and save the file with a meaningful name. +% This will help you identifying the appropriate defaults to be used for +% each protocol, and will improve the readability of the file by pointing +% to the modified parameters only. +% +% WARNING +% Modification of the defaults parameters may impair the integrity of the +% toolbox, leading to unexpected behaviour. ONLY RECOMMENDED FOR ADVANCED +% USERS - i.e. who have a good knowledge of the underlying algorithms and +% implementation. The SAME SET OF DEFAULT PARAMETERS must be used to +% process uniformly all the data from a given study. +% +% HOW DOES IT WORK? +% The modified defaults file can be selected using the "Configure toolbox" +% branch of the hMRI-Toolbox. For customization of B1 processing +% parameters, type "help hmri_b1_standard_defaults.m". +% +% DOCUMENTATION +% A brief description of each parameter is provided together with +% guidelines and recommendations to modify these parameters. With few +% exceptions, parameters should ONLY be MODIFIED and customized BY ADVANCED +% USERS, having a good knowledge of the underlying algorithms and +% implementation. +% Please refer to the documentation in the github WIKI for more details. +%__________________________________________________________________________ +% Written by E. Balteau, 2017. +% Cyclotron Research Centre, University of Liege, Belgium + +% Global hmri_def variable used across the whole toolbox +global hmri_def + +% Specify the research centre & scanner. Not mandatory. +hmri_def.centre = 'centre' ; % 'fil', 'lren', 'crc', 'sciz', 'cbs', ... +hmri_def.scanner = 'scanner name' ; % e.g. 'prisma', 'allegra', 'terra', 'achieva', ... + +%========================================================================== +% Common processing parameters +%========================================================================== + +% cleanup temporary directories. If set to true, all temporary directories +% are deleted at the end of map creation, only the "Results" directory and +% "Supplementary" subdirectory are kept. Setting "cleanup" to "false" might +% be convenient if one desires to have a closer look at intermediate +% processing steps. Otherwise "cleanup = true" is recommended for saving +% disk space. +hmri_def.cleanup = true; +% settings for JSON metadata: by default, separate JSON files are used to +% store the metadata (information on data acquisition and processing, +% tracking of input and output files), as JSON-formatted, tab-indented +% text. The following settings are recommended. No modification currently +% foreseen as useful... +hmri_def.json = struct('extended',false,'separate',true,'anonym','none',... + 'overwrite',true, 'indent','\t'); +% recommended TPM for segmentation and spatial processing. The hMRI toolbox +% provides a series of tissue probability maps. These TPMs could be +% replaced by other TPMs, to better match the population studied. +% ADVANCED USER ONLY. +hmri_def.TPM = fullfile(fileparts(fileparts(fileparts(mfilename('fullpath')))),'etpm','eTPM.nii'); +% default template for auto-reorientation. The template can be selected +% within the Auto-reorient module. The following is the default suggested +% for T1w images. Please refer to the Auto-reorient documentation for an +% appropriate choice of the template. +hmri_def.autoreorient_template = {fullfile(spm('dir'),'canonical','avg152T1.nii')}; + +%========================================================================== +% Default parameters for segmentation +% ADVANCED USERS ONLY! +% hmri_def.segment is effectively the job to be handed to spm_preproc_run +% By default, parameters are set to +% - create tissue class images (c*) in the native space of the source image +% (tissue(*).native = [1 0]) for tissue classes 1-5 +% - save both BiasField and BiasCorrected volume (channel.write = [1 1]) +% - recommended values from SPM12 (October 2017) +%========================================================================== + +% hmri_def.segment.channel.vols = cell array of file names, +% must be defined before calling spm_preproc_run +hmri_def.segment.channel.biasreg = 0.001; +hmri_def.segment.channel.biasfwhm = 60; +% hmri_def.segment.channel.write = [0 0]; % save nothing +% hmri_def.segment.channel.write = [1 0]; % save BiasField +% hmri_def.segment.channel.write = [0 1]; % save BiasCorrected volume +hmri_def.segment.channel.write = [1 1]; % save BiasField and BiasCorrected volume + +hmri_def.segment.tissue(1).tpm = {[hmri_def.TPM ',1']}; +hmri_def.segment.tissue(1).ngaus = 1; +hmri_def.segment.tissue(1).native = [1 0]; +hmri_def.segment.tissue(1).warped = [0 0]; +hmri_def.segment.tissue(2).tpm = {[hmri_def.TPM ',2']}; +hmri_def.segment.tissue(2).ngaus = 1; +hmri_def.segment.tissue(2).native = [1 0]; +hmri_def.segment.tissue(2).warped = [0 0]; +hmri_def.segment.tissue(3).tpm = {[hmri_def.TPM ',3']}; +hmri_def.segment.tissue(3).ngaus = 2; +hmri_def.segment.tissue(3).native = [1 0]; +hmri_def.segment.tissue(3).warped = [0 0]; +hmri_def.segment.tissue(4).tpm = {[hmri_def.TPM ',4']}; +hmri_def.segment.tissue(4).ngaus = 3; +hmri_def.segment.tissue(4).native = [1 0]; +hmri_def.segment.tissue(4).warped = [0 0]; +hmri_def.segment.tissue(5).tpm = {[hmri_def.TPM ',5']}; +hmri_def.segment.tissue(5).ngaus = 4; +hmri_def.segment.tissue(5).native = [1 0]; +hmri_def.segment.tissue(5).warped = [0 0]; +hmri_def.segment.tissue(6).tpm = {[hmri_def.TPM ',6']}; +hmri_def.segment.tissue(6).ngaus = 2; +hmri_def.segment.tissue(6).native = [0 0]; +hmri_def.segment.tissue(6).warped = [0 0]; +hmri_def.segment.warp.mrf = 1; +hmri_def.segment.warp.cleanup = 1; +hmri_def.segment.warp.reg = [0 0.001 0.5 0.05 0.2]; +hmri_def.segment.warp.affreg = 'mni'; +hmri_def.segment.warp.fwhm = 0; +hmri_def.segment.warp.samp = 3; +hmri_def.segment.warp.write = [0 0]; + +%========================================================================== +% R1/PD/R2s/MT map creation parameters +%========================================================================== + +%-------------------------------------------------------------------------- +% Coregistration of all input images to the average (or TE=0 fit) PDw image +%-------------------------------------------------------------------------- +% The coregistration step can be disabled using the following flag (not +% recommended). ADVANCED USER ONLY. +hmri_def.coreg2PDw = 1; + +%-------------------------------------------------------------------------- +% Coregistration parameters for weighted images to the average +% (or TE=0 fit) PDw image; see spm_coreg.m for details +%-------------------------------------------------------------------------- +% For high resolution (< 0.8 mm) data, we have found that adding extra +% steps to the registration can improve map quality, e.g. for 0.5 mm data +% we have found +% hmri_def.coreg_flags.sep=[4 2 1 0.5]; +% to be needed to avoid registration artefacts, but this comes at the +% expense of slowing down map creation +hmri_def.coreg_flags.sep=[4 2]; +hmri_def.coreg_flags.cost_fun = 'nmi'; +hmri_def.coreg_flags.fwhm = [7 7]; + +%-------------------------------------------------------------------------- +% Coregistration parameters for B1 bias maps to the average (or TE=0 fit) +% PDw image; see spm_coreg.m for details +%-------------------------------------------------------------------------- +hmri_def.coreg_bias_flags.sep=[4 2]; +hmri_def.coreg_bias_flags.cost_fun = 'nmi'; +hmri_def.coreg_bias_flags.fwhm = [7 7]; + +%-------------------------------------------------------------------------- +% Ordinary Least Squares & fit at TE=0 +%-------------------------------------------------------------------------- +% create an Ordinary Least Squares R2* map. The ESTATICS model is applied +% to calculate R2* map from all available contrasts. +% ADVANCED USER ONLY. +hmri_def.R2sOLS = 1; + +% Minimum number of echoes to calculate R2s map. Strictly speaking, the +% minimum is 2. For a robust estimation, the minimum number of echoes +% required depends on many factors, amongst which: +% - SNR/resolution +% - distribution/spacing between TEs: note that early echoes are more +% affected by the specific contrast, violating the assumption of a common +% decay between contrasts. +% - number of contrasts available (fewer echoes per contrast required for 3 +% (PDw, T1w, MTw) contrasts as compared to 2 or even 1) +% To be on the safe side, a minimum of 6 echoes is recommended (ESTATICS +% paper). Further studies are required to come up with more detailed and +% informed guidelines. Use fewer echoes at your own risk...! +hmri_def.neco4R2sfit = 4; + +% Define a coherent interpolation factor used all through the map creation +% process. Default is 3, but if you want to keep SNR and resolution as far +% as possible the same, it is recommended to use sinc interpolation (at +% least -4, in Siawoosh's experience -7 gives decent results). +% ADVANCED USER ONLY. +hmri_def.interp = 3; + +% Define the OLS fit as default. OLS fit at TE=0 for each contrast is used +% instead of averaged contrast images for the map calculation. +% ADVANCED USER ONLY. +hmri_def.fullOLS = true; + +%-------------------------------------------------------------------------- +% Usage of UNICORT-derived B1 maps for PD and/or MT maps calculation +% ADVANCED USER ONLY. +% WARNING: this method has not been validated for PD and MT calculation! +%-------------------------------------------------------------------------- +hmri_def.UNICORT.PD = false; +hmri_def.UNICORT.MT = false; + +%-------------------------------------------------------------------------- +% PD maps processing parameters +% ADVANCED USER ONLY. +%-------------------------------------------------------------------------- +hmri_def.PDproc.calibr = 1; % Calibration of the PD map (if PDw, T1w, + % B1map available and RF sensitivity bias correction applied somehow) + % based on PD(WM) = 69% [Tofts 2003]. +hmri_def.PDproc.WBMaskTh = 0.1; % Threshold for calculation of whole-brain mask from TPMs +hmri_def.PDproc.WMMaskTh = 0.95; % Threshold for calculation of white-matter mask from TPMs +hmri_def.PDproc.biasreg = 10^(-5); +hmri_def.PDproc.biasfwhm = 50; +hmri_def.PDproc.nr_echoes_forA = 6; % NOTE: in order to minimize R2* bias + % on the PD estimates and gain in robustness for bias-field + % correction, the number of echoes should be minimum ("average" + % calculated over the first echo only) for PD calculation. However, + % with T2*-weighting bias correction (see below), a higher number of + % echoes is preferred in order to provide good SNR. Note that when + % "fullOLS = true", this parameter has no impact whatsovever. +hmri_def.PDproc.T2scorr = 1; % to correct A map for T2*-weighting bias + % before PD map calculation. The correction is not required when + % "fullOLS = true" (TE=0 fit) and will be automatically disabled. + +%-------------------------------------------------------------------------- +% RF sensitivity bias correction: kernel for sensitivity map smoothing. +% ADVANCED USER ONLY. +%-------------------------------------------------------------------------- +hmri_def.RFsens.smooth_kernel = 12; + +%-------------------------------------------------------------------------- +% quantitative maps: quality evaluation and realignment to MNI +%-------------------------------------------------------------------------- +% creates a matlab structure containing markers of data quality +hmri_def.qMRI_maps.QA = false; +% realigns qMRI maps to MNI: the following parameter corresponds to the +% realignment implemented as part of the map calculation (see +% hmri_create_MTProt.m). Left here for backward compatibility. It is +% STRONGLY RECOMMENDED to reorient all images prior any processing using +% the Auto-Reorient module provided with the toolbox (type "help +% hmri_autoreorient" for details or open the SPM > Tools > hMRI Tools > +% Auto-Reorient module in the Batch GUI). ADVANCED USER ONLY. +hmri_def.qMRI_maps.ACPCrealign = 0; + +%-------------------------------------------------------------------------- +% Threshold values for qMRI maps +% The thresholds are meant to discard outliers generally due to low SNR in +% some brain areas, leading to physical non-sense values. Thresholding is +% required to process further the maps generated, when e.g. used +% segmentation algorithms make assumptions incompatible with existing +% outliers. +% NOTE that thresholding modifies the signal distribution and may alter +% the statistical results. +% ADVANCED USER ONLY. +%-------------------------------------------------------------------------- +hmri_def.qMRI_maps_thresh.R1 = 2000; % 1000*[s-1] +hmri_def.qMRI_maps_thresh.A = 10^5; % [a.u.] based on input images with intensities ranging approx. [0 4096]. +hmri_def.qMRI_maps_thresh.R2s = 10; % 1000*[s-1] +hmri_def.qMRI_maps_thresh.MTR = 50; +hmri_def.qMRI_maps_thresh.MTR_synt = 50; +hmri_def.qMRI_maps_thresh.MT = 5; % [p.u.] + +%-------------------------------------------------------------------------- +% MPM acquisition parameters and RF spoiling correction parameters +%-------------------------------------------------------------------------- +% ACQUISITION PARAMETERS: these values are initialised with defaults (v2k +% protocol - Prisma) and are updated at run-time with actual acquisition +% values (see hmri_create_MTProt.m). If TR/TE/FA cannot be determined from +% the input images, the following values will be used. If they don't match +% your own protocol values and if no TR/TE/FA values can be retrieved by +% the toolbox from your data, the following values should be adapted in the +% local defaults file. +% ADVANCED USER ONLY +hmri_def.MPMacq.TE_mtw = [2.34:2.34:14.04]'; %#ok<*NBRAK> % [ms] defined as column vector! +hmri_def.MPMacq.TE_t1w = [2.34:2.34:18.72]'; % [ms] +hmri_def.MPMacq.TE_pdw = [2.34:2.34:18.72]'; % [ms] +hmri_def.MPMacq.TR_mtw = 24.5; % [ms] +hmri_def.MPMacq.TR_t1w = 24.5; % [ms] +hmri_def.MPMacq.TR_pdw = 24.5; % [ms] +hmri_def.MPMacq.fa_mtw = 6; % [deg] +hmri_def.MPMacq.fa_t1w = 21; % [deg] +hmri_def.MPMacq.fa_pdw = 6; % [deg] +hmri_def.MPMacq.tag = 'v2k'; + +% IMPREFECT RF SPOILING CORRECTION +% The RF- and gradient-spoiled gradient echo sequences used to acquire the +% multiparametric mapping (MPM) data apply RF and gradient spoiling to +% destroy unwanted transverse magnetisation. Imperfect spoiling can leave +% residual bias in the apparent R1 map if no further correction is applied +% [Preibisch and Deichmann, MRM 61:125-135, 2009]. +% Correction coefficients are sequence-specific and can be determined by +% simulation [Preibisch and Deichmann 2009] to account for deviation from +% the Ernst equation. +% Since no run-time calculation of the correction parametes is currently +% performed in the toolbox, they've been calculated and stored for a +% series of standard protocols in a separate folder: +% config/local/ISC_parameters +% Additional protocol types and correction factors can be computed on-line +% (see more details below) and added to the local defaults file to apply +% proper RF correction, respectively. +% For the most standard MPM protocols using customised FLASH sequences on Siemens +% scanners, the above mentioned files folder contains files with spoiling correction coefficients +% The correction values P2_a and P2_b were obtained using the code supplied by R. Deichmann with the experimental +% parameters used for the standard MPM protocol and assuming T2 = 64 ms at 3T. +% When enabling the imperfect spoiling correction, make sure the +% coefficients used by choosing the appropriate file are definitely calculated for +% the protocol used! +% WARNING: although the TR and FA values are key parameters in simulating the +% imperfect spoiling, correction coefficients cannot be chosen based on TR +% and FA only. Amplitude and duration of gradient spoilers, phase increment +% for RF spoiling and diffusion effects must be taken into account. +% Therefore, two different sequences using identical TR anf FA are unlikely to +% use identical correction coefficients. + +% Example setting for choosing the classic FIL protocol (uncomment for usage): +hmri_def.imperfectSpoilCorr.enabled = true; +hmri_def.imperfectSpoilCorr.iscfile = fullfile(fileparts(mfilename('fullpath')),'ISC_parameters','protocol_01_classic_FIL.m'); + +% Instead, the coefficients can be calculated seperately for other protocols +% by the additional hMRI module 'Imperfect Spoiling Calc.' that efficiently +% calculates the protocol-specific correction parameters required to account for +% imperfect spoiling. The produced output file can then be set in the local +% defaults file to be incorporated in the processing. +% If not all acquisition parameters are known and vendor sequences at 3 T MRI are +% used, an alternative correction according to Baudrexel et al., 79:3082–3092, +% 2018 can be applied by specifying the RF spoiling phase increment, which +% differs per vendor. Please consult the paper for assumptions and full set +% of correction parameters. Following parameter values (in degree) are possible: +% 50 (Siemens), 117 (General Electrics), and 150 (Philips). + + +% Example setting for choosing the phase increment at Philips scanners: +% hmri_def.imperfectSpoilCorr.enabled = true; +% hmri_def.imperfectSpoilCorr.iscphase = 150; + +% ADVANCED USER ONLY. + + + +%========================================================================== +% Maps processing parameters +%========================================================================== + +%-------------------------------------------------------------------------- +% US segmentation parameters +%-------------------------------------------------------------------------- + +% recommended TPM for segmentation +hmri_def.proc.TPM = hmri_def.TPM ; +% Use the same as for the maps creation but one could (want to) use another +% one at some point. +% Map creation works with "standard" weighted-MR images to build the +% parametric maps. In the end these parametric maps taken together for a +% multichannel-segmention could show more details (for example subcortical +% nuclei?) and would therefore require a specific TPM. This TPM is of +% course still to be built at the moment... + +% Flags to write out posterior tissue classes in native & warped space +% - GM/WM/CSF -> write warped, mod+unmod, and native, native+dartelImp. +% - others -> nothing +hmri_def.proc.w_native = [[1 1];[1 1];[1 1];[0 0];[0 0];[0 0]]; +hmri_def.proc.w_warped = [[1 1];[1 1];[1 1];[0 0];[0 0];[0 0]]; +% Number of Gaussians per tissue class +hmri_def.proc.nGauss = [2 2 2 3 4 2]; % originally in SPM [1 1 2 3 4 2] + +end diff --git a/config/local/hmri_local_defaults_iscphase.m b/config/local/hmri_local_defaults_iscphase.m new file mode 100644 index 00000000..27c7bf44 --- /dev/null +++ b/config/local/hmri_local_defaults_iscphase.m @@ -0,0 +1,366 @@ +function hmri_local_defaults +% PURPOSE +% To set user-defined (site- or protocol-specific) defaults parameters +% which are used by the hMRI toolbox. Customized processing parameters can +% be defined, overwriting defaults from hmri_defaults. Acquisition +% protocols can be specified here as a fallback solution when no metadata +% are available. Note that the use of metadata is strongly recommended. +% +% RECOMMENDATIONS +% Parameters defined in this file are identical, initially, to the ones +% defined in hmri_defaults.m. It is recommended, when modifying this file, +% to remove all unchanged entries and save the file with a meaningful name. +% This will help you identifying the appropriate defaults to be used for +% each protocol, and will improve the readability of the file by pointing +% to the modified parameters only. +% +% WARNING +% Modification of the defaults parameters may impair the integrity of the +% toolbox, leading to unexpected behaviour. ONLY RECOMMENDED FOR ADVANCED +% USERS - i.e. who have a good knowledge of the underlying algorithms and +% implementation. The SAME SET OF DEFAULT PARAMETERS must be used to +% process uniformly all the data from a given study. +% +% HOW DOES IT WORK? +% The modified defaults file can be selected using the "Configure toolbox" +% branch of the hMRI-Toolbox. For customization of B1 processing +% parameters, type "help hmri_b1_standard_defaults.m". +% +% DOCUMENTATION +% A brief description of each parameter is provided together with +% guidelines and recommendations to modify these parameters. With few +% exceptions, parameters should ONLY be MODIFIED and customized BY ADVANCED +% USERS, having a good knowledge of the underlying algorithms and +% implementation. +% Please refer to the documentation in the github WIKI for more details. +%__________________________________________________________________________ +% Written by E. Balteau, 2017. +% Cyclotron Research Centre, University of Liege, Belgium + +% Global hmri_def variable used across the whole toolbox +global hmri_def + +% Specify the research centre & scanner. Not mandatory. +hmri_def.centre = 'centre' ; % 'fil', 'lren', 'crc', 'sciz', 'cbs', ... +hmri_def.scanner = 'scanner name' ; % e.g. 'prisma', 'allegra', 'terra', 'achieva', ... + +%========================================================================== +% Common processing parameters +%========================================================================== + +% cleanup temporary directories. If set to true, all temporary directories +% are deleted at the end of map creation, only the "Results" directory and +% "Supplementary" subdirectory are kept. Setting "cleanup" to "false" might +% be convenient if one desires to have a closer look at intermediate +% processing steps. Otherwise "cleanup = true" is recommended for saving +% disk space. +hmri_def.cleanup = true; +% settings for JSON metadata: by default, separate JSON files are used to +% store the metadata (information on data acquisition and processing, +% tracking of input and output files), as JSON-formatted, tab-indented +% text. The following settings are recommended. No modification currently +% foreseen as useful... +hmri_def.json = struct('extended',false,'separate',true,'anonym','none',... + 'overwrite',true, 'indent','\t'); +% recommended TPM for segmentation and spatial processing. The hMRI toolbox +% provides a series of tissue probability maps. These TPMs could be +% replaced by other TPMs, to better match the population studied. +% ADVANCED USER ONLY. +hmri_def.TPM = fullfile(fileparts(fileparts(fileparts(mfilename('fullpath')))),'etpm','eTPM.nii'); +% default template for auto-reorientation. The template can be selected +% within the Auto-reorient module. The following is the default suggested +% for T1w images. Please refer to the Auto-reorient documentation for an +% appropriate choice of the template. +hmri_def.autoreorient_template = {fullfile(spm('dir'),'canonical','avg152T1.nii')}; + +%========================================================================== +% Default parameters for segmentation +% ADVANCED USERS ONLY! +% hmri_def.segment is effectively the job to be handed to spm_preproc_run +% By default, parameters are set to +% - create tissue class images (c*) in the native space of the source image +% (tissue(*).native = [1 0]) for tissue classes 1-5 +% - save both BiasField and BiasCorrected volume (channel.write = [1 1]) +% - recommended values from SPM12 (October 2017) +%========================================================================== + +% hmri_def.segment.channel.vols = cell array of file names, +% must be defined before calling spm_preproc_run +hmri_def.segment.channel.biasreg = 0.001; +hmri_def.segment.channel.biasfwhm = 60; +% hmri_def.segment.channel.write = [0 0]; % save nothing +% hmri_def.segment.channel.write = [1 0]; % save BiasField +% hmri_def.segment.channel.write = [0 1]; % save BiasCorrected volume +hmri_def.segment.channel.write = [1 1]; % save BiasField and BiasCorrected volume + +hmri_def.segment.tissue(1).tpm = {[hmri_def.TPM ',1']}; +hmri_def.segment.tissue(1).ngaus = 1; +hmri_def.segment.tissue(1).native = [1 0]; +hmri_def.segment.tissue(1).warped = [0 0]; +hmri_def.segment.tissue(2).tpm = {[hmri_def.TPM ',2']}; +hmri_def.segment.tissue(2).ngaus = 1; +hmri_def.segment.tissue(2).native = [1 0]; +hmri_def.segment.tissue(2).warped = [0 0]; +hmri_def.segment.tissue(3).tpm = {[hmri_def.TPM ',3']}; +hmri_def.segment.tissue(3).ngaus = 2; +hmri_def.segment.tissue(3).native = [1 0]; +hmri_def.segment.tissue(3).warped = [0 0]; +hmri_def.segment.tissue(4).tpm = {[hmri_def.TPM ',4']}; +hmri_def.segment.tissue(4).ngaus = 3; +hmri_def.segment.tissue(4).native = [1 0]; +hmri_def.segment.tissue(4).warped = [0 0]; +hmri_def.segment.tissue(5).tpm = {[hmri_def.TPM ',5']}; +hmri_def.segment.tissue(5).ngaus = 4; +hmri_def.segment.tissue(5).native = [1 0]; +hmri_def.segment.tissue(5).warped = [0 0]; +hmri_def.segment.tissue(6).tpm = {[hmri_def.TPM ',6']}; +hmri_def.segment.tissue(6).ngaus = 2; +hmri_def.segment.tissue(6).native = [0 0]; +hmri_def.segment.tissue(6).warped = [0 0]; +hmri_def.segment.warp.mrf = 1; +hmri_def.segment.warp.cleanup = 1; +hmri_def.segment.warp.reg = [0 0.001 0.5 0.05 0.2]; +hmri_def.segment.warp.affreg = 'mni'; +hmri_def.segment.warp.fwhm = 0; +hmri_def.segment.warp.samp = 3; +hmri_def.segment.warp.write = [0 0]; + +%========================================================================== +% R1/PD/R2s/MT map creation parameters +%========================================================================== + +%-------------------------------------------------------------------------- +% Coregistration of all input images to the average (or TE=0 fit) PDw image +%-------------------------------------------------------------------------- +% The coregistration step can be disabled using the following flag (not +% recommended). ADVANCED USER ONLY. +hmri_def.coreg2PDw = 1; + +%-------------------------------------------------------------------------- +% Coregistration parameters for weighted images to the average +% (or TE=0 fit) PDw image; see spm_coreg.m for details +%-------------------------------------------------------------------------- +% For high resolution (< 0.8 mm) data, we have found that adding extra +% steps to the registration can improve map quality, e.g. for 0.5 mm data +% we have found +% hmri_def.coreg_flags.sep=[4 2 1 0.5]; +% to be needed to avoid registration artefacts, but this comes at the +% expense of slowing down map creation +hmri_def.coreg_flags.sep=[4 2]; +hmri_def.coreg_flags.cost_fun = 'nmi'; +hmri_def.coreg_flags.fwhm = [7 7]; + +%-------------------------------------------------------------------------- +% Coregistration parameters for B1 bias maps to the average (or TE=0 fit) +% PDw image; see spm_coreg.m for details +%-------------------------------------------------------------------------- +hmri_def.coreg_bias_flags.sep=[4 2]; +hmri_def.coreg_bias_flags.cost_fun = 'nmi'; +hmri_def.coreg_bias_flags.fwhm = [7 7]; + +%-------------------------------------------------------------------------- +% Ordinary Least Squares & fit at TE=0 +%-------------------------------------------------------------------------- +% create an Ordinary Least Squares R2* map. The ESTATICS model is applied +% to calculate R2* map from all available contrasts. +% ADVANCED USER ONLY. +hmri_def.R2sOLS = 1; + +% Minimum number of echoes to calculate R2s map. Strictly speaking, the +% minimum is 2. For a robust estimation, the minimum number of echoes +% required depends on many factors, amongst which: +% - SNR/resolution +% - distribution/spacing between TEs: note that early echoes are more +% affected by the specific contrast, violating the assumption of a common +% decay between contrasts. +% - number of contrasts available (fewer echoes per contrast required for 3 +% (PDw, T1w, MTw) contrasts as compared to 2 or even 1) +% To be on the safe side, a minimum of 6 echoes is recommended (ESTATICS +% paper). Further studies are required to come up with more detailed and +% informed guidelines. Use fewer echoes at your own risk...! +hmri_def.neco4R2sfit = 4; + +% Define a coherent interpolation factor used all through the map creation +% process. Default is 3, but if you want to keep SNR and resolution as far +% as possible the same, it is recommended to use sinc interpolation (at +% least -4, in Siawoosh's experience -7 gives decent results). +% ADVANCED USER ONLY. +hmri_def.interp = 3; + +% Define the OLS fit as default. OLS fit at TE=0 for each contrast is used +% instead of averaged contrast images for the map calculation. +% ADVANCED USER ONLY. +hmri_def.fullOLS = true; + +%-------------------------------------------------------------------------- +% Usage of UNICORT-derived B1 maps for PD and/or MT maps calculation +% ADVANCED USER ONLY. +% WARNING: this method has not been validated for PD and MT calculation! +%-------------------------------------------------------------------------- +hmri_def.UNICORT.PD = false; +hmri_def.UNICORT.MT = false; + +%-------------------------------------------------------------------------- +% PD maps processing parameters +% ADVANCED USER ONLY. +%-------------------------------------------------------------------------- +hmri_def.PDproc.calibr = 1; % Calibration of the PD map (if PDw, T1w, + % B1map available and RF sensitivity bias correction applied somehow) + % based on PD(WM) = 69% [Tofts 2003]. +hmri_def.PDproc.WBMaskTh = 0.1; % Threshold for calculation of whole-brain mask from TPMs +hmri_def.PDproc.WMMaskTh = 0.95; % Threshold for calculation of white-matter mask from TPMs +hmri_def.PDproc.biasreg = 10^(-5); +hmri_def.PDproc.biasfwhm = 50; +hmri_def.PDproc.nr_echoes_forA = 6; % NOTE: in order to minimize R2* bias + % on the PD estimates and gain in robustness for bias-field + % correction, the number of echoes should be minimum ("average" + % calculated over the first echo only) for PD calculation. However, + % with T2*-weighting bias correction (see below), a higher number of + % echoes is preferred in order to provide good SNR. Note that when + % "fullOLS = true", this parameter has no impact whatsovever. +hmri_def.PDproc.T2scorr = 1; % to correct A map for T2*-weighting bias + % before PD map calculation. The correction is not required when + % "fullOLS = true" (TE=0 fit) and will be automatically disabled. + +%-------------------------------------------------------------------------- +% RF sensitivity bias correction: kernel for sensitivity map smoothing. +% ADVANCED USER ONLY. +%-------------------------------------------------------------------------- +hmri_def.RFsens.smooth_kernel = 12; + +%-------------------------------------------------------------------------- +% quantitative maps: quality evaluation and realignment to MNI +%-------------------------------------------------------------------------- +% creates a matlab structure containing markers of data quality +hmri_def.qMRI_maps.QA = false; +% realigns qMRI maps to MNI: the following parameter corresponds to the +% realignment implemented as part of the map calculation (see +% hmri_create_MTProt.m). Left here for backward compatibility. It is +% STRONGLY RECOMMENDED to reorient all images prior any processing using +% the Auto-Reorient module provided with the toolbox (type "help +% hmri_autoreorient" for details or open the SPM > Tools > hMRI Tools > +% Auto-Reorient module in the Batch GUI). ADVANCED USER ONLY. +hmri_def.qMRI_maps.ACPCrealign = 0; + +%-------------------------------------------------------------------------- +% Threshold values for qMRI maps +% The thresholds are meant to discard outliers generally due to low SNR in +% some brain areas, leading to physical non-sense values. Thresholding is +% required to process further the maps generated, when e.g. used +% segmentation algorithms make assumptions incompatible with existing +% outliers. +% NOTE that thresholding modifies the signal distribution and may alter +% the statistical results. +% ADVANCED USER ONLY. +%-------------------------------------------------------------------------- +hmri_def.qMRI_maps_thresh.R1 = 2000; % 1000*[s-1] +hmri_def.qMRI_maps_thresh.A = 10^5; % [a.u.] based on input images with intensities ranging approx. [0 4096]. +hmri_def.qMRI_maps_thresh.R2s = 10; % 1000*[s-1] +hmri_def.qMRI_maps_thresh.MTR = 50; +hmri_def.qMRI_maps_thresh.MTR_synt = 50; +hmri_def.qMRI_maps_thresh.MT = 5; % [p.u.] + +%-------------------------------------------------------------------------- +% MPM acquisition parameters and RF spoiling correction parameters +%-------------------------------------------------------------------------- +% ACQUISITION PARAMETERS: these values are initialised with defaults (v2k +% protocol - Prisma) and are updated at run-time with actual acquisition +% values (see hmri_create_MTProt.m). If TR/TE/FA cannot be determined from +% the input images, the following values will be used. If they don't match +% your own protocol values and if no TR/TE/FA values can be retrieved by +% the toolbox from your data, the following values should be adapted in the +% local defaults file. +% ADVANCED USER ONLY +hmri_def.MPMacq.TE_mtw = [2.34:2.34:14.04]'; %#ok<*NBRAK> % [ms] defined as column vector! +hmri_def.MPMacq.TE_t1w = [2.34:2.34:18.72]'; % [ms] +hmri_def.MPMacq.TE_pdw = [2.34:2.34:18.72]'; % [ms] +hmri_def.MPMacq.TR_mtw = 24.5; % [ms] +hmri_def.MPMacq.TR_t1w = 24.5; % [ms] +hmri_def.MPMacq.TR_pdw = 24.5; % [ms] +hmri_def.MPMacq.fa_mtw = 6; % [deg] +hmri_def.MPMacq.fa_t1w = 21; % [deg] +hmri_def.MPMacq.fa_pdw = 6; % [deg] +hmri_def.MPMacq.tag = 'v2k'; + +% IMPREFECT RF SPOILING CORRECTION +% The RF- and gradient-spoiled gradient echo sequences used to acquire the +% multiparametric mapping (MPM) data apply RF and gradient spoiling to +% destroy unwanted transverse magnetisation. Imperfect spoiling can leave +% residual bias in the apparent R1 map if no further correction is applied +% [Preibisch and Deichmann, MRM 61:125-135, 2009]. +% Correction coefficients are sequence-specific and can be determined by +% simulation [Preibisch and Deichmann 2009] to account for deviation from +% the Ernst equation. +% Since no run-time calculation of the correction parametes is currently +% performed in the toolbox, they've been calculated and stored for a +% series of standard protocols in a separate folder: +% config/local/ISC_parameters +% Additional protocol types and correction factors can be computed on-line +% (see more details below) and added to the local defaults file to apply +% proper RF correction, respectively. +% For the most standard MPM protocols using customised FLASH sequences on Siemens +% scanners, the above mentioned files folder contains files with spoiling correction coefficients +% The correction values P2_a and P2_b were obtained using the code supplied by R. Deichmann with the experimental +% parameters used for the standard MPM protocol and assuming T2 = 64 ms at 3T. +% When enabling the imperfect spoiling correction, make sure the +% coefficients used by choosing the appropriate file are definitely calculated for +% the protocol used! +% WARNING: although the TR and FA values are key parameters in simulating the +% imperfect spoiling, correction coefficients cannot be chosen based on TR +% and FA only. Amplitude and duration of gradient spoilers, phase increment +% for RF spoiling and diffusion effects must be taken into account. +% Therefore, two different sequences using identical TR anf FA are unlikely to +% use identical correction coefficients. + +% Example setting for choosing the classic FIL protocol (uncomment for usage): +% hmri_def.imperfectSpoilCorr.enabled = true; +% hmri_def.imperfectSpoilCorr.iscfile = fullfile(fileparts(fileparts(mfilename('fullpath'))),'ISC_parameters','01_classic_FIL_protocol.m'); + +% Instead, the coefficients can be calculated seperately for other protocols +% by the additional hMRI module 'Imperfect Spoiling Calc.' that efficiently +% calculates the protocol-specific correction parameters required to account for +% imperfect spoiling. The produced output file can then be set in the local +% defaults file to be incorporated in the processing. +% If not all acquisition parameters are known and vendor sequences at 3 T MRI are +% used, an alternative correction according to Baudrexel et al., 79:3082–3092, +% 2018 can be applied by specifying the RF spoiling phase increment, which +% differs per vendor. Please consult the paper for assumptions and full set +% of correction parameters. Following parameter values (in degree) are possible: +% 50 (Siemens), 117 (General Electrics), and 150 (Philips). + + +% Example setting for choosing the phase increment at Philips scanners: +hmri_def.imperfectSpoilCorr.enabled = true; +hmri_def.imperfectSpoilCorr.iscphase = 150; + +% ADVANCED USER ONLY. + + + +%========================================================================== +% Maps processing parameters +%========================================================================== + +%-------------------------------------------------------------------------- +% US segmentation parameters +%-------------------------------------------------------------------------- + +% recommended TPM for segmentation +hmri_def.proc.TPM = hmri_def.TPM ; +% Use the same as for the maps creation but one could (want to) use another +% one at some point. +% Map creation works with "standard" weighted-MR images to build the +% parametric maps. In the end these parametric maps taken together for a +% multichannel-segmention could show more details (for example subcortical +% nuclei?) and would therefore require a specific TPM. This TPM is of +% course still to be built at the moment... + +% Flags to write out posterior tissue classes in native & warped space +% - GM/WM/CSF -> write warped, mod+unmod, and native, native+dartelImp. +% - others -> nothing +hmri_def.proc.w_native = [[1 1];[1 1];[1 1];[0 0];[0 0];[0 0]]; +hmri_def.proc.w_warped = [[1 1];[1 1];[1 1];[0 0];[0 0];[0 0]]; +% Number of Gaussians per tissue class +hmri_def.proc.nGauss = [2 2 2 3 4 2]; % originally in SPM [1 1 2 3 4 2] + +end diff --git a/hmri_corr_imperf_spoil.m b/hmri_corr_imperf_spoil.m old mode 100644 new mode 100755 index 3c3a8a7d..b00081a6 --- a/hmri_corr_imperf_spoil.m +++ b/hmri_corr_imperf_spoil.m @@ -1,4 +1,4 @@ -function hmri_corr_imperf_spoil(job) +function out = hmri_corr_imperf_spoil(job) %========================================================================== % PURPOSE % Compute coefficients to correct for the effect of imperfect spoiling on @@ -130,16 +130,24 @@ function hmri_corr_imperf_spoil(job) Results.Output.RMSE_percent.T1app=round(RMSE_App,3); Results.Output.RMSE_percent.T1corr=round(RMSE_Corr,3); -Results.ToCopy{1} =['hmri_def.MPMacq_set.names{NN} = ''' job.prot_name ''';' ]; -Results.ToCopy{end+1}=['hmri_def.MPMacq_set.tags{NN} = ''' strrep(job.prot_name,' ','') ''';']; -Results.ToCopy{end+1}=['hmri_def.MPMacq_set.vals{NN} = [' num2str([TR FA]) '];']; +Results.ToCopy{1} =['hmri_def.MPMacq_set.name = ''' job.prot_name ''';' ]; +Results.ToCopy{end+1}=['hmri_def.MPMacq_set.tag = ''' strrep(job.prot_name,' ','') ''';']; +Results.ToCopy{end+1}=['hmri_def.MPMacq_set.val = [' num2str([TR FA]) '];']; Results.ToCopy{end+1}=['hmri_def.imperfectSpoilCorr.' strrep(job.prot_name,' ','') '.tag = ''' strrep(job.prot_name,' ','') ''';' ]; Results.ToCopy{end+1}=['hmri_def.imperfectSpoilCorr.' strrep(job.prot_name,' ','') '.P2_a = [' num2str(round(polyCoeffA,4)) '];']; Results.ToCopy{end+1}=['hmri_def.imperfectSpoilCorr.' strrep(job.prot_name,' ','') '.P2_b = [' num2str(round(polyCoeffB,4)) '];']; -Results.ToCopy{end+1}=['hmri_def.imperfectSpoilCorr.' strrep(job.prot_name,' ','') '.enabled = hmri_def.imperfectSpoilCorr.enabled;']; results_filename = fullfile(job.outdir,[strrep(job.prot_name,' ',''),'.json']); spm_jsonwrite(results_filename{1},Results,struct('indent','\t')); +%% 6. write m-file for faster reading of parameters +out.ISC_file = {strrep(results_filename{1},'.json','.m')}; +results_mfile = fopen(out.ISC_file{1},'w'); +fprintf(results_mfile,'global hmri_def\n'); +for n = 1:numel(Results.ToCopy) + fprintf(results_mfile,'%s\n',Results.ToCopy{n}); +end +fclose(results_mfile); + end \ No newline at end of file diff --git a/hmri_create_MTProt.m b/hmri_create_MTProt.m index f47acb62..35da0a06 100644 --- a/hmri_create_MTProt.m +++ b/hmri_create_MTProt.m @@ -250,10 +250,10 @@ Ni.dat = file_array(PT1w_forA,dm,dt, 0,1,0); create(Ni); spm_progress_bar('Init',dm(3),Ni.descrip,'planes completed'); - for p = 1:dm(3), + for p = 1:dm(3) M = spm_matrix([0 0 p]); Y = zeros(dm(1:2)); - for nr = 1:PDproc.nr_echoes_forA, + for nr = 1:PDproc.nr_echoes_forA M1 = V(nr).mat\V(1).mat*M; Y = Y + spm_slice_vol(V(nr),M1,dm(1:2),mpm_params.interp); end @@ -270,13 +270,13 @@ % NOTE: coregistration can be disabled using the hmri_def.coreg2PDw flag x_MT2PD = [0 0 0 0 0 0]; -if MTwidx +if MTwidx if mpm_params.coreg x_MT2PD = coreg_mt(Pavg{PDwidx}, Pavg{MTwidx}, mpm_params.coreg_flags); end end x_T12PD = [0 0 0 0 0 0]; -if T1widx +if T1widx if mpm_params.coreg x_T12PD = coreg_mt(Pavg{PDwidx}, Pavg{T1widx}, mpm_params.coreg_flags); coreg_mt(Pavg{PDwidx}, PT1w_forA, mpm_params.coreg_flags); @@ -348,7 +348,7 @@ W = (reg'*reg)\reg'; spm_progress_bar('Init',dm(3),'multi-contrast R2* fit','planes completed'); - for p = 1:dm(3), + for p = 1:dm(3) M = spm_matrix([0 0 p 0 0 0 1 1 1]); data = zeros([numel(TE) dm(1:2)]); @@ -472,7 +472,7 @@ W = (reg'*reg)\reg'; spm_progress_bar('Init',dm(3),'OLS R2* fit','planes completed'); - for p = 1:dm(3), + for p = 1:dm(3) M = spm_matrix([0 0 p 0 0 0 1 1 1]); data = zeros([sum(nechoes) dm(1:2)]); @@ -595,6 +595,47 @@ if MTwidx; fa_mtw_rad = fa_mtw * pi / 180; end if T1widx; fa_t1w_rad = fa_t1w * pi / 180; end +% create output files of FA correction factors for debugging +if isfield(ISC,'iscphase') + for ccon=1:mpm_params.ncon + ctag = mpm_params.input(ccon).tag; + Ni = nifti; + Ni.mat = V_pdw(1).mat; + Ni.mat0 = V_pdw(1).mat; + Ni.descrip = sprintf('Flip Angle map for %s-weighted contrast corrected for imperfect spoiling',ctag); + Ni.dat = file_array(fullfile(calcpath,[outbasename '_' ctag '_FAcorr.nii']),dm,dt, 0,1,0); + create(Ni); + FAcmap.(ctag) = Ni; + end + spm_progress_bar('Init',dm(3),'Calculating FA maps (ISC)','planes completed'); + for p = 1:dm(3) + M = M0*spm_matrix([0 0 p]); + if isempty(V_trans) + f_T = 1; + else + f_T = spm_slice_vol(V_trans(2,:),V_trans(2,:).mat\M,dm(1:2),mpm_params.interp)/100; % divide by 100, since p.u. maps + end + for ccon=1:mpm_params.ncon + ctag = mpm_params.input(ccon).tag; + eval(sprintf('f_FA_%s = f_T * fa_%sw;',ctag,lower(ctag))) + eval(sprintf('ISC_corr = zeros(size(f_FA_%s));',ctag)) + for k = 0:5 + for l = 0:5 + ISC_corr = ISC_corr + ISC.P(k+1,l+1) * eval(sprintf('f_FA_%s.^k * TR_%sw^l',ctag,lower(ctag))); + end + end + eval(sprintf('f_FA_%s = f_FA_%s .* ISC_corr * pi / 180;',ctag,ctag)) + FAcmap.(ctag).dat(:,:,p) = eval(sprintf('f_FA_%s',ctag)); + end + spm_progress_bar('Set',p); + end + for ccon=1:mpm_params.ncon + ctag = mpm_params.input(ccon).tag; + eval(strcat('V_FAcmap_',ctag,' = spm_vol(fullfile(calcpath,[''',outbasename,'_',ctag,'_FAcorr.nii'']));')) + end + spm_progress_bar('Clear'); +end + spm_progress_bar('Init',dm(3),'Calculating maps','planes completed'); % First calculate R1 & MTR @@ -637,7 +678,7 @@ % correct T1 for transmit bias f_T with fa_true = f_T * fa_nom % T1corr = T1 / f_T / f_T - if ISC.enabled + if isfield(ISC,'iscfile') % MFC: We do have P2_a and P2_b parameters for this sequence % => T1 = A(B1) + B(B1)*T1app (see Preibisch 2009) T1 = ISC.P2_a(1)*f_T.^2 + ... @@ -646,6 +687,11 @@ (ISC.P2_b(1)*f_T.^2+ISC.P2_b(2)*f_T+ISC.P2_b(3)) .* ... ((((PDw / fa_pdw_rad) - (T1w / fa_t1w_rad)+eps) ./ ... max((T1w * fa_t1w_rad / 2 / TR_t1w) - (PDw * fa_pdw_rad / 2 / TR_pdw),eps))./f_T.^2); + elseif isfield(ISC,'iscphase') + f_FA_T1 = spm_slice_vol(V_FAcmap_T1,V_FAcmap_T1.mat\M,dm(1:2),mpm_params.interp); + f_FA_PD = spm_slice_vol(V_FAcmap_PD,V_FAcmap_PD.mat\M,dm(1:2),mpm_params.interp); + T1 = ((((PDw ./ f_FA_PD) - (T1w ./ f_FA_T1)+eps) ./ ... + max((T1w .* f_FA_T1 / 2 / TR_t1w) - (PDw .* f_FA_PD/ 2 / TR_pdw),eps))); else % MFC: We do not have P2_a or P2_b parameters for this sequence % => T1 = T1app @@ -725,7 +771,12 @@ % - f_T has been calculated using UNICORT *AND* the UNICORT.PD flag % is enabled (advanced user only! method not validated yet!) if(~isempty(f_T))&&(~mpm_params.UNICORT.R1 || mpm_params.UNICORT.PD) - A = T1 .* (T1w_forA .*(fa_t1w_rad*f_T) / 2 / TR_t1w) + (T1w_forA ./ (fa_t1w_rad*f_T)); + if isfield(ISC,'iscphase') + f_FA_T1 = spm_slice_vol(V_FAcmap_T1,V_FAcmap_T1.mat\M,dm(1:2),mpm_params.interp); + A = T1 .* (T1w_forA .* f_FA_T1 / 2 / TR_t1w) + (T1w_forA ./ f_FA_T1); + else + A = T1 .* (T1w_forA .*(fa_t1w_rad*f_T) / 2 / TR_t1w) + (T1w_forA ./ (fa_t1w_rad*f_T)); + end else % semi-quantitative A A = T1 .* (T1w_forA * fa_t1w_rad / 2 / TR_t1w) + (T1w_forA / fa_t1w_rad); @@ -740,18 +791,29 @@ % and PDw ones... if MTwidx MTw = spm_slice_vol(Vavg(MTwidx),Vavg(MTwidx).mat\M,dm(1:2),3); - T1_forMT = ((PDw / fa_pdw_rad) - (T1w / fa_t1w_rad)) ./ ... - max((T1w * (fa_t1w_rad / 2 / TR_t1w)) - (PDw * fa_pdw_rad / 2 / TR_pdw),eps); - A_forMT = T1_forMT .* (T1w * fa_t1w_rad / 2 / TR_t1w) + (T1w / fa_t1w_rad); - - % MT in [p.u.]; offset by - famt * famt / 2 * 100 where MT_w = 0 (outside mask) - MT = ( (A_forMT * fa_mtw_rad - MTw) ./ (MTw+eps) ./ (T1_forMT + eps) * TR_mtw - fa_mtw_rad^2 / 2 ) * 100; + if isfield(ISC,'iscphase') + f_FA_T1 = spm_slice_vol(V_FAcmap_T1,V_FAcmap_T1.mat\M,dm(1:2),mpm_params.interp); + f_FA_PD = spm_slice_vol(V_FAcmap_PD,V_FAcmap_PD.mat\M,dm(1:2),mpm_params.interp); + f_FA_MT = spm_slice_vol(V_FAcmap_MT,V_FAcmap_MT.mat\M,dm(1:2),mpm_params.interp); + + T1_forMT = ((PDw ./ f_FA_PD) - (T1w ./ f_FA_T1)) ./ ... + max((T1w .* (f_FA_T1 / 2 / TR_t1w)) - (PDw .* f_FA_PD / 2 / TR_pdw),eps); + A_forMT = T1_forMT .* (T1w .* f_FA_T1 / 2 / TR_t1w) + (T1w ./ f_FA_T1); + MT = ( (A_forMT .* f_FA_MT - MTw) ./ (MTw+eps) ./ (T1_forMT + eps) * TR_mtw - f_FA_MT.^2 / 2 ) * 100; + else + T1_forMT = ((PDw / fa_pdw_rad) - (T1w / fa_t1w_rad)) ./ ... + max((T1w * (fa_t1w_rad / 2 / TR_t1w)) - (PDw * fa_pdw_rad / 2 / TR_pdw),eps); + A_forMT = T1_forMT .* (T1w * fa_t1w_rad / 2 / TR_t1w) + (T1w / fa_t1w_rad); + + % MT in [p.u.]; offset by - famt * famt / 2 * 100 where MT_w = 0 (outside mask) + MT = ( (A_forMT * fa_mtw_rad - MTw) ./ (MTw+eps) ./ (T1_forMT + eps) * TR_mtw - fa_mtw_rad^2 / 2 ) * 100; + 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) + if (~isempty(f_T)&&~isfield(ISC,'iscphase'))&&(~mpm_params.UNICORT.R1 || mpm_params.UNICORT.MT) MT = MT .* (1 - 0.4) ./ (1 - 0.4 * f_T); end @@ -822,7 +884,7 @@ Vsave = spm_vol(ACPC_images(i,:)); Vsave.descrip = [Vsave.descrip ' - AC-PC realigned']; spm_write_vol(Vsave,spm_read_vols(spm_vol(ACPC_images(i,:)))); - end; + end % Save transformation matrix spm_jsonwrite(fullfile(supplpath,'hMRI_map_creation_ACPCrealign_transformation_matrix.json'),R,struct('indent','\t')); @@ -840,7 +902,7 @@ % if no MT map available. Therefore, we must at least have R1 available, % i.e. both PDw and T1w inputs... if (mpm_params.QA.enable||(PDproc.calibr)) && (PDwidx && T1widx) - if ~isempty(fMT); + if ~isempty(fMT) Vsave = spm_vol(fMT); else % ~isempty(fR1); Vsave = spm_vol(fR1); @@ -1312,51 +1374,61 @@ function PDcalculation(fA, fTPM, mpm_params) mpm_params.nr_echoes4avg = min(length(find(mpm_params.input(1).TE + coregsensfnam{i} = fullfile(calcpath, tag, spm_file(output_list{1}.rfiles{1},'filename')); %#ok end coregsensfnam = coregsensfnam'; % must be a nx1 cellstr clear matlabbatch diff --git a/hmri_create_b1map.m b/hmri_create_b1map.m index 1ce2ab0d..c657e776 100755 --- a/hmri_create_b1map.m +++ b/hmri_create_b1map.m @@ -55,6 +55,10 @@ % processing B1 map from rf_map data P_trans = calc_rf_map(jobsubj, b1map_params); + case 'SDAM' + % processing B1 map from SDAM data + P_trans = calc_SDAM_b1map(jobsubj, b1map_params); + case 'pre_processed_B1' if b1map_params.scafac ~= 1 % rescale if scaling factor other than 1 is provided @@ -154,15 +158,7 @@ % save output images VB1 = V1; -% VB1.pinfo = [max(B1map(:))/16384;0;0]; -% VB1.fname = fullfile(outpath, [sname '_B1map.nii']); -% spm_write_vol(VB1,B1map); - -% VB1.pinfo = [max(B1map_norm(:))/16384;0;0]; -% VB1.fname = fullfile(outpath, [sname '_B1map_norm.nii']); -% spm_write_vol(VB1,B1map_norm); - -VB1.pinfo = [max(smB1map_norm(:))/16384;0;0]; +VB1.pinfo(1) = max(smB1map_norm(:))/spm_type(VB1.dt(1),'maxval'); VB1.descrip = 'B1+ map - smoothed and normalised (p.u.) - AFI protocol'; VB1.fname = fullfile(outpath, [sname '_B1map.nii']); spm_write_vol(VB1,smB1map_norm); @@ -187,6 +183,75 @@ end +%% =======================================================================% +% B1 map calculation - SDAM protocol % TL copied and adapted from AFI_b1map +%=========================================================================% +function P_trans = calc_SDAM_b1map(jobsubj, b1map_params) + +% default format specifications for the output metadata +json = hmri_get_defaults('json'); + +% define output dir +outpath = jobsubj.path.b1path; +b1map_params.outpath = outpath; + +% First image = alpha2 (2*alpha1) and second image = alpha1. +file_alpha1 = b1map_params.b1input(2,:); +file_alpha2 = b1map_params.b1input(1,:); +V1 = spm_vol(file_alpha1); +V2 = spm_vol(file_alpha2); +Y1 = spm_read_vols(V1); +Y2 = spm_read_vols(V2); + +alphanom = b1map_params.b1acq.alphanom; + +% Mask = squeeze(Vol1); +% threshold = (prctile(Mask(:),98)-prctile(Mask(:),2))*0.1+prctile(Mask(:),2); +% Mask = (Mask>threshold); + +B1map = (1/alphanom).*acosd(Y2./(2*Y1)); +B1map_norm = abs(B1map)*100; + +% smoothed map +smB1map_norm = zeros(size(B1map_norm)); +pxs = sqrt(sum(V1.mat(1:3,1:3).^2)); % Voxel resolution +smth = 8./pxs; +spm_smooth(B1map_norm,smB1map_norm,smth); + +% masking +% B1map = B1map.*Mask; +% B1map_norm = B1map_norm.*Mask; +% smB1map_norm = smB1map_norm.*Mask; + +sname = spm_file(V1.fname,'filename'); + +% save output images +VB1 = V1; +VB1.pinfo(1) = max(smB1map_norm(:))/spm_type(VB1.dt(1),'maxval'); +VB1.descrip = 'B1+ map - smoothed and normalised (p.u.) - SDAM protocol'; +VB1.fname = fullfile(outpath, [sname '_B1map.nii']); +spm_write_vol(VB1,smB1map_norm); + +% set and write metadata +input_files = b1map_params.b1input; +Output_hdr = init_b1_output_metadata(input_files, b1map_params); +Output_hdr.history.procstep.descrip = [Output_hdr.history.procstep.descrip ' (SDAM protocol)']; +Output_hdr.history.output.imtype = 'B1+ map (SDAM protocol)'; +set_metadata(VB1.fname,Output_hdr,json); + +% Rename anatomical reference for uniformity between protocols +B1ref = fullfile(outpath, [sname '_B1ref.nii']); +copyfile(char(file_alpha1),B1ref); +try copyfile([spm_str_manip(char(file_alpha1),'r') '.json'],[spm_str_manip(B1ref,'r') '.json']); end %#ok<*TRYNC> + +% requires anatomic image + map +P_trans = char(B1ref,char(VB1.fname)); + +% VB1.fname = fullfile(outpath, [sname '_B1map_mask.nii']); +% spm_write_vol(VB1,Mask); + +end + %% =======================================================================% % B1 map calculation - SE/STE EPI protocol %=========================================================================% @@ -463,7 +528,7 @@ sname = spm_file(V1.fname,'basename'); VB1 = V1; -VB1.pinfo = [max(smB1map_norm(:))/16384;0;0]; % what is this for? (TL) +VB1.pinfo(1) = max(smB1map_norm(:))/spm_type(VB1.dt(1),'maxval'); VB1.fname = fullfile(outpath, [sname '_B1map.nii']); VB1.descrip = 'Smoothed & normalised (p.u.) B1 bias map - TFL B1map protocol'; spm_write_vol(VB1,smB1map_norm); @@ -522,7 +587,7 @@ sname = spm_file(V1.fname,'basename'); VB1 = V1; -VB1.pinfo = [max(smB1map_norm(:))/16384;0;0]; % what is this for? (TL) +VB1.pinfo(1) = max(smB1map_norm(:))/spm_type(VB1.dt(1),'maxval'); VB1.fname = fullfile(outpath, [sname '_B1map.nii']); VB1.descrip = 'Smoothed & normalised (p.u.) B1 bias map - TFL B1map protocol'; spm_write_vol(VB1,smB1map_norm); @@ -788,6 +853,26 @@ end end + case 'SDAM' + if ~isempty(b1map_params.b1input) + hmri_log(sprintf('SDAM protocol selected ...'),b1map_params.nopuflags); + b1hdr = get_metadata(b1map_params.b1input(1,:)); + + try + fa = get_metadata_val(b1hdr{2},'FlipAngle'); + if isempty(fa) + hmri_log(sprintf('WARNING: using defaults values for flip angle\n(FA = %d deg) instead of metadata', ... + b1map_params.b1acq.alphanom),b1map_params.defflags); + else + b1map_params.b1acq.alphanom = fa; + end + + catch + hmri_log(sprintf(['WARNING: possibly no metadata associated to the input images. \n' ... + 'Default acquisition and processing parameters will be used.']),b1map_params.defflags); + end + end + case 'tfl_b1_map' if ~isempty(b1map_params.b1input) hmri_log(sprintf('SIEMENS tfl_b1map protocol selected ...'),b1map_params.nopuflags); diff --git a/hmri_quality_display.m b/hmri_quality_display.m new file mode 100644 index 00000000..b3e25438 --- /dev/null +++ b/hmri_quality_display.m @@ -0,0 +1,129 @@ +function h = hmri_quality_display(disp_list) +%========================================================================== +% PURPOSE: to display a series of volumes using spm_check_registration. +% USAGE: h = hmri_quality_display(disp_list) +% where +% - disp_list is an array of structures with fields: +% - fnam: name of the image file name +% - title: title for the image +% - range: intensity range [min max] +% - h is the figure handle +%========================================================================== +% Written by Evelyne Balteau +% Cyclotron Research Centre, University of Liege +%========================================================================== + +global st; +% st is a global variable storing the display parameters. +% st is a structure with fields: +% n: 0 +% vols: {24x1 cell} +% bb: [2x3 double] +% Space: [4x4 double] +% centre: [-2.4096 30.6205 -6.1265] +% callback: ';' +% xhairs: 1 +% hld: 1 +% fig: [1x1 Figure] +% mode: 1 +% plugins: {'browser' 'contour' 'display' 'goto_max' 'mesh' 'movie' 'reorient' 'rgb' 'roi' 'save' 'dwprofile' 'pli' 'quiver' 'quiver3d' 'tracepath'} +% snap: [] +% st.vols is a cell array whose elements (st.vols{i}) are structures with fields: +% fname: 'D:\home\data\truc.nii' +% dim: [240 256 176] +% dt: [512 0] +% pinfo: [3x1 double] +% mat: [4x4 double] +% n: [1 1] +% descrip: 'NIfTI description field' +% private: [1x1 nifti] +% ax: {3x1 cell} +% premul: [4x4 double] +% window: 'auto' +% mapping: 'linear' +% area: [0.0100 0.5100 0.4800 0.4800] + + + +% gather all file names +fnamlist = char(disp_list.fnam); +% display all images +spm_check_registration(fnamlist); +% to get the section shifted in x-direction +% (to avoid displaying the plane between left and right hemispheres) +% NB: depending on the history of manipulating st.centre, the latter can +% either be a column or a line vector. So... +st.centre(1) = st.centre(1) + 10; +st.centre(2) = st.centre(2) + 0; +st.centre(3) = st.centre(3) + 0; +% st.centre = st.centre + [10 0 0]; + +% to zoom in +%st.Space = [0.25 0 0 0;0 0.25 0 0;0 0 0.25 0;0 0 0 1]; + +% to give a title to each orthviews and adjust intensity range: +for cim=1:length(disp_list) + txt = disp_list(cim).title; + % rearrange title if too long + % (split over several lines of maximum maxchar characters length) + maxchar = 20; + if length(txt)>maxchar + try + txtsplit = strsplit(txt); % problem when using Matlab < R2013a (when strsplit and strjoin were added) + catch %#ok<*CTCH> + txtsplit = regexp(txt,regexptranslate('escape',' '),'split'); + end + clear linsplit; + if length(txtsplit)> 1 % text naturally splittable + linsplit{1} = txtsplit{1}; + cw = 2; + cl = 1; + while cw + else + cl = cl+1; + linsplit{cl} = txtsplit{cw}; + end + cw = cw+1; + end + else % force split + cl = 1; + while length(txt)>maxchar + linsplit{cl} = txt(1:maxchar); + txt = txt(maxchar+1:end); + cl = cl+1; + end + linsplit{cl} = txt; + end + + try + txt = strjoin(linsplit,'\n'); % problem when using Matlab < R2013a (when strsplit and strjoin were added) + catch + txt = linsplit{1}; + for cc = 2:length(linsplit) + txt = sprintf('%s\n%s', txt, linsplit{cc}); + end + end + end + + %htit = get(st.vols{cim}.ax{3}.ax,'Title'); + hxlab = get(st.vols{cim}.ax{3}.ax,'Xlabel'); + if any( disp_list(cim).range - round( disp_list(cim).range)) + set(hxlab, 'String',sprintf('%s\n\nIntensity range\n[%.1f %.1f]', txt, disp_list(cim).range),'FontSize',10); + % set(st.vols{cim}.ax{3}.ax, 'Xlabel',text('String',sprintf('Intensity range\n[%.1f %.1f]', disp_list(cim).range))); + else + set(hxlab, 'String',sprintf('%s\n\nIntensity range\n[%d %d]', txt, disp_list(cim).range),'FontSize',10); + % set(st.vols{cim}.ax{3}.ax, 'Xlabel',text('String',sprintf('Intensity range\n[%d %d]', disp_list(cim).range))); + end + %set(htit, 'String',sprintf('%s',disp_list(cim).title)); + % set(st.vols{cim}.ax{3}.ax, 'Title',text('String',sprintf('%s',disp_list(cim).title))); + st.vols{cim}.window = disp_list(cim).range; +end + +spm_orthviews('Redraw'); + +h = st.fig; + +end \ No newline at end of file diff --git a/tbx_cfg_hmri.m b/tbx_cfg_hmri.m old mode 100644 new mode 100755 index 6b87a5d5..3145926d --- a/tbx_cfg_hmri.m +++ b/tbx_cfg_hmri.m @@ -43,6 +43,6 @@ 'and will include a number of (as yet unspecified) extensions in ',... 'future updates. Please report any bugs or problems to lreninfo@gmail.com.'] }'; -hmri.values = {tbx_scfg_hmri_config tbx_scfg_hmri_dicom_import tbx_scfg_hmri_autoreorient tbx_scfg_hmri_imperf_spoil tbx_scfg_hmri_create tbx_scfg_hmri_proc}; +hmri.values = {tbx_scfg_hmri_config tbx_scfg_hmri_dicom_import tbx_scfg_hmri_autoreorient tbx_scfg_hmri_imperf_spoil tbx_scfg_hmri_create tbx_scfg_hmri_quality tbx_scfg_hmri_proc}; end diff --git a/tbx_scfg_hmri_create.m b/tbx_scfg_hmri_create.m index b5bfc567..140e99b3 100644 --- a/tbx_scfg_hmri_create.m +++ b/tbx_scfg_hmri_create.m @@ -222,6 +222,20 @@ b1_input_tfl.val = {b1raw}; +% --------------------------------------------------------------------- +% SDAM B1 protocol +% --------------------------------------------------------------------- +b1_input_SDAM = cfg_branch; +b1_input_SDAM.tag = 'SDAM'; +b1_input_SDAM.name = 'Saturated Double Angle Method'; +b1_input_SDAM.help = {'Saturated Double Angle Method (SDAM) protocol.', ... + 'As B1 input, please select a 2*alpha/alpha (e.g. 120°/60°) pair of images in that order.', ... + ['Regarding processing parameters, you can either stick with metadata and standard ' ... + 'defaults parameters (recommended) or select your own [hmri_b1_local_defaults_*.m] customised defaults file ' ... + '(fallback for situations where no metadata are available).']}; +b1_input_SDAM.val = {b1raw b1parameters}; + + % --------------------------------------------------------------------- % i3D_AFI B1 protocol % --------------------------------------------------------------------- @@ -271,6 +285,7 @@ 'PLoS One 2012;7(3):e32379].'] [' - 3D AFI: 3D actual flip angle imaging (AFI) method based on [Yarnykh VL, ' ... 'Magn Reson Med 2007;57:192-200].'] + [' - Saturated Double Angle Method (SDAM)'] [' - tfl_b1_map: Siemens product sequence for B1 mapping based on turbo FLASH.'] [' - rf_map: Siemens product sequence for B1 mapping based on SE/STE.'] [' - no B1 correction: if selected no B1 bias correction will be applied.'] @@ -281,7 +296,7 @@ 'using the approach described in [Weiskopf et al., NeuroImage 2011; 54:2116-2124]. ' ... 'WARNING: the correction only applies to R1 maps.'] }; %#ok<*NBRAK> -b1_type.values = {b1_input_3DEPI b1_input_3DAFI b1_input_tfl b1_input_rfmap b1_input_preproc b1_input_UNICORT b1_input_noB1}; +b1_type.values = {b1_input_3DEPI b1_input_3DAFI b1_input_SDAM b1_input_tfl b1_input_rfmap b1_input_preproc b1_input_UNICORT b1_input_noB1}; b1_type.val = {b1_input_3DEPI}; % --------------------------------------------------------------------- @@ -513,4 +528,4 @@ dep = cdep; end -%_______________________________________________________________________ \ No newline at end of file +%_______________________________________________________________________ diff --git a/tbx_scfg_hmri_imperf_spoil.m b/tbx_scfg_hmri_imperf_spoil.m index 49d25924..ca11d701 100644 --- a/tbx_scfg_hmri_imperf_spoil.m +++ b/tbx_scfg_hmri_imperf_spoil.m @@ -10,7 +10,7 @@ % %_______________________________________________________________________ % Wellcome Centre for Human Neuroimaging -% Nadège Corbin - May 2020 +% Nad�ge Corbin - May 2020 % ====================================================================== % --------------------------------------------------------------------- @@ -174,3 +174,4 @@ 'this module computes coefficients required to correct for imperfect spoiling in the FLASH volumes ' ... 'using the method proposed by Preibisch & Deichmann, MRM 2009, 61(1):125'}; imperf_spoil.prog = @hmri_corr_imperf_spoil; +end diff --git a/tbx_scfg_hmri_quality.m b/tbx_scfg_hmri_quality.m new file mode 100755 index 00000000..be39082e --- /dev/null +++ b/tbx_scfg_hmri_quality.m @@ -0,0 +1,246 @@ +function quality = tbx_scfg_hmri_quality +%========================================================================== +% Configuration (sub)file for the hMRI-toolbox +% PURPOSE +% Implementation of a series of (basic) quality assessment tools, +% mainly based on visual inspection of results, as a first step. +% +% LIST OF TOOLS +% - check reorientation: display (in CheckReg) the reoriented image +% (preferably the one used for reorientation) and the template used for +% reorientation -> saved as PNG file for quick visual inspection. +% - display all maps: with standard widowing, and saved as PNG file for +% quick wisual inspection. +% - display input images: per image type (e.g. all B1 input images, of all +% FLASH images, etc...), adjust windowing to visualise all images with +% identical windowing. Visual inspection, search for obvious motion +% artefacts, ... +% +% WARNING +% So far, these tools are not meant to automatically flag anomalies. The +% assessment still relies on the expertise of the person visualizing the +% data! Moreover, the PNG files saved don't show everything, and anomalies +% can still keep unnoticed... +%========================================================================== +% Written by Evelyne Balteau +% Cyclotron Research Centre, University of Liege +%========================================================================== + +%-------------------------------------------------------------------------- +% res_dir Output directory +%-------------------------------------------------------------------------- +res_dir = cfg_files; +res_dir.tag = 'res_dir'; +res_dir.name = 'Output directory'; +res_dir.val{1} = {''}; +res_dir.help = {['Files produced by this function will be written into ' ... + 'this output directory. If no directory is given, results images will ' ... + 'be written to the current working directory.']}; +res_dir.filter = 'dir'; +res_dir.ufilter = '.*'; +res_dir.num = [0 1]; + +% --------------------------------------------------------------------- +% check input images +% --------------------------------------------------------------------- +input_images = cfg_files; +input_images.tag = 'input_images'; +input_images.name = 'Input images'; +input_images.help = {'Input images.' + 'Select images according to the description of the specific quality check you are about to run. Maximum number of images = 15.'}; +input_images.filter = 'image'; +input_images.ufilter = '.*'; +input_images.num = [0 15]; +input_images.val = {''}; + +% --------------------------------------------------------------------- +% check_type - Menu listing the types of visual check available +% --------------------------------------------------------------------- +check_type = cfg_menu; +check_type.tag = 'check_type'; +check_type.name = 'Quality check'; +check_type.help = {'Specify the type of quality check you want to run. ' + 'You can select either:' + ['- Input images: Select images all from a single image type ' ... + '(e.g. all B1 input images, or all spoiled gradient echo (FLASH) images, etc...). ' ... + 'The images will be displayed and the intensity scaling will be ' ... + 'automatically adjusted to visualise all images with ' ... + 'identical intensity scaling. Purpose: visual inspection, search for obvious motion ' ... + 'artefacts, ...'] + ['- Check (re)orientation: Select images that are expected to be realigned ' ... + '(e.g. FLASH image and template used for auto-reorientation). ' ... + 'Images are displayed with intensity scaling adjusted for each image.'] + ['- qMRI maps: Select B1+ maps, RF sensitivity maps, R1, R2*, MT or PD maps. ' ... + 'The intensity scaling will automatically adjust to use standard and ' ... + 'comparable scaling for each qMRI map. ']}; +check_type.labels = {'Input images' + 'Check (re)orientation' + 'qMRI maps'}'; +check_type.values = {'check_inputs' + 'check_reorient' + 'check_maps'}'; +check_type.val = {'check_inputs'}; + +% --------------------------------------------------------------------- +% visual_check - Visual inspection of various type of data +% --------------------------------------------------------------------- +visual_check = cfg_exbranch; +visual_check.tag = 'visual_check'; +visual_check.name = 'Visual check'; +visual_check.val = { check_type input_images res_dir }; +visual_check.help = {'Visual inspection of the data.' + ['Select type of data, images (according to data type) and output directory ' ... + 'where the results will be saved. A PNG image is saved for later review.']}; +visual_check.prog = @hmri_quality_visual_check; +% visual_check.vout = @vout_quality_visual_check; % not implemented + +% --------------------------------------------------------------------- +% quality - list of quality tools +% --------------------------------------------------------------------- +quality = cfg_choice; +quality.tag = 'quality'; +quality.name = 'Quality tools'; +quality.help = { + ['Implementation of a series of (basic) quality assessment tools, ' ... + 'mainly based on visual inspection of results, as a first step. ']}; +quality.values = {visual_check}; + +end + + +%========================================================================== +% For all visual check cases... +%========================================================================== +function hmri_quality_visual_check(job) +% job is a structure with fields: +% check_type: 'check_inputs'/'check_reorient'/'check_maps' +% input_images: {filenames} +% res_dir: {dirname} + +% flags for hmri_log +log_flags = struct('ComWin',true,'PopUp',false,'LogFile',struct('Enabled',false)); + +if isempty(job.res_dir) + job.res_dir = {pwd}; +end +job.res_dir = char(job.res_dir); + +for cim = 1:length(job.input_images) + Y{cim} = spm_read_vols(spm_vol(char(job.input_images{cim}))); %#ok<*AGROW> +end + +switch(job.check_type) + case 'check_inputs' + % define min and max threshold based on percentiles of the ensemble + % of input images for identical intensity scaling for all images + tmp = []; + for cim = 1:length(Y) + tmp = [tmp; Y{cim}(:)]; + end + mini = prctile(tmp,5); + maxi = prctile(tmp,95); +% threshold = (myprctile(mask(:),98)-myprctile(mask(:),2))*0.02+myprctile(mask(:),2); + for cim=1:length(Y) + disp_list(cim).fnam = char(job.input_images{cim}); + disp_list(cim).title = spm_file(disp_list(cim).fnam,'basename'); + disp_list(cim).title(strcmp(disp_list(cim).title, '_')) = ' '; + disp_list(cim).range = [mini maxi]; + end + + case 'check_reorient' + % define min and max threshold based on percentiles for each image + % separately + for cim = 1:length(Y) + mini = prctile(Y{cim}(:),5); + maxi = prctile(Y{cim}(:),95); + disp_list(cim).fnam = char(job.input_images{cim}); + disp_list(cim).title = spm_file(disp_list(cim).fnam,'basename'); + disp_list(cim).title(strcmp(disp_list(cim).title, '_')) = ' '; + disp_list(cim).range = [mini maxi]; + end + + case 'check_maps' + % define min and max threshold based on data type: + % - B1: [75-125] p.u. + % - RF sensitivity: [0 1] + % - R2*: [0 70] s-1 + % - R1: [0 1.4] s-1 + % - MT saturation: [0 2] p.u. + % - PD: [50 120] p.u. + for cim = 1:length(Y) + notok = false; + disp_list(cim).fnam = char(job.input_images{cim}); + tmphdr = get_metadata(disp_list(cim).fnam); + if isempty(tmphdr) + hmri_log(sprintf(['WARNING (hmri_quality): No metadata associated with qMRI map.' ... + '\n%s \nImage displayed with automatic (non-standard) intensity scaling.'], ... + disp_list(cim).fnam), log_flags); + notok = true; + end + try + tmp = tmphdr{1}.history.output.imtype; + if iscell(tmp) + disp_list(cim).title = tmphdr{1}.history.output.imtype{1}; + else + disp_list(cim).title = tmphdr{1}.history.output.imtype; + end + catch %#ok + hmri_log(sprintf(['WARNING (hmri_quality): No imtype defined for metadata associated with qMRI map.' ... + '\n%s \nImage displayed with automatic (non-standard) intensity scaling.'], ... + disp_list(cim).fnam), log_flags); + notok = true; + end + if notok + mini = prctile(Y{cim}(:),5); + maxi = prctile(Y{cim}(:),95); + disp_list(cim).title = spm_file(disp_list(cim).fnam,'basename'); + disp_list(cim).title(strcmp(disp_list(cim).title, '_')) = ' '; + disp_list(cim).range = [mini maxi]; + else + if strfind(disp_list(cim).title,'MT') % strfind is case sensitive! + disp_list(cim).range = [0 2]; + elseif strfind(disp_list(cim).title,'PD') + disp_list(cim).range = [50 120]; + elseif strfind(disp_list(cim).title,'R1') + disp_list(cim).range = [0 1.4]; + elseif strfind(disp_list(cim).title,'R2') + disp_list(cim).range = [0 70]; + elseif strfind(disp_list(cim).title,'B1') + disp_list(cim).range = [75 125]; + elseif strfind(disp_list(cim).title,'RF sensitivity map') + disp_list(cim).range = [0 2]; + else % includes case of A maps + hmri_log(sprintf(['WARNING (hmri_quality): Non standard scaling for qMRI map.' ... + '\n%s \nImage displayed with automatic (non-standard) intensity scaling.'], ... + disp_list(cim).fnam), log_flags); + mini = prctile(Y{cim}(:),5); + maxi = prctile(Y{cim}(:),95); + disp_list(cim).range = [mini maxi]; + end + end + end + + otherwise + error('Unknown visual check'); +end + +if ~isempty(disp_list) + h = hmri_quality_display(disp_list); + + colormap(h,'gray'); + spm_orthviews('Redraw'); + res = 1; + savenam = fullfile(job.res_dir, [job.check_type sprintf('_%0.3d.png', res)]); + while exist(spm_file(savenam,'suffix','_gray'),'file') + res = res+1; + savenam = fullfile(job.res_dir, [job.check_type sprintf('_%0.3d.png', res)]); + end + set_metadata(strrep(savenam,'png','json'),disp_list); + print(h,'-dpng',spm_file(savenam,'suffix','_gray')); + + colormap(h,'jet'); + spm_orthviews('Redraw'); + print(h,'-dpng',spm_file(savenam,'suffix','_jet')); +end + +end \ No newline at end of file