From 9bc846efd5d218d72ecec9c9af3cf595a148e223 Mon Sep 17 00:00:00 2001 From: ebalteau Date: Tue, 31 Jul 2018 23:40:30 +0200 Subject: [PATCH 01/22] Module for visual inspection of the results Three types of visualisation are possible: - automated intensity scaling for all images (to check e.g. reorientation comparing input image and template) - identical intensity scaling for all images (to check e.g. series of echoes from a T1w acquisition) - standard scaling as defined for each output map (R1, R2*, MT, PD, B1, RFsens). Scaling is automated (between percentiles 5 and 95) if image type cannot be determined. Allows for quick visual check of any obvious anomaly in the output maps. --- hmri_quality_display.m | 74 +++++++++++++ tbx_cfg_hmri.m | 2 +- tbx_scfg_hmri_quality.m | 240 ++++++++++++++++++++++++++++++++++++++++ 3 files changed, 315 insertions(+), 1 deletion(-) create mode 100644 hmri_quality_display.m create mode 100644 tbx_scfg_hmri_quality.m diff --git a/hmri_quality_display.m b/hmri_quality_display.m new file mode 100644 index 00000000..e8b753a7 --- /dev/null +++ b/hmri_quality_display.m @@ -0,0 +1,74 @@ +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) +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 intesity range: +for cim=1:length(disp_list) + if any( disp_list(cim).range - round( disp_list(cim).range)) + set(st.vols{cim}.ax{3}.ax, 'Xlabel',text('String',sprintf('Intensity range\n[%.1f %.1f]', disp_list(cim).range))); + else + set(st.vols{cim}.ax{3}.ax, 'Xlabel',text('String',sprintf('Intensity range\n[%d %d]', disp_list(cim).range))); + end + 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 index 1078cece..ad888f9c 100644 --- a/tbx_cfg_hmri.m +++ b/tbx_cfg_hmri.m @@ -31,6 +31,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_create tbx_scfg_hmri_proc }; +hmri.values = {tbx_scfg_hmri_config tbx_scfg_hmri_dicom_import tbx_scfg_hmri_autoreorient tbx_scfg_hmri_create tbx_scfg_hmri_quality tbx_scfg_hmri_proc }; end diff --git a/tbx_scfg_hmri_quality.m b/tbx_scfg_hmri_quality.m new file mode 100644 index 00000000..f82626b5 --- /dev/null +++ b/tbx_scfg_hmri_quality.m @@ -0,0 +1,240 @@ +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 + disp_list(cim).title = tmphdr{1}.history.output.imtype{1}; + catch + 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 1]; + 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 + 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 From 5697df46402214265a34846496310ef6db1ffcdd Mon Sep 17 00:00:00 2001 From: ebalteau Date: Wed, 1 Aug 2018 00:52:16 +0200 Subject: [PATCH 02/22] Few adjustments - range for RF sensitivity maps [0 2] - backward/forward compatibility between very last version of SPM and previous ones - to handle B1+ and RFsens maps imtype --- hmri_quality_display.m | 11 ++++++++--- tbx_scfg_hmri_quality.m | 11 ++++++++--- 2 files changed, 16 insertions(+), 6 deletions(-) diff --git a/hmri_quality_display.m b/hmri_quality_display.m index e8b753a7..8ed57a5e 100644 --- a/hmri_quality_display.m +++ b/hmri_quality_display.m @@ -58,12 +58,17 @@ % to give a title to each orthviews and adjust intesity range: for cim=1:length(disp_list) + 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(st.vols{cim}.ax{3}.ax, 'Xlabel',text('String',sprintf('Intensity range\n[%.1f %.1f]', disp_list(cim).range))); + set(hxlab, 'String',sprintf('Intensity range\n[%.1f %.1f]', disp_list(cim).range)); + % set(st.vols{cim}.ax{3}.ax, 'Xlabel',text('String',sprintf('Intensity range\n[%.1f %.1f]', disp_list(cim).range))); else - set(st.vols{cim}.ax{3}.ax, 'Xlabel',text('String',sprintf('Intensity range\n[%d %d]', disp_list(cim).range))); + set(hxlab, 'String',sprintf('Intensity range\n[%d %d]', disp_list(cim).range)); + % set(st.vols{cim}.ax{3}.ax, 'Xlabel',text('String',sprintf('Intensity range\n[%d %d]', disp_list(cim).range))); end - set(st.vols{cim}.ax{3}.ax, 'Title',text('String',sprintf('%s',disp_list(cim).title))); + 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 diff --git a/tbx_scfg_hmri_quality.m b/tbx_scfg_hmri_quality.m index f82626b5..14c6a6b6 100644 --- a/tbx_scfg_hmri_quality.m +++ b/tbx_scfg_hmri_quality.m @@ -178,8 +178,13 @@ function hmri_quality_visual_check(job) notok = true; end try - disp_list(cim).title = tmphdr{1}.history.output.imtype{1}; - catch + 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); @@ -203,7 +208,7 @@ function hmri_quality_visual_check(job) 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 1]; + 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.'], ... From 6e64805c7250b4b7e9769f43d949b7bc2b54a68a Mon Sep 17 00:00:00 2001 From: ebalteau Date: Wed, 1 Aug 2018 14:42:37 +0200 Subject: [PATCH 03/22] Visual quality check continued... Bit of cosmetic changes... --- hmri_quality_display.m | 44 +++++++++++++++++++++++++++++++++++++----- 1 file changed, 39 insertions(+), 5 deletions(-) diff --git a/hmri_quality_display.m b/hmri_quality_display.m index 8ed57a5e..c3edbd43 100644 --- a/hmri_quality_display.m +++ b/hmri_quality_display.m @@ -56,18 +56,52 @@ % 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 intesity range: +% to give a title to each orthviews and adjust intensity range: for cim=1:length(disp_list) - htit = get(st.vols{cim}.ax{3}.ax,'Title'); + 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 + txtsplit = strsplit(txt); + 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 + + txt = strjoin(linsplit,'\n'); + 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('Intensity range\n[%.1f %.1f]', 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('Intensity range\n[%d %d]', disp_list(cim).range)); + 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(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 From bc7a507dcc75bc5737bebac7c3bfaf3da164a7f5 Mon Sep 17 00:00:00 2001 From: ebalteau Date: Sat, 25 Aug 2018 23:00:44 +0200 Subject: [PATCH 04/22] Bug fixes and backward compatibility - Depending on the history of manipulating the st.centre field, the latter can either be a column or a line vector. Elements processed one by one to avoid problems. - strsplit and strjoin were added to Matlab from version R2013a. Alternative implementation for that specific case. --- hmri_quality_display.m | 22 +++++++++++++++++++--- 1 file changed, 19 insertions(+), 3 deletions(-) diff --git a/hmri_quality_display.m b/hmri_quality_display.m index c3edbd43..b3e25438 100644 --- a/hmri_quality_display.m +++ b/hmri_quality_display.m @@ -51,7 +51,12 @@ spm_check_registration(fnamlist); % to get the section shifted in x-direction % (to avoid displaying the plane between left and right hemispheres) -st.centre = st.centre + [10 0 0]; +% 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]; @@ -63,7 +68,11 @@ % (split over several lines of maximum maxchar characters length) maxchar = 20; if length(txt)>maxchar - txtsplit = strsplit(txt); + 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}; @@ -89,7 +98,14 @@ linsplit{cl} = txt; end - txt = strjoin(linsplit,'\n'); + 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'); From 83c3182f28887fa3c35ec9eea5741f50df39cad5 Mon Sep 17 00:00:00 2001 From: ebalteau Date: Wed, 3 Oct 2018 07:41:46 +0200 Subject: [PATCH 05/22] Batch GUI - Imperfect spoiling correction A start... --- tbx_scfg_hmri_create.m | 44 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 44 insertions(+) diff --git a/tbx_scfg_hmri_create.m b/tbx_scfg_hmri_create.m index b5bfc567..687c1649 100644 --- a/tbx_scfg_hmri_create.m +++ b/tbx_scfg_hmri_create.m @@ -7,6 +7,50 @@ % Christophe Phillips + +%-------------------------------------------------------------------------- +% Imperfect spoiling correction +%-------------------------------------------------------------------------- +isc = cfg_menu; +isc.tag = 'isc'; +isc.name = 'Imperfect spoiling correction'; +isc.help = {['The RF- and gradient-spoiled gradient echo sequences used to ' ... + 'acquire the multiparametric mapping (MPM) data use 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 used ' ... + '[Preibisch and Deichmann 2009]. For specific MPM protocols using the customised ' ... + 'sequences, the hMRI-toolbox provides spoiling correction coefficients to ' ... + 'account for deviation from the Ernst equation. The coefficients are ' ... + 'sequence-specific and can be determined by simulation [Preibisch and Deichmann 2009].'], ... + 'The following options are available:', ... + '- None: by default the correction is disabled. ', ... + '- Standard: if your MPM protocol is one of the standard, ' ... + 'predefined ones, the correction parameters can be automaticallyselect the "Standard" option. ' ... + 'You can also provide the toolbox with a customized set of correction parametersSelect the "Customized" ](DefaultsAndCustomization). An additional module that efficiently calculates the protocol-specific correction parameters required to account for imperfect spoiling is planned. + + + + + +The user can review and keep track of all the information ' ... + 'collected, including warnings and other messages coming up during ' ... + 'the creation of the maps. By default, the information is logged in ' ... + 'the Matlab Command Window, in a log file saved in the "Results/Supplementary" ' ... + 'directory, and when more critical, displayed as a pop-up message.'], ... + ['The latter must be disabled for processing series of datasets (since it ' ... + 'blocks the execution of the code) but it is strongly recommended to ' ... + 'leave it enabled when piloting the data processing (single subject) ' ... + 'to read through and acknowledge every message and make sure ' ... + 'everything is set up properly before running the processing on a ' ... + 'whole group.'], ... + ['More information about the various messages and action to be taken ' ... + '(or not) accordingly can be found on the hMRI-Toolbox WIKI (http://hmri.info). ' ... + 'In particular, see the "Debug tips & tricks" section.']}; +isc.labels = {'Disable' 'Enable'}; +isc.values = {false true}; +isc.val = {true}; + + %-------------------------------------------------------------------------- % To enable/disable pop-up messages for all warnings - recommended when % piloting the data processing. From 6b97a61c9cbe4a2769f2ae166b91274e68fed6a7 Mon Sep 17 00:00:00 2001 From: ebalteau Date: Tue, 9 Oct 2018 17:48:59 +0200 Subject: [PATCH 06/22] ISC - batch GUI continued... Two options can be selected: no correction (default) or select a json file containing the correction coefficients. --- tbx_scfg_hmri_create.m | 88 ++++++++++++++++++++++++------------------ 1 file changed, 51 insertions(+), 37 deletions(-) diff --git a/tbx_scfg_hmri_create.m b/tbx_scfg_hmri_create.m index 687c1649..77829f88 100644 --- a/tbx_scfg_hmri_create.m +++ b/tbx_scfg_hmri_create.m @@ -7,48 +7,62 @@ % Christophe Phillips +%-------------------------------------------------------------------------- +% Imperfect spoiling correction - no correction +%-------------------------------------------------------------------------- +iscnone = cfg_entry; +iscnone.tag = 'iscnone'; +iscnone.name = 'None'; +iscnone.help = {'No correction will be applied.'}; +iscnone.strtype = 's'; +iscnone.num = [1 Inf]; +iscnone.val = {'none'}; + +%-------------------------------------------------------------------------- +% Imperfect spoiling correction - select coefficients file +%-------------------------------------------------------------------------- +iscfile = cfg_files; +iscfile.tag = 'iscfile'; +iscfile.name = 'Select ISC file'; +iscfile.help = {['Select a *.json file containing the correction coefficients ' ... + 'corresponding to the current acquisition protocol.']}; +iscfile.filter = 'json'; +iscfile.dir = fullfile(fileparts(mfilename('fullpath')),'config'); +iscfile.ufilter = '.*'; +iscfile.num = [1 1]; %-------------------------------------------------------------------------- -% Imperfect spoiling correction +% Imperfect spoiling correction - select coefficients %-------------------------------------------------------------------------- -isc = cfg_menu; +isc = cfg_choice; isc.tag = 'isc'; -isc.name = 'Imperfect spoiling correction'; -isc.help = {['The RF- and gradient-spoiled gradient echo sequences used to ' ... - 'acquire the multiparametric mapping (MPM) data use RF and gradient spoiling to ' ... +isc.name = 'Imperfect spoiling correction (ISC)'; +isc.help = {['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 used ' ... - '[Preibisch and Deichmann 2009]. For specific MPM protocols using the customised ' ... - 'sequences, the hMRI-toolbox provides spoiling correction coefficients to ' ... - 'account for deviation from the Ernst equation. The coefficients are ' ... - 'sequence-specific and can be determined by simulation [Preibisch and Deichmann 2009].'], ... + 'residual bias in the apparent R1 map if no further correction is applied ' ... + '[Preibisch and Deichmann 2009]. Correction coefficients are sequence-' ... + 'specific and can be determined by simulation [Preibisch and Deichmann 2009] to ' ... + 'account for deviation from the Ernst equation. '],[''], ... 'The following options are available:', ... - '- None: by default the correction is disabled. ', ... - '- Standard: if your MPM protocol is one of the standard, ' ... - 'predefined ones, the correction parameters can be automaticallyselect the "Standard" option. ' ... - 'You can also provide the toolbox with a customized set of correction parametersSelect the "Customized" ](DefaultsAndCustomization). An additional module that efficiently calculates the protocol-specific correction parameters required to account for imperfect spoiling is planned. - - - - - -The user can review and keep track of all the information ' ... - 'collected, including warnings and other messages coming up during ' ... - 'the creation of the maps. By default, the information is logged in ' ... - 'the Matlab Command Window, in a log file saved in the "Results/Supplementary" ' ... - 'directory, and when more critical, displayed as a pop-up message.'], ... - ['The latter must be disabled for processing series of datasets (since it ' ... - 'blocks the execution of the code) but it is strongly recommended to ' ... - 'leave it enabled when piloting the data processing (single subject) ' ... - 'to read through and acknowledge every message and make sure ' ... - 'everything is set up properly before running the processing on a ' ... - 'whole group.'], ... - ['More information about the various messages and action to be taken ' ... - '(or not) accordingly can be found on the hMRI-Toolbox WIKI (http://hmri.info). ' ... - 'In particular, see the "Debug tips & tricks" section.']}; -isc.labels = {'Disable' 'Enable'}; -isc.values = {false true}; -isc.val = {true}; + '- None [default]: no correction applied. ', ... + ['- Select ISC file: you must choose a set of correction coefficients. ' ... + 'For the most standard MPM protocols using customised FLASH sequences on Siemens ' ... + 'scanners, the hMRI-toolbox provides spoiling correction coefficients. ' ... + 'For other protocols, the coefficients can be calculated seperately ' ... + 'and loaded here.'],[''], ... + ['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.'],[''], ... + ['ISC-MODULE: an additional hMRI module that efficiently calculates the ' ... + 'protocol-specific correction parameters required to account for ' ... + 'imperfect spoiling is planned. To be continued...']}; +isc.values = {iscnone iscfile}; +isc.val = {iscnone}; %-------------------------------------------------------------------------- @@ -461,7 +475,7 @@ subj.tag = 'subj'; subj.name = 'Subject'; subj.help = {'Specify a subject for maps calculation.'}; -subj.val = {output sensitivity b1_type raws popup}; +subj.val = {output sensitivity b1_type raws isc popup}; % --------------------------------------------------------------------- % data Data From 57c94f2c3942277e998ffcebeecf585cddab7a84 Mon Sep 17 00:00:00 2001 From: Tobias Leutritz Date: Tue, 1 Jan 2019 22:12:52 +0100 Subject: [PATCH 07/22] implementation of Baudrexel's RF spoiling correction --- hmri_create_MTProt.m | 204 ++++++++++++++++++++++++++++++----------- tbx_scfg_hmri_create.m | 19 +++- 2 files changed, 168 insertions(+), 55 deletions(-) mode change 100644 => 100755 hmri_create_MTProt.m mode change 100644 => 100755 tbx_scfg_hmri_create.m diff --git a/hmri_create_MTProt.m b/hmri_create_MTProt.m old mode 100644 new mode 100755 index 0b4d014a..d6096d50 --- 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}); 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}); coreg_mt(Pavg{PDwidx}, PT1w_forA); @@ -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,43 @@ 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]); + 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 + 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 +674,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 +683,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 +767,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 +787,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 +880,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 +898,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); @@ -1322,46 +1380,80 @@ function PDcalculation(fA, fTPM, mpm_params) % if T1w and PDw data available, identify the protocol to define imperfect % spoiling correction parameters (for T1 map calculation) -ISC = hmri_get_defaults('imperfectSpoilCorr.enabled'); -if ~ISC +if isfield(jobsubj.isc,'iscnone') hmri_log(sprintf(['INFO: Imperfect spoiling correction is disabled.' ... '\nIf your data were acquired with one of the standard MPM ' ... '\nprotocols (customized MT-FLASH sequence) for which the correction ' ... '\ncoefficients are available, it is recommended to enable that option.']),mpm_params.defflags); -end -if mpm_params.PDwidx && mpm_params.T1widx && ISC - % retrieve all available protocols: - MPMacq_sets = hmri_get_defaults('MPMacq_set'); - % current protocol is defined by [TR_pdw TR_t1w fa_pdw fa_t1w]: - MPMacq_prot = [mpm_params.input(mpm_params.PDwidx).TR; - mpm_params.input(mpm_params.T1widx).TR; - mpm_params.input(mpm_params.PDwidx).fa; - mpm_params.input(mpm_params.T1widx).fa]'; - % then match the values and find protocol tag - nsets = numel(MPMacq_sets.vals); - ii = 0; mtch = false; - while ~mtch && ii < nsets - ii = ii+1; - if all(MPMacq_prot == MPMacq_sets.vals{ii}) - mtch = true; - prot_tag = MPMacq_sets.tags{ii}; - hmri_log(sprintf(['INFO: MPM acquisition protocol = %s.' ... - '\n\tThe coefficients corresponding to this protocol will be applied' ... - '\n\tto correct for imperfect spoiling. Please check carefully that' ... - '\n\tthe protocol used is definitely the one for which the ' ... - '\n\tcoefficients have been calculated.'],prot_tag),mpm_params.defflags); +elseif mpm_params.PDwidx && mpm_params.T1widx + if isfield(jobsubj.isc,'iscfile') + % old method here - has to be changed! (TODO) + % retrieve all available protocols: + MPMacq_sets = hmri_get_defaults('MPMacq_set'); + % current protocol is defined by [TR_pdw TR_t1w fa_pdw fa_t1w]: + MPMacq_prot = [mpm_params.input(mpm_params.PDwidx).TR; + mpm_params.input(mpm_params.T1widx).TR; + mpm_params.input(mpm_params.PDwidx).fa; + mpm_params.input(mpm_params.T1widx).fa]'; + % then match the values and find protocol tag + nsets = numel(MPMacq_sets.vals); + ii = 0; mtch = false; + while ~mtch && ii < nsets + ii = ii+1; + if all(MPMacq_prot == MPMacq_sets.vals{ii}) + mtch = true; + prot_tag = MPMacq_sets.tags{ii}; + hmri_log(sprintf(['INFO: MPM acquisition protocol = %s.' ... + '\n\tThe coefficients corresponding to this protocol will be applied' ... + '\n\tto correct for imperfect spoiling. Please check carefully that' ... + '\n\tthe protocol used is definitely the one for which the ' ... + '\n\tcoefficients have been calculated.'],prot_tag),mpm_params.defflags); + end + end + if ~mtch + prot_tag = 'Unknown'; + hmri_log(sprintf(['WARNING: MPM protocol unknown. ' ... + '\n\tCorrection for imperfect spoiling will not be applied.']),mpm_params.defflags); + end + % now retrieve imperfect spoiling correction coefficients + mpm_params.proc.ISC = hmri_get_defaults(['imperfectSpoilCorr.',prot_tag]); + elseif isfield(jobsubj.isc,'iscphase') + P = []; + % set parameters for spoiling correction due to Baudrexel et al. MRM 2018 (DOI 10.1002/mrm.26979) + if jobsubj.isc.iscphase == 50 + P = [9.639e-1, 4.989e-3, -1.254e-4, -3.180e-6, 1.527e-7, -1.462e-9; ... + 5.880e-3, -1.056e-3, 4.801e-5, -8.549e-7, 5.382e-9, 0; ... + 4.143e-4, -4.920e-6, -1.560e-7, 2.282e-9, 0, 0; ... + -1.5059e-5, 2.334e-7, -1.189e-9, 0, 0, 0; ... + 9.449e-8, -1.025e-9, 0, 0, 0, 0; ... + -4.255e-10, 0, 0, 0, 0, 0]; + elseif jobsubj.isc.iscphase == 117 + P = [9.381e-1, 4.266e-3, 2.535e-4, -2.289e-5, 5.402e-7, -4.146e-9; ... + 1.653e-2, -2.172e-3, 7.491e-5, -1.051e-6, 5.331e-9, 0; ... + 3.145e-4, 3.704e-5, -1.123e-6, 8.369e-9, 0, 0; ... + -3.848e-5, 2.773e-7, 1.662e-9, 0, 0, 0; ... + 6.230e-7, -4.019e-9, 0, 0, 0, 0; ... + -2.988e-9, 0, 0, 0, 0, 0]; + elseif jobsubj.isc.iscphase == 150 + P = [6.678e-1, 9.131e-2, -7.728e-3, 2.863e-4, -4.869e-6, 3.112e-8; ... + -3.710e-2, 2.845e-3, -7.786e-5, 8.546e-7, 2.837e-9, 0; ... + 1.448e-3, -7.537e-5, 1.403e-6, -8.865e-9, 0, 0; ... + -2.181e-5, 6.141e-7, -5.141e-9, 0, 0, 0; ... + 1.990e-7, -1.978e-9, 0, 0, 0, 0; ... + -8.617e-10, 0, 0, 0, 0, 0]; + end + if isempty(P) + hmri_log(sprintf(['WARNING: Imperfect Spoiling Correction can not be applied \n' ... + ' with the given Phase Increment of %i°, which is only available for 50°, 117°, and 150°.' ... + ' \nPlease check inputs!'],jobsubj.isc.iscphase), mpm_params.defflags); + mpm_params.proc.ISC.enabled = false; + else + mpm_params.proc.ISC = struct('iscphase',jobsubj.isc,'P',P); + hmri_log(sprintf(['INFO: Imperfect Spoiling Correction applied \n' ... + ' according to the given Phase Increment of %i°.'],jobsubj.isc.iscphase), mpm_params.nopuflags); end end - if ~mtch - prot_tag = 'Unknown'; - hmri_log(sprintf(['WARNING: MPM protocol unknown. ' ... - '\n\tCorrection for imperfect spoiling will not be applied.']),mpm_params.defflags); - end -else - prot_tag = 'Unknown'; end -% now retrieve imperfect spoiling correction coefficients -mpm_params.proc.ISC = hmri_get_defaults(['imperfectSpoilCorr.',prot_tag]); % RF sensitivity bias correction mpm_params.proc.RFsenscorr = jobsubj.sensitivity; @@ -1456,8 +1548,10 @@ function PDcalculation(fA, fTPM, mpm_params) mpm_params.output(coutput).suffix = 'R1'; mpm_params.output(coutput).units = 's-1'; mpm_params.output(coutput).descrip{1} = 'R1 map [s-1]'; - if mpm_params.proc.ISC.enabled - mpm_params.output(coutput).descrip{end+1} = '- Imperfect Spoiling Correction applied'; + if isfield(mpm_params.proc.ISC,'iscfile') + mpm_params.output(coutput).descrip{end+1} = '- Imperfect Spoiling Correction applied (a la Preibisch and Deichmann)'; + elseif isfield(mpm_params.proc.ISC,'iscphase') + mpm_params.output(coutput).descrip{end+1} = '- Imperfect Spoiling Correction applied (a la Baudrexel et al.)'; else mpm_params.output(coutput).descrip{end+1} = '- no Imperfect Spoiling Correction applied'; end @@ -1511,6 +1605,9 @@ function PDcalculation(fA, fTPM, mpm_params) end otherwise mpm_params.output(coutput).descrip{end+1} = sprintf('- B1+ bias correction using provided B1 map (%s)',B1transcorr{1}); + if isfield(mpm_params.proc.ISC,'iscphase') + mpm_params.output(coutput).descrip{end+1} = '- Imperfect Spoiling Correction applied (a la Baudrexel et al.)'; + end end switch RFsenscorr{1} case 'RF_none' @@ -1555,6 +1652,9 @@ function PDcalculation(fA, fTPM, mpm_params) end otherwise mpm_params.output(coutput).descrip{end+1} = sprintf('- B1+ bias correction using provided B1 map (%s)',B1transcorr{1}); + if isfield(mpm_params.proc.ISC,'iscphase') + mpm_params.output(coutput).descrip{end+1} = '- Imperfect Spoiling Correction applied (a la Baudrexel et al.)'; + end end switch RFsenscorr{1} case 'RF_none' diff --git a/tbx_scfg_hmri_create.m b/tbx_scfg_hmri_create.m old mode 100644 new mode 100755 index 77829f88..d2d6d3a9 --- a/tbx_scfg_hmri_create.m +++ b/tbx_scfg_hmri_create.m @@ -31,6 +31,18 @@ iscfile.ufilter = '.*'; iscfile.num = [1 1]; +%-------------------------------------------------------------------------- +% Imperfect spoiling correction - according to Phase Increment +%-------------------------------------------------------------------------- +iscphase = cfg_entry; +iscphase.tag = 'iscphase'; +iscphase.name = 'Specify phase increment'; +iscphase.help = {['Specify the phase increment of the spoiled ' ... + ' gradient echo sequence to correct imperfect spoiling. '... + ' Possible values (in degree): 50 (Siemens), 117 (General Electrics), and 150 (Philips).']}; +%scphase.values = {50 117 150}; +iscphase.val = {50}; + %-------------------------------------------------------------------------- % Imperfect spoiling correction - select coefficients %-------------------------------------------------------------------------- @@ -50,7 +62,8 @@ 'For the most standard MPM protocols using customised FLASH sequences on Siemens ' ... 'scanners, the hMRI-toolbox provides spoiling correction coefficients. ' ... 'For other protocols, the coefficients can be calculated seperately ' ... - 'and loaded here.'],[''], ... + 'and loaded here.'], ... + '- Specify RF spoiling phase increment for correction according to Baudrexel et al. 2018',[''], ... ['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. ' ... @@ -61,7 +74,7 @@ ['ISC-MODULE: an additional hMRI module that efficiently calculates the ' ... 'protocol-specific correction parameters required to account for ' ... 'imperfect spoiling is planned. To be continued...']}; -isc.values = {iscnone iscfile}; +isc.values = {iscnone iscfile iscphase}; isc.val = {iscnone}; @@ -571,4 +584,4 @@ dep = cdep; end -%_______________________________________________________________________ \ No newline at end of file +%_______________________________________________________________________ From b3b915910aabbd7d390ac529fdd9850c12fb9977 Mon Sep 17 00:00:00 2001 From: Tobias Leutritz Date: Fri, 24 Jan 2020 12:56:53 +0100 Subject: [PATCH 08/22] fix missing bias correction on PD map --- hmri_create_MTProt.m | 38 +++++++++++++++++--------------------- 1 file changed, 17 insertions(+), 21 deletions(-) diff --git a/hmri_create_MTProt.m b/hmri_create_MTProt.m index d6096d50..f0674067 100755 --- a/hmri_create_MTProt.m +++ b/hmri_create_MTProt.m @@ -1150,27 +1150,23 @@ function PDcalculation(fA, fTPM, mpm_params) % Bias-field correction of masked A map % use unified segmentation with uniform defaults across the toolbox: -if isfield(mpm_params.proc.RFsenscorr,'RF_us') - job_bfcorr = hmri_get_defaults('segment'); - job_bfcorr.channel.vols = {V_maskedA.fname}; - job_bfcorr.channel.biasreg = PDproc.biasreg; - job_bfcorr.channel.biasfwhm = PDproc.biasfwhm; - job_bfcorr.channel.write = [1 0]; % need the BiasField, obviously! - for ctis=1:length(job_bfcorr.tissue) - job_bfcorr.tissue(ctis).native = [0 0]; % no need to write c* volumes - end - output_list = spm_preproc_run(job_bfcorr); - - % Bias field correction of A map. - % Bias field calculation is based on the masked A map, while correction - % must be applied to the unmasked A map. The BiasField is therefore - % retrieved from previous step and applied onto the original A map. - BFfnam = output_list.channel.biasfield{1}; - BF = double(spm_read_vols(spm_vol(BFfnam))); - Y = BF.*spm_read_vols(spm_vol(fA)); -else - Y = spm_read_vols(spm_vol(fA)); -end +job_bfcorr = hmri_get_defaults('segment'); +job_bfcorr.channel.vols = {V_maskedA.fname}; +job_bfcorr.channel.biasreg = PDproc.biasreg; +job_bfcorr.channel.biasfwhm = PDproc.biasfwhm; +job_bfcorr.channel.write = [1 0]; % need the BiasField, obviously! +for ctis=1:length(job_bfcorr.tissue) + job_bfcorr.tissue(ctis).native = [0 0]; % no need to write c* volumes +end +output_list = spm_preproc_run(job_bfcorr); + +% Bias field correction of A map. +% Bias field calculation is based on the masked A map, while correction +% must be applied to the unmasked A map. The BiasField is therefore +% retrieved from previous step and applied onto the original A map. +BFfnam = output_list.channel.biasfield{1}; +BF = double(spm_read_vols(spm_vol(BFfnam))); +Y = BF.*spm_read_vols(spm_vol(fA)); % Calibration of flattened A map to % water content using typical white % matter value from the litterature (69%) From 41da8f772264ad3d7b53e4228b34d3c87d3df66f Mon Sep 17 00:00:00 2001 From: Tobias Leutritz Date: Sat, 11 Jul 2020 23:58:46 +0200 Subject: [PATCH 09/22] move parameters for correcting imperfect RF spoiling to external files in config/local/ISC_parameters --- config/hmri_defaults.m | 110 ------------------ .../ISC_parameters/01_classic_FIL_protocol.m | 9 ++ .../02_new_FIL_Helms_protocol.m | 9 ++ .../03_Siemens_product_sequence_Lausanne.m | 9 ++ .../04_high_res_0.8mm_FIL_protocol.m | 9 ++ .../05_NEW_High_res_0.8mm_FIL_protocol | 9 ++ .../ISC_parameters/06_NEW_1mm_protocol_v2k.m | 10 ++ .../ISC_parameters/07_800um_protocol_v3star.m | 13 +++ config/local/ISC_parameters/README.txt | 14 +++ 9 files changed, 82 insertions(+), 110 deletions(-) mode change 100644 => 100755 config/hmri_defaults.m create mode 100644 config/local/ISC_parameters/01_classic_FIL_protocol.m create mode 100644 config/local/ISC_parameters/02_new_FIL_Helms_protocol.m create mode 100644 config/local/ISC_parameters/03_Siemens_product_sequence_Lausanne.m create mode 100644 config/local/ISC_parameters/04_high_res_0.8mm_FIL_protocol.m create mode 100644 config/local/ISC_parameters/05_NEW_High_res_0.8mm_FIL_protocol create mode 100644 config/local/ISC_parameters/06_NEW_1mm_protocol_v2k.m create mode 100644 config/local/ISC_parameters/07_800um_protocol_v3star.m create mode 100644 config/local/ISC_parameters/README.txt diff --git a/config/hmri_defaults.m b/config/hmri_defaults.m old mode 100644 new mode 100755 index 4cbdb6d9..4b6d8f3b --- a/config/hmri_defaults.m +++ b/config/hmri_defaults.m @@ -258,116 +258,6 @@ 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. -% When enabling the imperfect spoiling correction, make sure the -% coefficients retrieved in the list below 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]; -% 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/01_classic_FIL_protocol.m b/config/local/ISC_parameters/01_classic_FIL_protocol.m new file mode 100644 index 00000000..e6a0b2cd --- /dev/null +++ b/config/local/ISC_parameters/01_classic_FIL_protocol.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/02_new_FIL_Helms_protocol.m b/config/local/ISC_parameters/02_new_FIL_Helms_protocol.m new file mode 100644 index 00000000..4965176d --- /dev/null +++ b/config/local/ISC_parameters/02_new_FIL_Helms_protocol.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/03_Siemens_product_sequence_Lausanne.m b/config/local/ISC_parameters/03_Siemens_product_sequence_Lausanne.m new file mode 100644 index 00000000..b8d7f1fb --- /dev/null +++ b/config/local/ISC_parameters/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/04_high_res_0.8mm_FIL_protocol.m b/config/local/ISC_parameters/04_high_res_0.8mm_FIL_protocol.m new file mode 100644 index 00000000..ebb45aef --- /dev/null +++ b/config/local/ISC_parameters/04_high_res_0.8mm_FIL_protocol.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/05_NEW_High_res_0.8mm_FIL_protocol b/config/local/ISC_parameters/05_NEW_High_res_0.8mm_FIL_protocol new file mode 100644 index 00000000..cbf9a936 --- /dev/null +++ b/config/local/ISC_parameters/05_NEW_High_res_0.8mm_FIL_protocol @@ -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/06_NEW_1mm_protocol_v2k.m b/config/local/ISC_parameters/06_NEW_1mm_protocol_v2k.m new file mode 100644 index 00000000..85629883 --- /dev/null +++ b/config/local/ISC_parameters/06_NEW_1mm_protocol_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/07_800um_protocol_v3star.m b/config/local/ISC_parameters/07_800um_protocol_v3star.m new file mode 100644 index 00000000..73a604b5 --- /dev/null +++ b/config/local/ISC_parameters/07_800um_protocol_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/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. From aa5839e5a9c7e754bf0dda39c754dbed3b255216 Mon Sep 17 00:00:00 2001 From: Tobias Leutritz Date: Sun, 12 Jul 2020 00:00:29 +0200 Subject: [PATCH 10/22] create external file for imperfect RF spoiling correction parameters --- EPG_for_hmri_tbx/Unit_Test_Protocol.m | 8 ++++++++ hmri_corr_imperf_spoil.m | 18 +++++++++++++----- 2 files changed, 21 insertions(+), 5 deletions(-) create mode 100755 EPG_for_hmri_tbx/Unit_Test_Protocol.m mode change 100644 => 100755 hmri_corr_imperf_spoil.m 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/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 From 94cdbe80621dfee6f7a64dc5e3812e2eaab95a9e Mon Sep 17 00:00:00 2001 From: Tobias Leutritz Date: Sun, 12 Jul 2020 00:01:36 +0200 Subject: [PATCH 11/22] provide dependency for loading imperfect RF spoiling correction parameters in batch system --- tbx_scfg_hmri_create.m | 4 ++-- tbx_scfg_hmri_imperf_spoil.m | 13 ++++++++++++- 2 files changed, 14 insertions(+), 3 deletions(-) mode change 100644 => 100755 tbx_scfg_hmri_imperf_spoil.m diff --git a/tbx_scfg_hmri_create.m b/tbx_scfg_hmri_create.m index d2d6d3a9..dbddf8c3 100755 --- a/tbx_scfg_hmri_create.m +++ b/tbx_scfg_hmri_create.m @@ -24,9 +24,9 @@ iscfile = cfg_files; iscfile.tag = 'iscfile'; iscfile.name = 'Select ISC file'; -iscfile.help = {['Select a *.json file containing the correction coefficients ' ... +iscfile.help = {['Select an *.m file containing the correction coefficients ' ... 'corresponding to the current acquisition protocol.']}; -iscfile.filter = 'json'; +iscfile.filter = 'm'; iscfile.dir = fullfile(fileparts(mfilename('fullpath')),'config'); iscfile.ufilter = '.*'; iscfile.num = [1 1]; diff --git a/tbx_scfg_hmri_imperf_spoil.m b/tbx_scfg_hmri_imperf_spoil.m old mode 100644 new mode 100755 index 49d25924..c3bd75a4 --- 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,14 @@ '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; +imperf_spoil.vout = @vout_imperf_spoil; +end + +% output m-file filename as dependency +function dep = vout_imperf_spoil(job) +out.ISC_file = fullfile(job.outdir,[strrep(job.prot_name,' ',''),'.m']); %#ok +dep(1) = cfg_dep; +dep(1).sname = 'ISC file'; +dep(1).src_output = substruct('.','ISC_file');spm fmri +dep(1).tgt_spec = cfg_findspec({{'filter','file'}}); +end From d81b15873d3d4ade73b62a067b3a2217fb31aec3 Mon Sep 17 00:00:00 2001 From: Tobias Leutritz Date: Sun, 12 Jul 2020 00:03:05 +0200 Subject: [PATCH 12/22] load parameters for correction of imperfect spoiling correction from external m file --- hmri_create_MTProt.m | 82 +++++++++++++++++++++++++------------------- 1 file changed, 47 insertions(+), 35 deletions(-) diff --git a/hmri_create_MTProt.m b/hmri_create_MTProt.m index a1c4301c..97c6f1f8 100755 --- a/hmri_create_MTProt.m +++ b/hmri_create_MTProt.m @@ -674,7 +674,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 + ... @@ -1376,47 +1376,59 @@ function PDcalculation(fA, fTPM, mpm_params) % if T1w and PDw data available, identify the protocol to define imperfect % spoiling correction parameters (for T1 map calculation) -ISC = hmri_get_defaults('imperfectSpoilCorr.enabled'); -if ~ISC +if isfield(jobsubj.isc,'iscnone') hmri_log(sprintf(['INFO: Imperfect spoiling correction is disabled.' ... '\nIf your data were acquired with one of the standard MPM ' ... '\nprotocols (customized MT-FLASH sequence) for which the correction ' ... '\ncoefficients are available, it is recommended to enable that option.']),mpm_params.defflags); -end -if mpm_params.PDwidx && mpm_params.T1widx && ISC - % retrieve all available protocols: - MPMacq_sets = hmri_get_defaults('MPMacq_set'); - % current protocol is defined by [TR_pdw TR_t1w fa_pdw fa_t1w]: - MPMacq_prot = [mpm_params.input(mpm_params.PDwidx).TR; - mpm_params.input(mpm_params.T1widx).TR; - mpm_params.input(mpm_params.PDwidx).fa; - mpm_params.input(mpm_params.T1widx).fa]'; - % then match the values and find protocol tag - nsets = numel(MPMacq_sets.vals); - ii = 0; mtch = false; - while ~mtch && ii < nsets - ii = ii+1; - if all(MPMacq_prot == MPMacq_sets.vals{ii}) - mtch = true; - prot_tag = MPMacq_sets.tags{ii}; - hmri_log(sprintf(['INFO: MPM acquisition protocol = %s.' ... - '\n\tThe coefficients corresponding to this protocol will be applied' ... - '\n\tto correct for imperfect spoiling. Please check carefully that' ... - '\n\tthe protocol used is definitely the one for which the ' ... - '\n\tcoefficients have been calculated.'],prot_tag),mpm_params.defflags); +elseif mpm_params.PDwidx && mpm_params.T1widx + if isfield(jobsubj.isc,'iscfile') + eval(sprintf('run(''%s'')',jobsubj.isc.iscfile{1})); + MPMacq_set = hmri_get_defaults('MPMacq_set'); + prot_tag = MPMacq_set.tag; + hmri_log(sprintf(['INFO: MPM acquisition protocol = %s.' ... + '\n\tThe coefficients corresponding to this protocol will be applied' ... + '\n\tto correct for imperfect spoiling.'],prot_tag),mpm_params.defflags); + % now retrieve imperfect spoiling correction coefficients + mpm_params.proc.ISC = hmri_get_defaults(sprintf('imperfectSpoilCorr.%s',prot_tag)); + mpm_params.proc.ISC.iscfile = 'iscfile'; + elseif isfield(jobsubj.isc,'iscphase') + P = []; + % set parameters for spoiling correction due to Baudrexel et al. MRM 2018 (DOI 10.1002/mrm.26979) + if jobsubj.isc.iscphase == 50 + P = [9.639e-1, 4.989e-3, -1.254e-4, -3.180e-6, 1.527e-7, -1.462e-9; ... + 5.880e-3, -1.056e-3, 4.801e-5, -8.549e-7, 5.382e-9, 0; ... + 4.143e-4, -4.920e-6, -1.560e-7, 2.282e-9, 0, 0; ... + -1.5059e-5, 2.334e-7, -1.189e-9, 0, 0, 0; ... + 9.449e-8, -1.025e-9, 0, 0, 0, 0; ... + -4.255e-10, 0, 0, 0, 0, 0]; + elseif jobsubj.isc.iscphase == 117 + P = [9.381e-1, 4.266e-3, 2.535e-4, -2.289e-5, 5.402e-7, -4.146e-9; ... + 1.653e-2, -2.172e-3, 7.491e-5, -1.051e-6, 5.331e-9, 0; ... + 3.145e-4, 3.704e-5, -1.123e-6, 8.369e-9, 0, 0; ... + -3.848e-5, 2.773e-7, 1.662e-9, 0, 0, 0; ... + 6.230e-7, -4.019e-9, 0, 0, 0, 0; ... + -2.988e-9, 0, 0, 0, 0, 0]; + elseif jobsubj.isc.iscphase == 150 + P = [6.678e-1, 9.131e-2, -7.728e-3, 2.863e-4, -4.869e-6, 3.112e-8; ... + -3.710e-2, 2.845e-3, -7.786e-5, 8.546e-7, 2.837e-9, 0; ... + 1.448e-3, -7.537e-5, 1.403e-6, -8.865e-9, 0, 0; ... + -2.181e-5, 6.141e-7, -5.141e-9, 0, 0, 0; ... + 1.990e-7, -1.978e-9, 0, 0, 0, 0; ... + -8.617e-10, 0, 0, 0, 0, 0]; + end + if isempty(P) + hmri_log(sprintf(['WARNING: Imperfect Spoiling Correction can not be applied \n' ... + ' with the given Phase Increment of %i°, which is only available for 50°, 117°, and 150°.' ... + ' \nPlease check inputs!'],jobsubj.isc.iscphase), mpm_params.defflags); + mpm_params.proc.ISC.enabled = false; + else + mpm_params.proc.ISC = struct('iscphase',jobsubj.isc,'P',P); + hmri_log(sprintf(['INFO: Imperfect Spoiling Correction applied \n' ... + ' according to the given Phase Increment of %i°.'],jobsubj.isc.iscphase), mpm_params.nopuflags); end end - if ~mtch - prot_tag = 'Unknown'; - hmri_log(sprintf(['WARNING: MPM protocol unknown. ' ... - '\n\tCorrection for imperfect spoiling will not be applied.']),mpm_params.defflags); - end -else - prot_tag = 'Unknown'; end -% now retrieve imperfect spoiling correction coefficients -mpm_params.proc.ISC = hmri_get_defaults(['imperfectSpoilCorr.',prot_tag]); -mpm_params.proc.ISC.enabled = ISC; % RF sensitivity bias correction mpm_params.proc.RFsenscorr = jobsubj.sensitivity; From ee12423b2fafd1d17ab1f873ec8d70fb0c06cc63 Mon Sep 17 00:00:00 2001 From: Tobias Leutritz Date: Sun, 12 Jul 2020 00:18:12 +0200 Subject: [PATCH 13/22] add SDAM B1 mapping protocol --- hmri_create_b1map.m | 93 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 93 insertions(+) diff --git a/hmri_create_b1map.m b/hmri_create_b1map.m index 1ce2ab0d..bd202a05 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 @@ -187,6 +191,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 %=========================================================================% @@ -788,6 +861,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); From ae7967bce50512607d5bc4c2237032c2cb497cf1 Mon Sep 17 00:00:00 2001 From: Tobias Leutritz Date: Sun, 12 Jul 2020 00:19:51 +0200 Subject: [PATCH 14/22] correct data scaling for B1 maps --- hmri_create_b1map.m | 14 +++----------- 1 file changed, 3 insertions(+), 11 deletions(-) diff --git a/hmri_create_b1map.m b/hmri_create_b1map.m index bd202a05..c657e776 100755 --- a/hmri_create_b1map.m +++ b/hmri_create_b1map.m @@ -158,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); @@ -536,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); @@ -595,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); From 0354b71c57dc34f2e2f1f8c8b48f6966533cc1db Mon Sep 17 00:00:00 2001 From: Tobias Leutritz Date: Sun, 12 Jul 2020 00:22:44 +0200 Subject: [PATCH 15/22] add SDAM B1 mapping protocol --- config/hmri_b1_standard_defaults.m | 8 +++++++- tbx_scfg_hmri_create.m | 17 ++++++++++++++++- 2 files changed, 23 insertions(+), 2 deletions(-) mode change 100644 => 100755 config/hmri_b1_standard_defaults.m diff --git a/config/hmri_b1_standard_defaults.m b/config/hmri_b1_standard_defaults.m old mode 100644 new mode 100755 index 3e66d692..46b5199d --- 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/tbx_scfg_hmri_create.m b/tbx_scfg_hmri_create.m index dbddf8c3..4589da4b 100755 --- a/tbx_scfg_hmri_create.m +++ b/tbx_scfg_hmri_create.m @@ -293,6 +293,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 % --------------------------------------------------------------------- @@ -342,6 +356,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.'] @@ -352,7 +367,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}; % --------------------------------------------------------------------- From cf93f6fcb6025358c1b2cfcc3cfc8809af64d65b Mon Sep 17 00:00:00 2001 From: Tobias Leutritz Date: Mon, 13 Jul 2020 00:05:19 +0200 Subject: [PATCH 16/22] add basic quality tool by EB (visual check) --- tbx_cfg_hmri.m | 2 +- tbx_scfg_hmri_quality.m | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) mode change 100644 => 100755 tbx_cfg_hmri.m mode change 100644 => 100755 tbx_scfg_hmri_quality.m 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_quality.m b/tbx_scfg_hmri_quality.m old mode 100644 new mode 100755 index 14c6a6b6..be39082e --- a/tbx_scfg_hmri_quality.m +++ b/tbx_scfg_hmri_quality.m @@ -235,6 +235,7 @@ function hmri_quality_visual_check(job) 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'); From 83643b4cedb6a2181910b8acbea83a957ad125ed Mon Sep 17 00:00:00 2001 From: Tobias Leutritz Date: Mon, 13 Jul 2020 00:08:10 +0200 Subject: [PATCH 17/22] change help texts according to current implementation --- tbx_scfg_hmri_create.m | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/tbx_scfg_hmri_create.m b/tbx_scfg_hmri_create.m index 4589da4b..edb5ff6a 100755 --- a/tbx_scfg_hmri_create.m +++ b/tbx_scfg_hmri_create.m @@ -60,9 +60,12 @@ '- None [default]: no correction applied. ', ... ['- Select ISC file: you must choose a set of correction coefficients. ' ... 'For the most standard MPM protocols using customised FLASH sequences on Siemens ' ... - 'scanners, the hMRI-toolbox provides spoiling correction coefficients. ' ... + 'scanners, the hMRI-toolbox provides spoiling correction coefficients in '... + 'the folder config/local/ISC_parameters. ' ... 'For other protocols, the coefficients can be calculated seperately ' ... - 'and loaded here.'], ... + 'by the additional hMRI module ''Imperfect Spoiling Calc.'' that efficiently ' ... + 'calculates the protocol-specific correction parameters required to account for ' ... + 'imperfect spoiling, which can be loaded here (also via dependency).'], ... '- Specify RF spoiling phase increment for correction according to Baudrexel et al. 2018',[''], ... ['WARNING: although the TR ' ... 'and FA values are key parameters in simulating the imperfect spoiling, ' ... @@ -70,10 +73,7 @@ '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.'],[''], ... - ['ISC-MODULE: an additional hMRI module that efficiently calculates the ' ... - 'protocol-specific correction parameters required to account for ' ... - 'imperfect spoiling is planned. To be continued...']}; + 'use identical correction coefficients.'],['']}; isc.values = {iscnone iscfile iscphase}; isc.val = {iscnone}; From dca4a8e3550149fc9e0bd3f5cce25419880d002a Mon Sep 17 00:00:00 2001 From: Tobias Leutritz Date: Mon, 13 Jul 2020 00:10:33 +0200 Subject: [PATCH 18/22] fix typo --- tbx_scfg_hmri_imperf_spoil.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tbx_scfg_hmri_imperf_spoil.m b/tbx_scfg_hmri_imperf_spoil.m index c3bd75a4..a9280d1e 100755 --- a/tbx_scfg_hmri_imperf_spoil.m +++ b/tbx_scfg_hmri_imperf_spoil.m @@ -182,6 +182,6 @@ out.ISC_file = fullfile(job.outdir,[strrep(job.prot_name,' ',''),'.m']); %#ok dep(1) = cfg_dep; dep(1).sname = 'ISC file'; -dep(1).src_output = substruct('.','ISC_file');spm fmri +dep(1).src_output = substruct('.','ISC_file'); dep(1).tgt_spec = cfg_findspec({{'filter','file'}}); end From c28fd1dade12670ca8af057555ed2ebb0afeb177 Mon Sep 17 00:00:00 2001 From: Tobias Leutritz Date: Wed, 9 Sep 2020 21:40:01 +0200 Subject: [PATCH 19/22] fix unset ISC structure when not enabled at all --- hmri_create_MTProt.m | 1 + 1 file changed, 1 insertion(+) diff --git a/hmri_create_MTProt.m b/hmri_create_MTProt.m index 97c6f1f8..356a551f 100755 --- a/hmri_create_MTProt.m +++ b/hmri_create_MTProt.m @@ -1381,6 +1381,7 @@ function PDcalculation(fA, fTPM, mpm_params) '\nIf your data were acquired with one of the standard MPM ' ... '\nprotocols (customized MT-FLASH sequence) for which the correction ' ... '\ncoefficients are available, it is recommended to enable that option.']),mpm_params.defflags); + mpm_params.proc.ISC.enabled = false; elseif mpm_params.PDwidx && mpm_params.T1widx if isfield(jobsubj.isc,'iscfile') eval(sprintf('run(''%s'')',jobsubj.isc.iscfile{1})); From 3098800a33b469c7c1db022a50da6a586d967506 Mon Sep 17 00:00:00 2001 From: Tobias Leutritz Date: Sat, 21 Nov 2020 22:30:22 +0100 Subject: [PATCH 20/22] reinsert coreg_flags after faulty merge --- hmri_create_MTProt.m | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/hmri_create_MTProt.m b/hmri_create_MTProt.m index 10ba6417..2e3d4a44 100755 --- a/hmri_create_MTProt.m +++ b/hmri_create_MTProt.m @@ -1507,6 +1507,16 @@ function PDcalculation(fA, fTPM, mpm_params) % coregistration of all images to the PDw average (or TE=0 fit): mpm_params.coreg = hmri_get_defaults('coreg2PDw'); +% coregistration flags for weighted images +mpm_params.coreg_flags = hmri_get_defaults('coreg_flags'); +hmri_log(sprintf('=== Registration Settings for weighted images ==='),mpm_params.nopuflags); +hmri_log(sprintf('Method: %s, Sampling: %s, Smoothing: %s', mpm_params.coreg_flags.cost_fun, mat2str(mpm_params.coreg_flags.sep), mat2str(mpm_params.coreg_flags.fwhm)),mpm_params.nopuflags); + +% coregistration flags for B1 to PDw +mpm_params.coreg_bias_flags = hmri_get_defaults('coreg_bias_flags'); +hmri_log(sprintf('=== Registration Settings for B1 bias images ==='),mpm_params.nopuflags); +hmri_log(sprintf('Method: %s, Sampling: %s, Smoothing: %s', mpm_params.coreg_bias_flags.cost_fun, mat2str(mpm_params.coreg_bias_flags.sep), mat2str(mpm_params.coreg_bias_flags.fwhm)),mpm_params.nopuflags); + % Prepare output for R1, PD, MT and R2* maps RFsenscorr = fieldnames(mpm_params.proc.RFsenscorr); B1transcorr = fieldnames(jobsubj.b1_type); From ca43ec071d93b025cfe1b24daf20d0a44c2441d6 Mon Sep 17 00:00:00 2001 From: Tobias Leutritz Date: Fri, 27 Nov 2020 21:27:49 +0100 Subject: [PATCH 21/22] rebuild imperfect spoiling correction to get rid of extra GUI entry and fix few bugs in that context --- config/hmri_b1_standard_defaults.m | 0 config/hmri_defaults.m | 46 ++- ...L_protocol.m => protocol_01_classic_FIL.m} | 0 ...protocol.m => protocol_02_new_FIL_Helms.m} | 0 ...ol_03_Siemens_product_sequence_Lausanne.m} | 0 ...col.m => protocol_04_high_res_0.8mm_FIL.m} | 0 ...l => protocol_05_NEW_High_res_0.8mm_FIL.m} | 0 ...otocol_v2k.m => protocol_06_NEW_1mm_v2k.m} | 0 ...ol_v3star.m => protocol_07_800um_v3star.m} | 0 config/local/hmri_local_defaults.m | 168 +++----- config/local/hmri_local_defaults_iscfile.m | 366 ++++++++++++++++++ config/local/hmri_local_defaults_iscphase.m | 366 ++++++++++++++++++ hmri_create_MTProt.m | 29 +- tbx_scfg_hmri_create.m | 73 +--- tbx_scfg_hmri_imperf_spoil.m | 10 - 15 files changed, 846 insertions(+), 212 deletions(-) mode change 100755 => 100644 config/hmri_b1_standard_defaults.m mode change 100755 => 100644 config/hmri_defaults.m rename config/local/ISC_parameters/{01_classic_FIL_protocol.m => protocol_01_classic_FIL.m} (100%) rename config/local/ISC_parameters/{02_new_FIL_Helms_protocol.m => protocol_02_new_FIL_Helms.m} (100%) rename config/local/ISC_parameters/{03_Siemens_product_sequence_Lausanne.m => protocol_03_Siemens_product_sequence_Lausanne.m} (100%) rename config/local/ISC_parameters/{04_high_res_0.8mm_FIL_protocol.m => protocol_04_high_res_0.8mm_FIL.m} (100%) rename config/local/ISC_parameters/{05_NEW_High_res_0.8mm_FIL_protocol => protocol_05_NEW_High_res_0.8mm_FIL.m} (100%) rename config/local/ISC_parameters/{06_NEW_1mm_protocol_v2k.m => protocol_06_NEW_1mm_v2k.m} (100%) rename config/local/ISC_parameters/{07_800um_protocol_v3star.m => protocol_07_800um_v3star.m} (100%) create mode 100644 config/local/hmri_local_defaults_iscfile.m create mode 100644 config/local/hmri_local_defaults_iscphase.m mode change 100755 => 100644 hmri_create_MTProt.m mode change 100755 => 100644 tbx_scfg_hmri_create.m mode change 100755 => 100644 tbx_scfg_hmri_imperf_spoil.m diff --git a/config/hmri_b1_standard_defaults.m b/config/hmri_b1_standard_defaults.m old mode 100755 new mode 100644 diff --git a/config/hmri_defaults.m b/config/hmri_defaults.m old mode 100755 new mode 100644 index b99e2d3c..d1147189 --- 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,6 +280,50 @@ 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. +% 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; + %-------------------------------------------------------------------------- % B1 mapping processing parameters %-------------------------------------------------------------------------- diff --git a/config/local/ISC_parameters/01_classic_FIL_protocol.m b/config/local/ISC_parameters/protocol_01_classic_FIL.m similarity index 100% rename from config/local/ISC_parameters/01_classic_FIL_protocol.m rename to config/local/ISC_parameters/protocol_01_classic_FIL.m diff --git a/config/local/ISC_parameters/02_new_FIL_Helms_protocol.m b/config/local/ISC_parameters/protocol_02_new_FIL_Helms.m similarity index 100% rename from config/local/ISC_parameters/02_new_FIL_Helms_protocol.m rename to config/local/ISC_parameters/protocol_02_new_FIL_Helms.m diff --git a/config/local/ISC_parameters/03_Siemens_product_sequence_Lausanne.m b/config/local/ISC_parameters/protocol_03_Siemens_product_sequence_Lausanne.m similarity index 100% rename from config/local/ISC_parameters/03_Siemens_product_sequence_Lausanne.m rename to config/local/ISC_parameters/protocol_03_Siemens_product_sequence_Lausanne.m diff --git a/config/local/ISC_parameters/04_high_res_0.8mm_FIL_protocol.m b/config/local/ISC_parameters/protocol_04_high_res_0.8mm_FIL.m similarity index 100% rename from config/local/ISC_parameters/04_high_res_0.8mm_FIL_protocol.m rename to config/local/ISC_parameters/protocol_04_high_res_0.8mm_FIL.m diff --git a/config/local/ISC_parameters/05_NEW_High_res_0.8mm_FIL_protocol b/config/local/ISC_parameters/protocol_05_NEW_High_res_0.8mm_FIL.m similarity index 100% rename from config/local/ISC_parameters/05_NEW_High_res_0.8mm_FIL_protocol rename to config/local/ISC_parameters/protocol_05_NEW_High_res_0.8mm_FIL.m diff --git a/config/local/ISC_parameters/06_NEW_1mm_protocol_v2k.m b/config/local/ISC_parameters/protocol_06_NEW_1mm_v2k.m similarity index 100% rename from config/local/ISC_parameters/06_NEW_1mm_protocol_v2k.m rename to config/local/ISC_parameters/protocol_06_NEW_1mm_v2k.m diff --git a/config/local/ISC_parameters/07_800um_protocol_v3star.m b/config/local/ISC_parameters/protocol_07_800um_v3star.m similarity index 100% rename from config/local/ISC_parameters/07_800um_protocol_v3star.m rename to config/local/ISC_parameters/protocol_07_800um_v3star.m 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_create_MTProt.m b/hmri_create_MTProt.m old mode 100755 new mode 100644 index 2e3d4a44..35da0a06 --- a/hmri_create_MTProt.m +++ b/hmri_create_MTProt.m @@ -610,7 +610,11 @@ spm_progress_bar('Init',dm(3),'Calculating FA maps (ISC)','planes completed'); for p = 1:dm(3) M = M0*spm_matrix([0 0 p]); - 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 + 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))) @@ -1370,17 +1374,16 @@ function PDcalculation(fA, fTPM, mpm_params) mpm_params.nr_echoes4avg = min(length(find(mpm_params.input(1).TE -dep(1) = cfg_dep; -dep(1).sname = 'ISC file'; -dep(1).src_output = substruct('.','ISC_file'); -dep(1).tgt_spec = cfg_findspec({{'filter','file'}}); end From 2d8fcfbfafd6ca6df028358fd256766a9caf1426 Mon Sep 17 00:00:00 2001 From: Tobias Leutritz Date: Sun, 2 May 2021 21:04:20 +0200 Subject: [PATCH 22/22] used sub-folders for contrasts to avoid mixup in case of same names of multi-echo data per contrast --- hmri_create_RFsens.m | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) mode change 100644 => 100755 hmri_create_RFsens.m diff --git a/hmri_create_RFsens.m b/hmri_create_RFsens.m old mode 100644 new mode 100755 index 0d50dc35..c8c8b7cb --- a/hmri_create_RFsens.m +++ b/hmri_create_RFsens.m @@ -58,7 +58,8 @@ %========================================================================== % Calculate quantitative RF sensitivity maps (HC/BC division) %========================================================================== - qsensmap = spm_imcalc(smoothedmaps, fullfile(calcpath, sprintf('sensMap_HC_over_BC_division_%s.nii',rfsens_params.input(ccon).tag)), 'i1./i2'); + qsensmap = spm_imcalc(smoothedmaps, fullfile(calcpath, rfsens_params.input(ccon).tag, ... + sprintf('sensMap_HC_over_BC_division_%s.nii',rfsens_params.input(ccon).tag)), 'i1./i2'); qsensmap = qsensmap.fname; % set metadata input_files = sensmaps; @@ -73,7 +74,8 @@ corrected_structurals = cell(nSTRUCT,1); for i=1:nSTRUCT - corrected_structurals{i} = fullfile(calcpath, spm_file(spm_file(structurals(i,:),'filename'),'suffix','_RFSC')); + corrected_structurals{i} = fullfile(calcpath, rfsens_params.input(ccon).tag, ... + spm_file(spm_file(structurals(i,:),'filename'),'suffix','_RFSC')); spm_imcalc({structurals(i,:), qsensmap}, corrected_structurals{i}, 'i1./i2'); % set metadata (relates only to original inputs to keep it % readable and trackable since intermediate calculation directories @@ -137,6 +139,7 @@ ccon = ccon+1; rfsens_params.input(ccon).fnam = tmpfnam; rfsens_params.input(ccon).tag = 'MT'; + mkdir(fullfile(rfsens_params.calcpath,'MT')) rfsens_params.MTidx = ccon; end % 2) try PDw contrast: @@ -147,6 +150,7 @@ ccon = ccon+1; rfsens_params.input(ccon).fnam = tmpfnam; rfsens_params.input(ccon).tag = 'PD'; + mkdir(fullfile(rfsens_params.calcpath,'PD')) rfsens_params.PDidx = ccon; end % 3) try T1w contrast: @@ -157,6 +161,7 @@ ccon = ccon+1; rfsens_params.input(ccon).fnam = tmpfnam; rfsens_params.input(ccon).tag = 'T1'; + mkdir(fullfile(rfsens_params.calcpath,'T1')) rfsens_params.T1idx = ccon; % zero index means no contrast available end rfsens_params.ncon = ccon; % number of contrasts available @@ -193,7 +198,8 @@ end for csens=1:length(raw_sens_input) tmprawfnam = spm_file(raw_sens_input{csens},'number',''); - tmpfnam{csens} = fullfile(rfsens_params.calcpath,spm_file(tmprawfnam,'filename')); + tmpfnam{csens} = fullfile(rfsens_params.calcpath, ... + rfsens_params.input(ccon).tag, spm_file(tmprawfnam,'filename')); copyfile(tmprawfnam, tmpfnam{csens}); try copyfile([spm_str_manip(tmprawfnam,'r') '.json'],[spm_str_manip(tmpfnam{csens},'r') '.json']); end end @@ -256,7 +262,7 @@ matlabbatch{1}.spm.spatial.coreg.estwrite.roptions.wrap = [0 0 0]; matlabbatch{1}.spm.spatial.coreg.estwrite.roptions.mask = 0; output_list = spm_jobman('run',matlabbatch); - coregsensfnam{i} = fullfile(calcpath, spm_file(output_list{1}.rfiles{1},'filename')); %#ok + coregsensfnam{i} = fullfile(calcpath, tag, spm_file(output_list{1}.rfiles{1},'filename')); %#ok end coregsensfnam = coregsensfnam'; % must be a nx1 cellstr clear matlabbatch