From 0010d28489b644e0724d2577761b5061715ab025 Mon Sep 17 00:00:00 2001 From: ftadel Date: Wed, 21 Apr 2021 11:07:20 +0200 Subject: [PATCH] Anatomy: SPM Segment: Enable T1+T2 segmentation --- doc/license.html | 2 +- doc/version.txt | 2 +- toolbox/anatomy/bst_normalize_mni.m | 128 ++++++++++-------- toolbox/anatomy/mri_normalize_segment.m | 44 ++++-- .../process/functions/process_mni_normalize.m | 89 ++++++++++-- 5 files changed, 184 insertions(+), 81 deletions(-) diff --git a/doc/license.html b/doc/license.html index 1b79be996..93933ec8a 100644 --- a/doc/license.html +++ b/doc/license.html @@ -5,7 +5,7 @@

THERE IS NO UNDO BUTTON - SET UP A BACKUP OF YOUR DATABASE

-Version: 3.210420 (20-Apr-2021)
+Version: 3.210421 (21-Apr-2021)
COPYRIGHT © 2000-2020 USC & McGill University.
 diff --git a/doc/version.txt b/doc/version.txt index c1450c8f5..ee17d3985 100644 --- a/doc/version.txt +++ b/doc/version.txt @@ -1,2 +1,2 @@ % Brainstorm -% v. 3.210420 (20-Apr-2021) \ No newline at end of file +% v. 3.210421 (21-Apr-2021) \ No newline at end of file diff --git a/toolbox/anatomy/bst_normalize_mni.m b/toolbox/anatomy/bst_normalize_mni.m index 6f6633048..bec845939 100644 --- a/toolbox/anatomy/bst_normalize_mni.m +++ b/toolbox/anatomy/bst_normalize_mni.m @@ -1,16 +1,18 @@ -function [sMri, errMsg] = bst_normalize_mni(MriFile, Method) +function [sMriT1, errMsg] = bst_normalize_mni(T1File, Method, T2File) % BST_NORMALIZE_MNI: Compute deformation fields to the MNI ICBM152 space. % -% USAGE: [sMri, errMsg] = bst_normalize_mni(MriFile, Method='maff8') -% [sMri, errMsg] = bst_normalize_mni(sMri, Method='maff8') -% bst_normalize_mni('install') % Only installs default SPM tpm.nii +% USAGE: [sMriT1, errMsg] = bst_normalize_mni(T1File, Method='maff8', T2File=[]) +% [sMriT1, errMsg] = bst_normalize_mni(sMriT1, Method='maff8', sMriT2=[]) +% bst_normalize_mni('install') % Only installs default SPM tpm.nii % % INPUTS: -% - MriFile : Relative path to a MRI file in the Brainstorm database -% - sMri : Brainstorm MRI structure -% - Method : String defining the method to use for the registration -% 'maff8' : SPM mutual information algorithm (affine transform) -% 'segment' : SPM12 segment +% - T1File : Relative path to a T1 MRI file in the Brainstorm database +% - sMriT1 : Brainstorm T1 MRI structure +% - Method : String defining the method to use for the registration +% 'maff8' : SPM mutual information algorithm (affine transform) +% 'segment' : SPM12 segment +% - T2File : Relative path to a T2 MRI file in the Brainstorm database +% - sMriT2 : Brainstorm T2 MRI structure % @============================================================================= % This function is part of the Brainstorm software: @@ -30,33 +32,45 @@ % For more information type "brainstorm license" at command prompt. % =============================================================================@ % -% Authors: Francois Tadel, 2015-2020 +% Authors: Francois Tadel, 2015-2021 %% ===== PARSE INPUTS ===== % Inializations global GlobalData; errMsg = []; % Usage: bst_normalize_mni('install') -if isequal(MriFile, 'install') +if isequal(T1File, 'install') isInstall = 1; - sMri = []; + sMriT1 = []; else isInstall = 0; - % Usage: bst_normalize_mni(sMri) - if ~ischar(MriFile) - sMri = MriFile; - MriFile = []; - % Usage: bst_normalize_mni(MriFile) + % Usage: bst_normalize_mni(sMriT1) + if ~ischar(T1File) + sMriT1 = T1File; + T1File = []; + % Usage: bst_normalize_mni(T1File) else - sMri = []; + sMriT1 = []; end end -% Other parameters +% Default method if (nargin < 2) || isempty(Method) Method = 'maff8'; end - - +% T2 MRI +if (nargin < 3) || isempty(T2File) + T2File = []; + sMriT2 = []; +else + if ~ischar(T2File) + sMriT2 = T2File; + T2File = []; + else + sMriT2 = []; + end +end + + %% ===== GET SPM TEMPLATE ===== % Open progress bar isProgress = bst_progress('isVisible'); @@ -113,18 +127,24 @@ %% ===== LOAD ANATOMY ===== -if isempty(sMri) +% T1 MRI +if isempty(sMriT1) % Progress bar bst_progress('text', 'Loading input MRI...'); % Check if it is loaded in memory - [sMri, iLoadedMri] = bst_memory('GetMri', MriFile); + [sMriT1, iLoadedMri] = bst_memory('GetMri', T1File); % If not: load it from the file - if isempty(sMri) - sMri = in_mri_bst(MriFile); + if isempty(sMriT1) + sMriT1 = in_mri_bst(T1File); end else iLoadedMri = []; end +% T2 MRI +if isempty(sMriT2) && ~isempty(T2File) + sMriT2 = in_mri_bst(T2File); +end + %% ===== MNI NORMALIZATION ===== TpmFiles = []; @@ -135,10 +155,10 @@ % Progress bar bst_progress('text', 'Resampling MRI...'); % Resample volume if needed - if any(abs(sMri.Voxsize - [1 1 1]) > 0.001) - [sMriRes, Tres] = mri_resample(sMri, [256 256 256], [1 1 1], 'linear'); + if any(abs(sMriT1.Voxsize - [1 1 1]) > 0.001) + [sMriRes, Tres] = mri_resample(sMriT1, [256 256 256], [1 1 1], 'linear'); else - sMriRes = sMri; + sMriRes = sMriT1; Tres = []; end % Compute affine transformation to MNI space @@ -148,8 +168,8 @@ Tmni = Tmni * Tres; end % Save results into the MRI structure - sMri.NCS.R = Tmni(1:3,1:3); - sMri.NCS.T = Tmni(1:3,4); + sMriT1.NCS.R = Tmni(1:3,1:3); + sMriT1.NCS.T = Tmni(1:3,4); % SPM12 SEGMENT case 'segment' @@ -161,15 +181,15 @@ % Progress bar bst_progress('text', 'Running SPM batch... (see command window)'); % Compute non-linear registration to MNI space - [sMri, TpmFiles] = mri_normalize_segment(sMri, TpmFile); - if isempty(sMri) + [sMriT1, TpmFiles] = mri_normalize_segment(sMriT1, TpmFile, sMriT2); + if isempty(sMriT1) errMsg = 'SPM Segment failed.'; return; end end catch errMsg = ['bst_normalize_mni/' Method ': ' lasterr()]; - sMri = []; + sMriT1 = []; return; end @@ -177,10 +197,10 @@ %% ===== SAVE RESULTS ===== bst_progress('text', 'Saving normalization...'); % Compute default fiducials positions based on MNI coordinates -sMri = mri_set_default_fid(sMri, 'maff8'); +sMriT1 = mri_set_default_fid(sMriT1, 'maff8'); % Save modifications in the MRI file -if ~isempty(MriFile) - bst_save(file_fullpath(MriFile), sMri, 'v6'); +if ~isempty(T1File) + bst_save(file_fullpath(T1File), sMriT1, 'v6'); end @@ -188,32 +208,32 @@ % If the MRI is currently loaded if ~isempty(iLoadedMri) % Update structures - GlobalData.Mri(iLoadedMri).NCS.R = sMri.NCS.R; - GlobalData.Mri(iLoadedMri).NCS.T = sMri.NCS.T; - GlobalData.Mri(iLoadedMri).NCS.AC = sMri.NCS.AC; - GlobalData.Mri(iLoadedMri).NCS.PC = sMri.NCS.PC; - GlobalData.Mri(iLoadedMri).NCS.IH = sMri.NCS.IH; - GlobalData.Mri(iLoadedMri).NCS.Origin = sMri.NCS.Origin; - if isfield(sMri.NCS,'y') && isfield(sMri.NCS,'iy') && isfield(sMri.NCS,'y_vox2ras') - GlobalData.Mri(iLoadedMri).NCS.y = sMri.NCS.y; - GlobalData.Mri(iLoadedMri).NCS.iy = sMri.NCS.iy; - GlobalData.Mri(iLoadedMri).NCS.y_vox2ras = sMri.NCS.y_vox2ras; + GlobalData.Mri(iLoadedMri).NCS.R = sMriT1.NCS.R; + GlobalData.Mri(iLoadedMri).NCS.T = sMriT1.NCS.T; + GlobalData.Mri(iLoadedMri).NCS.AC = sMriT1.NCS.AC; + GlobalData.Mri(iLoadedMri).NCS.PC = sMriT1.NCS.PC; + GlobalData.Mri(iLoadedMri).NCS.IH = sMriT1.NCS.IH; + GlobalData.Mri(iLoadedMri).NCS.Origin = sMriT1.NCS.Origin; + if isfield(sMriT1.NCS,'y') && isfield(sMriT1.NCS,'iy') && isfield(sMriT1.NCS,'y_vox2ras') + GlobalData.Mri(iLoadedMri).NCS.y = sMriT1.NCS.y; + GlobalData.Mri(iLoadedMri).NCS.iy = sMriT1.NCS.iy; + GlobalData.Mri(iLoadedMri).NCS.y_vox2ras = sMriT1.NCS.y_vox2ras; end - GlobalData.Mri(iLoadedMri).SCS.R = sMri.SCS.R; - GlobalData.Mri(iLoadedMri).SCS.T = sMri.SCS.T; - GlobalData.Mri(iLoadedMri).SCS.NAS = sMri.SCS.NAS; - GlobalData.Mri(iLoadedMri).SCS.LPA = sMri.SCS.LPA; - GlobalData.Mri(iLoadedMri).SCS.RPA = sMri.SCS.RPA; - GlobalData.Mri(iLoadedMri).SCS.Origin = sMri.SCS.Origin; + GlobalData.Mri(iLoadedMri).SCS.R = sMriT1.SCS.R; + GlobalData.Mri(iLoadedMri).SCS.T = sMriT1.SCS.T; + GlobalData.Mri(iLoadedMri).SCS.NAS = sMriT1.SCS.NAS; + GlobalData.Mri(iLoadedMri).SCS.LPA = sMriT1.SCS.LPA; + GlobalData.Mri(iLoadedMri).SCS.RPA = sMriT1.SCS.RPA; + GlobalData.Mri(iLoadedMri).SCS.Origin = sMriT1.SCS.Origin; end %% ===== TISSUE CLASSIFICATION ===== % Import tissue classification -if ~isempty(TpmFiles) && ~isempty(MriFile) +if ~isempty(TpmFiles) && ~isempty(T1File) bst_progress('text', 'Loading tissue segmentations...'); % Get subject - [sSubject, iSubject] = bst_get('MriFile', MriFile); + [sSubject, iSubject] = bst_get('T1File', T1File); % Import tissue classification import_mri(iSubject, TpmFiles, 'SPM-TPM', 0, 1, 'tissues_segment'); end diff --git a/toolbox/anatomy/mri_normalize_segment.m b/toolbox/anatomy/mri_normalize_segment.m index 42d5ca1c5..25c14e39e 100644 --- a/toolbox/anatomy/mri_normalize_segment.m +++ b/toolbox/anatomy/mri_normalize_segment.m @@ -1,7 +1,9 @@ -function [sMri, TpmFiles] = mri_normalize_segment(sMri, TpmFile) +function [sMriT1, TpmFiles] = mri_normalize_segment(sMriT1, TpmFile, sMriT2) % MRI_NORMALIZE_SEGMENT: Non-linear normalization to the MNI ICBM152 space % and tissue segmentation using SPM's Segment batch. % +% USAGE: [sMriT1, TpmFiles] = mri_normalize_segment(sMriT1, TpmFile, sMriT2=[]) +% % The MNI152 space depends on the TPM.nii file given in input: % - Default in SPM12 : IXI549 template @@ -23,25 +25,43 @@ % For more information type "brainstorm license" at command prompt. % =============================================================================@ % -% Authors: Francois Tadel, 2020 +% Authors: Francois Tadel, 2020-2021 +% === PARSE INPUTS === +if (nargin < 3) || isempty(sMriT2) + sMriT2 = []; +end % === SAVE FILES IN TMP FOLDER === % Output variables TpmFiles = []; -% Save source MRI in .nii format -baseName = 'spm_segment.nii'; -NiiFile = bst_fullfile(bst_get('BrainstormTmpDir'), baseName); -out_mri_nii(sMri, NiiFile); +% Save T1 MRI in .nii format +baseName = 'spm_segment_T1.nii'; +T1Nii = bst_fullfile(bst_get('BrainstormTmpDir'), baseName); +out_mri_nii(sMriT1, T1Nii); +% Save T2 MRI in .nii format +if ~isempty(sMriT2) + baseName = 'spm_segment_T2.nii'; + T2Nii = bst_fullfile(bst_get('BrainstormTmpDir'), baseName); + out_mri_nii(sMriT1, T2Nii); +else + T2Nii = []; +end % === RUN SPM SEGMENT === % Disable warnings warning('off', 'MATLAB:RandStream:ActivatingLegacyGenerators'); % Prepare SPM batch -matlabbatch{1}.spm.spatial.preproc.channel.vols = {[NiiFile ',1']}; -matlabbatch{1}.spm.spatial.preproc.channel.biasreg = 0.001; -matlabbatch{1}.spm.spatial.preproc.channel.biasfwhm = 60; -matlabbatch{1}.spm.spatial.preproc.channel.write = [0 0]; +matlabbatch{1}.spm.spatial.preproc.channel(1).vols = {[T1Nii ',1']}; +matlabbatch{1}.spm.spatial.preproc.channel(1).biasreg = 0.001; +matlabbatch{1}.spm.spatial.preproc.channel(1).biasfwhm = 60; +matlabbatch{1}.spm.spatial.preproc.channel(1).write = [0 0]; +if ~isempty(T2Nii) + matlabbatch{1}.spm.spatial.preproc.channel(2).vols = {[T2Nii ',1']}; + matlabbatch{1}.spm.spatial.preproc.channel(2).biasreg = 0.001; + matlabbatch{1}.spm.spatial.preproc.channel(2).biasfwhm = 60; + matlabbatch{1}.spm.spatial.preproc.channel(2).write = [0 0]; +end matlabbatch{1}.spm.spatial.preproc.tissue(1).tpm = {[TpmFile, ',1']}; matlabbatch{1}.spm.spatial.preproc.tissue(1).ngaus = 1; matlabbatch{1}.spm.spatial.preproc.tissue(1).native = [1 0]; @@ -88,11 +108,11 @@ RegInvFile = bst_fullfile(bst_get('BrainstormTmpDir'), ['iy_' baseName]); if ~file_exist(RegFile) || ~file_exist(RegInvFile) disp('BST> SPM Segment failed.'); - sMri = []; + sMriT1 = []; return; end % Import deformation fields -sMri = import_mnireg(sMri, RegFile, RegInvFile, 'segment'); +sMriT1 = import_mnireg(sMriT1, RegFile, RegInvFile, 'segment'); % === LOAD TISSUES === TpmFiles = {... diff --git a/toolbox/process/functions/process_mni_normalize.m b/toolbox/process/functions/process_mni_normalize.m index 08514c14c..da69cd91f 100644 --- a/toolbox/process/functions/process_mni_normalize.m +++ b/toolbox/process/functions/process_mni_normalize.m @@ -2,8 +2,8 @@ % PROCESS_MNI_NORMALIZE: Compute deformation fields to the MNI ICBM152 space. % % USAGE: OutputFiles = process_mni_normalize('Run', sProcess, sInputs) -% [isOk, errMsg] = process_mni_normalize('Compute', MriFile, Method) -% process_mni_normalize('ComputeInteractive', MriFile, Method=[ask], isUnload=1) +% [isOk, errMsg] = process_mni_normalize('Compute', T1File, Method) +% process_mni_normalize('ComputeInteractive', T1File, Method=[ask], isUnload=1) % @============================================================================= % This function is part of the Brainstorm software: @@ -23,7 +23,7 @@ % For more information type "brainstorm license" at command prompt. % =============================================================================@ % -% Authors: Francois Tadel, 2020 +% Authors: Francois Tadel, 2020-2021 eval(macro_method); end @@ -54,6 +54,10 @@ 'maff8', 'segment'}; sProcess.options.method.Type = 'radio_label'; sProcess.options.method.Value = 'maff8'; + % Use T2 when available + sProcess.options.uset2.Comment = 'Use T2 when available ("segment" only)'; + sProcess.options.uset2.Type = 'checkbox'; + sProcess.options.uset2.Value = 0; end @@ -76,6 +80,8 @@ end % Get method Method = sProcess.options.method.Value; + % Use T2 + UseT2 = sProcess.options.uset2.Value; % ===== GET SUBJECT ===== % Get subject @@ -89,12 +95,27 @@ bst_report('Error', sProcess, [], ['No MRI available for subject "' SubjectName '".']); return end - + % Get default MRI as T1 + T1File = sSubject.Anatomy(sSubject.iAnatomy).FileName; + % Look for T2 MRI + T2File = []; + if UseT2 && (length(sSubject.Anatomy) > 1) + iT2 = find(~cellfun(@(c)isempty(strfind(c,'t2')), lower({sSubject.Anatomy.Comment}))); + % Warning when multiple T2 images + if (length(iT2) > 1) + bst_report('Warning', sProcess, [], ['Subject "' sSubject.Name '" has multiple anatomy volumes with tag "T2": not using T2 to avoid confusion.' 10 ... + 'To use a T2 MRI for better volume segmentation, rename or delete the extra T2-labelled files from the subject anatomy.']); + iT2 = []; + end + % Get T2 filename + if ~isempty(iT2) + T2File = sSubject.Anatomy(iT2).FileName; + end + end + % ===== COMPUTE MNI TRANSFORMATION ===== - % Get default MRI - MriFile = sSubject.Anatomy(sSubject.iAnatomy).FileName; % Call normalize function - [sMri, errMsg] = Compute(MriFile, Method); + [sMriT1, errMsg] = Compute(T1File, Method, T2File); % Error handling if ~isempty(errMsg) bst_report('Error', sProcess, [], errMsg); @@ -105,14 +126,22 @@ %% ===== COMPUTE ===== -function [sMri, errMsg] = Compute(MriFile, Method) - [sMri, errMsg] = bst_normalize_mni(MriFile, Method); +function [sMriT1, errMsg] = Compute(T1File, Method, T2File) + % Parse inputs + if (nargin < 3) || isempty(T2File) + T2File = []; + end + % Compute normalization + [sMriT1, errMsg] = bst_normalize_mni(T1File, Method, T2File); end %% ===== COMPUTE/INTERACTIVE ===== -function sMri = ComputeInteractive(MriFile, Method, isUnload) +function sMriT1 = ComputeInteractive(T1File, Method, isUnload, T2File) % Parse inputs + if (nargin < 4) || isempty(T2File) + T2File = []; + end if (nargin < 3) || isempty(isUnload) isUnload = 1; end @@ -130,17 +159,51 @@ 'MNI normalization method', [], {sProcess.options.method.Comment{2,:}, 'Cancel'}, sProcess.options.method.Comment{2,1}); % Cancel if isempty(Method) || strcmpi(Method, 'Cancel') - sMri = []; + sMriT1 = []; return; end end + % For Segment method: look for T2 MRI in the subject + if strcmpi(Method, 'segment') && isempty(T2File) + % Get subject info + sSubject = bst_get('MriFile', T1File); + % Find any possible T2 + iT2 = []; + if (length(sSubject.Anatomy) > 1) + iT2 = find(~cellfun(@(c)isempty(strfind(c,'t2')), lower({sSubject.Anatomy.Comment}))); + % Warning when multiple T2 images + if (length(iT2) > 1) + res = java_dialog('question', ... + ['Subject "' sSubject.Name '" has multiple anatomy volumes with tag "T2".' 10 10 ... + 'If you want to use a T2 MRI for better volume segmentation, ' 10 ... + 'rename or delete the extra T2-labelled files from the subject anatomy.' 10 10 ... + 'If you click OK, no T2 image will be used in the segmentation.' 10 10], 'MNI normalization', [], {'OK', 'Cancel'}, 'OK'); + if isempty(res) || isequal(res, 'Cancel') + return; + end + iT2 = []; + % Confirmation for using T2 image + elseif (length(iT2) == 1) + isT2 = java_dialog('confirm', ... + ['Subject "' sSubject.Name '" seems to include one T2 MRI, named "' sSubject.Anatomy(iT2).Comment '".' 10 10 ... + 'Use this T2 image to refine the volume segmentation?' 10 10], 'MNI normalization'); + if ~isT2 + iT2 = []; + end + end + end + % Get T2 filename + if ~isempty(iT2) + T2File = sSubject.Anatomy(iT2).FileName; + end + end % Open progress bar bst_progress('start', 'MNI normalization', 'Initialization...'); % Call non-interactive function - [sMri, errMsg] = Compute(MriFile, Method); + [sMriT1, errMsg] = Compute(T1File, Method, T2File); % Error handling if ~isempty(errMsg) - sMri = []; + sMriT1 = []; bst_error(errMsg, 'MNI normalization', 0); end % Close progress bar