Updates brain parts masks

Adding the "subcortex wo BG" (#5) and "BG only" (#6) masks
Fixing small error in other mak generation.

Author: ChristophePhillips

eTPM/generate_corticalMask.m

@@ -13,30 +13,35 @@
 % OUTPUT
 % A single file is generated, msk_BrainParts.nii, with 4 different volumes 
 % covering the brain: 
-% - cortex without WM, 
-% - cortex with WM, 
-% - subcortex (pons + basal ganglia), and
-% - cerebellum.
+% 1. cortex without WM, 
+% 2. cortex with WM, 
+% 3. subcortex (Basal Ganglia, Amygdala, brainstem),
+% 4. cerebellum.
+% 5. subcortex without the Basal Ganglia
+% 6. Basal Ganglia only (pallidum, caudate, putamen, nucleaus accumbens,
+%   thamus proper)
 % 
 % NOTE
-% The mask generated include *both* GM ande WM! It is up to the user to
-% also provide a segmented GM or WM map, e.g. c1*.nii and c2*.nii, in order
-% to decide which tissue type to work with.
-% 
+% For an individual subject, the masks generated likely include *both* GM
+% and WM voxels!
+% It is up to the user to also provide a segmented GM or WM map, e.g. 
+% wc1*.nii and wc2*.nii images, in order to specific voxel values in the
+% brain part(s) of interest.
 %_______________________________________________________________________
 % Copyright (C) 2017 Cyclotron Research Centre
 
 % Written by C. Phillips.
 % Cyclotron Research Centre, University of Liege, Belgium
 
-%% DEFINE a few filenames
+%% DEFINE a few filenames & parameters
 fn_atlas = 'labels_Neuromorphometrics.nii';
 fn_parts = 'msk_BrainParts.nii';
 dr_TPM = fullfile(spm('dir'),'tpm'); % SPM's tpm folder
 dr_TPMuswl = fullfile(spm_file(which('tbx_cfg_USwLesion.m'),'path'), ...
     'eTPM'); % USwL's eTPM folder
 mask_smoothing = 4*ones(1,3); % in mm
 sval_thr = .1; % threshold for smoothed volume
+nVols = 6;
 
 %% LOAD the atlas
 Vatlas = spm_vol(fullfile(dr_TPM,fn_atlas));
@@ -59,7 +64,7 @@
 sval_cort_wWM = zeros(SZ);
 spm_smooth(val_cort_wWM,sval_cort_wWM,mask_smoothing./vx_sz); % extend a bit by smoothing
 
-%% SUBCORTICAL part
+%% SUBCORTICAL part: Basal Ganglial, Amygdala, brainstem,...
 l_rois = [ 23 30 31 32 35 36 37 47 48 55 56 57 58 59 60 61 62 75 76 ];
 val_subc = zeros(SZ);
 sval_subc = zeros(SZ);
@@ -82,25 +87,57 @@
 end
 spm_smooth(val_cereb,sval_cereb,mask_smoothing./vx_sz); % extend a bit by smoothing
 
+%% SUBCORTICAL part, without Basal Ganglia
+l_rois = [ 31 32 35 47 48 61 62 75 76 ];
+val_subc_woBG = zeros(SZ);
+sval_subc_woBG = zeros(SZ);
+for ii=l_rois
+    val_subc_woBG = val_subc_woBG +  double(val_atl == ii);
+end
+spm_smooth(val_subc_woBG,sval_subc_woBG,mask_smoothing./vx_sz); % extend a bit by smoothing
+
+%% SUBCORTICAL part, Basal Ganglia only
+% BG area includes these bits
+% - pallidum, #55 and #56
+% - caudate, #36 and #37
+% - putamen, #57 and #58
+% - nucleaus accumbens, #23 and #30
+% - thamus proper, #59 and #60
+l_rois = [55 56 36 37 57 58 23 30 59 60];
+val_subc_BG = zeros(SZ);
+sval_subc_BG = zeros(SZ);
+for ii=l_rois
+    val_subc_BG = val_subc_BG +  double(val_atl == ii);
+end
+spm_smooth(val_subc_BG,sval_subc_BG,mask_smoothing./vx_sz); % extend a bit by smoothing
+
 %% Build the 4D volume
 % Keep a voxel in a mask if sval>.1 and larger than the other 2
-val_parts  = zeros([SZ 3]);
+val_parts  = zeros([SZ nVols]);
 val_parts(:,:,:,1) = sval_cort>sval_thr & ...
     sval_cort>sval_cereb & sval_cort>sval_subc; % -> cortical part
 val_parts(:,:,:,3) = sval_subc>sval_thr & ...
-    sval_subc>sval_cereb & sval_subc>sval_cort; % -> sub-cortical part
+    sval_subc>=sval_cereb & sval_subc>=sval_cort; % -> sub-cortical part
 val_parts(:,:,:,4) = sval_cereb>sval_thr & ...
-    sval_cereb>sval_cort & sval_cereb>sval_subc; % -> cerebellum part
+    sval_cereb>=sval_cort & sval_cereb>=sval_subc; % -> cerebellum part
 
+% Keep a voxel if sval>.1 and not part of cerebelum or sub-cortex
 val_parts(:,:,:,2) = sval_cort_wWM>sval_thr & ...
-    ~val_parts(:,:,:,3) & ~val_parts(:,:,:,3); % -> cortical part with WM
+    ~val_parts(:,:,:,3) & ~val_parts(:,:,:,4); % -> cortical part with WM
+
+% Keep a voxel if sval>.1 and in sub-cortex but not in the other one
+val_parts(:,:,:,5) = val_subc_woBG>sval_thr & ...
+    val_subc_woBG>=val_subc_BG & val_parts(:,:,:,3); % -> sub-cortical wo BG
+val_parts(:,:,:,6) = val_subc_BG>sval_thr & ...
+    val_subc_BG>=val_subc_woBG & val_parts(:,:,:,3); % -> BG only
 
 %% save the 4D volume
 fn_mask = fullfile(dr_TPMuswl,fn_parts);
 % assuming that Vatlas is a uint8 file -> 1 voxel = 1 byte
 descrip = {'Cotical GM mask', 'Cotical GM+WM mask', ...
-            'Subcortical GM mask','Cerebellar GM mask'};
-for ii=1:4
+            'Subcortical GM mask','Cerebellar GM mask', ...
+            'Subcortical wo BG GM mask', 'Basal Ganglia GM mask'};
+for ii=1:nVols
     V_mask(ii) = Vatlas; %#ok<*SAGROW>
     V_mask(ii).pinfo(1) = 1/255;
     V_mask(ii).pinfo(3) = V_mask(ii).pinfo(3) + (ii-1)*prod(SZ);