fix conflict

analysis/beamformer.m

@@ -2,7 +2,7 @@
 clear;
 
 % Flag to indicate whether we update source and leadfield models
-update_models = false;
+update_models = true;
 
 % Add util and template dirs to path
 addpath('/project/3031004.01/meg-ahat/util')
@@ -303,13 +303,13 @@
 end
 
 %% Average over lateral contrasts
-close all
+
 lateral_dif_sources = cell(1,numel(subjects));
-conditions = ["con" "strobe"];
-title_contrast = ["Left minus right attention" "High minus low arithmetic difficulty"];
 deriv_anat_dir = fullfile(derivatives_dir, 'sub-030', '/ses-001/anat/'); % Use sub 30 mri for now
 mri_realigned_file = fullfile(deriv_anat_dir, 'mri_realigned.mat');
 load (mri_realigned_file)
+title_contrast = ["Left minus right attention" "High minus low arithmetic difficulty"];
+
 for task_no = 1:numel(tasks)
     task = tasks(task_no)
     for condition = conditions

analysis/make_head_model.m

@@ -0,0 +1,169 @@
+%% Add util dir to path
+clear;
+
+addpath('/project/3031004.01/meg-ahat/util')
+addpath('/project/3031004.01/meg-ahat/analysis')
+
+% Define directories
+data_dir = '/project/3031004.01/data/';
+raw2_dir = fullfile(data_dir, 'raw2');
+derivatives_dir = fullfile(data_dir, 'derivatives');
+
+% Start logging
+diaryfile = fullfile(data_dir, 'make_head_model.log');
+if (exist(diaryfile, 'file'))
+  delete(diaryfile);
+end
+diary (diaryfile)
+
+% Set up Fieldtrip
+configure_ft
+
+% Load data details
+data_details_cfg = get_data_details();
+
+% Define subjects - should eventually be centralised
+% subjects = data_details_cfg.new_trigger_subs; %[8 9 11 13 17 18 21:23 25 27:30];
+subjects = [24 26]
+
+for sub = subjects
+    close all
+    sub_str = sprintf('sub-%03d', sub)
+
+    deriv_anat_dir = fullfile(derivatives_dir, sub_str, '/ses-001/anat/');
+    raw2_meg_dir = fullfile(raw2_dir, sprintf('sub-%03d', sub), '/ses-001/meg/');
+
+    if not(exist(deriv_anat_dir, 'dir'))
+      mkdir(deriv_anat_dir);
+    end
+
+    % Check if fiducials/realignment has already been done
+    mri_file = fullfile(deriv_anat_dir, 'mri_realigned.mat');
+    if isfile(mri_file)
+        load (mri_file)
+    else
+
+        % Load anatomical MRI data
+        anat_dir = fullfile(raw2_dir, sub_str, 'ses-001', 'anat');
+        fname = sprintf('%s_ses-001_T1w.nii', sub_str);
+        mri = ft_read_mri(fullfile(anat_dir, fname));
+        mri = ft_determine_coordsys(mri);
+
+        % Re-align MRI to MEG
+        cfg = [];
+        cfg.method = 'interactive';
+        cfg.coordsys = 'ctf';
+        cfg.viewresult = 'no';
+        mri_realigned = ft_volumerealign(cfg, mri); 
+        save (mri_file, 'mri_realigned', '-v7.3')
+    end
+
+    % Segment MRI
+    cfg           = [];
+    cfg.output    = 'brain';
+    cfg.spmmethod = 'old';
+    cfg.spmversion = 'spm12';
+    segmentedmri  = ft_volumesegment(cfg, mri_realigned);
+
+    % Build head model
+    cfg = [];
+    cfg.method='singleshell';
+    mri_headmodel = ft_prepare_headmodel(cfg, segmentedmri);
+
+    % Visualise
+    mri_headmodel = ft_convert_units(mri_headmodel, 'mm');
+    sens = ft_read_sens(fullfile(raw2_meg_dir, ...
+       sprintf('%s_ses-001_task-flicker_meg.ds', sub_str)), 'senstype', 'meg');
+
+    figure
+    title(sub_str)
+    ft_plot_sens(sens, 'style', '*b');
+
+    hold on
+    ft_plot_headmodel(mri_headmodel);
+
+    % Save headmodel to disk
+    mri_headmodel_file = fullfile(deriv_anat_dir, 'mri_headmodel.mat');
+    save (mri_headmodel_file, 'mri_headmodel', '-v7.3')
+
+end
+
+diary off
+% 
+% %% ARCHIVE
+% %% Add util dir to path
+% addpath('/project/3031004.01/meg-ahat/util')
+% 
+% %% Metedata
+% % Sebject (defined as cell array for future compatibility)
+% subj = {
+%   '099'
+%   };
+% ses = 1;
+% 
+% %% Set pilot data directory
+% data_dir = '/project/3031004.01/pilot-data';
+% diaryfile = strcat(data_dir, '/make_head_model.log');
+% 
+% if (exist(diaryfile, 'file'))
+%   delete(diaryfile);
+% end
+% 
+% diary (diaryfile)
+% 
+% raw2_dir = '/project/3031004.01/pilot-data/raw2';
+% derivatives_dir = '/project/3031004.01/pilot-data/derivatives';
+% deriv_anat_dir = sprintf('%s/sub-%s/ses-00%d/anat/', derivatives_dir, subj{1}, ses);
+% 
+% %% Se tup Fieldtrip
+% configure_ft
+% 
+% %% Load anatomical MRI data
+% anat_dir = sprintf('%s/sub-%s/ses-00%d/anat/', raw2_dir, subj{1}, ses);
+% fname = sprintf('sub-%s_ses-00%d_T1w.nii', subj{1}, ses);
+% mri = ft_read_mri(strcat(anat_dir, fname));
+% %mri.coordsys = "ctf"; %% OBS: I added this, though it was implicit in the data in tutorial
+% mri = ft_determine_coordsys(mri);
+% %% Re-align MRI to MEG
+% cfg = [];
+% cfg.method = 'interactive';
+% cfg.coordsys = 'ctf';
+% cfg.viewresult = 'yes';
+% mri_realigned = ft_volumerealign(cfg, mri); 
+% fname = sprintf('mri_realigned.mat', subj{1}, ses);
+% % use fullfile instead
+% save (sprintf('%s%s', deriv_anat_dir, fname), 'mri_realigned', '-v7.3');
+% 
+% %% Segment MRI
+% cfg           = [];
+% cfg.output    = 'brain';
+% segmentedmri  = ft_volumesegment(cfg, mri_realigned);
+% 
+% %% Build head model
+% cfg = [];
+% cfg.method='singleshell';
+% mri_headmodel = ft_prepare_headmodel(cfg, segmentedmri);
+% 
+% %% Save segmented MRI to derivatives
+% 
+% if not(exist(deriv_anat_dir, 'dir'))
+%   mkdir(deriv_anat_dir);
+% end
+% 
+% fname = sprintf('sub-%s_ses-00%d_proc-segmented_T1w.mat', subj{1}, ses);
+% save (sprintf('%s%s', deriv_anat_dir, fname), 'mri_headmodel', '-v7.3');
+% 
+% 
+% diary off
+% 
+% %% Visualise
+% mri_headmodel = ft_convert_units(mri_headmodel, 'cm');
+% %meg_dir = sprintf('%s/sub-%s/ses-00%d/meg/', raw2_dir, subj{1}, ses);
+% %fname = sprintf('sub-%s_ses-00%d_T1w.nii', subj{1}, ses);
+% sens = ft_read_sens(sprintf('%s/sub-%s/ses-00%d/meg/sub-%s_ses-00%d_task-flicker_meg.ds', raw2_dir, subj{1}, ses, subj{1}, ses), 'senstype', 'meg');
+% 
+% figure
+% ft_plot_sens(sens, 'style', '*b');
+% 
+% hold on
+% ft_plot_headmodel(mri_headmodel);

util/get_data_details.m

@@ -0,0 +1,110 @@
+function data_details_cfg = get_data_details()
+    data_details_cfg = [];
+
+    all_subs = 1:30;
+
+    data_details_cfg.old_trigger_subs = [1:7 10 12 14:17 19:20];
+    data_details_cfg.new_trigger_subs = setdiff(all_subs, data_details_cfg.old_trigger_subs);
+    assert(length(data_details_cfg.new_trigger_subs) == length(data_details_cfg.old_trigger_subs))
+
+    generic_strct = [];
+    generic_strct.run1 = [];
+    generic_strct.run2 = [];
+
+    for sub = all_subs
+        data_details_cfg.(sprintf('sub%03d', sub)) = [];
+        data_details_cfg.(sprintf('sub%03d', sub)).run1 = [];
+        data_details_cfg.(sprintf('sub%03d', sub)).run2 = [];
+    end
+
+    data_details_cfg.sub001.run1.suffix = 'a042d115'; % commit ID
+    data_details_cfg.sub001.run2.suffix = 'a042d115';
+
+    data_details_cfg.sub002.run1.suffix = 'a042d115';
+    data_details_cfg.sub002.run2.suffix = 'a042d115';
+
+    data_details_cfg.sub003.run1.suffix = 'a042d115';
+    data_details_cfg.sub003.run2.suffix = 'a042d115';
+
+    data_details_cfg.sub004.run1.suffix = 'a042d115';
+    data_details_cfg.sub004.run2.suffix = 'a042d115';
+
+    data_details_cfg.sub005.run1.suffix = 'a042d115';
+    data_details_cfg.sub005.run2.suffix = 'a042d115';
+
+    data_details_cfg.sub006.run1.suffix = 'b995195c';
+    data_details_cfg.sub006.run2.suffix = 'b995195c';
+
+    data_details_cfg.sub007.run1.suffix = '234e5f40';
+    data_details_cfg.sub007.run2.suffix = 'b2aef43d';
+
+    data_details_cfg.sub008.run1.suffix = '00e480a9';
+    data_details_cfg.sub008.run2.suffix = '00e480a9';
+
+    data_details_cfg.sub009.run1.suffix = 'af7905d8';
+    data_details_cfg.sub009.run2.suffix = 'af7905d8';
+
+    data_details_cfg.sub010.run1.suffix = 'cc66c65b';
+    data_details_cfg.sub010.run2.suffix = 'cc66c65b';
+
+    data_details_cfg.sub011.run1.suffix = '00e480a9';
+    data_details_cfg.sub011.run2.suffix = '00e480a9';
+
+    data_details_cfg.sub012.run1.suffix = 'd70eacfb';
+    data_details_cfg.sub012.run2.suffix = 'd70eacfb';
+
+    data_details_cfg.sub013.run1.suffix = '00e480a9';
+    data_details_cfg.sub013.run2.suffix = '00e480a9';
+
+    data_details_cfg.sub014.eyetrack_missing = true;
+    data_details_cfg.sub014.run1.suffix = 'cc66c65b';
+    data_details_cfg.sub014.run2.suffix = 'cc66c65b';
+
+    data_details_cfg.sub015.run1.suffix = 'fbd05d5a';
+    data_details_cfg.sub015.run2.suffix = 'fbd05d5a';
+
+    data_details_cfg.sub016.run1.suffix = 'cc66c65b';
+    data_details_cfg.sub016.run2.suffix = 'cc66c65b';
+
+    data_details_cfg.sub017.run1.suffix = 'fbd05d5a';
+    data_details_cfg.sub017.run2.suffix = 'fbd05d5a';
+
+    data_details_cfg.sub018.run1.suffix = '00e480a9';
+    data_details_cfg.sub018.run2.suffix = '00e480a9';
+
+    data_details_cfg.sub019.run1.suffix = 'cc66c65b_00'; % commit ID + overwrite protection index
+    data_details_cfg.sub019.run2.suffix = 'cc66c65b';
+
+    data_details_cfg.sub020.run1.suffix = 'cc66c65b';
+    data_details_cfg.sub020.run2.suffix = 'cc66c65b';
+
+    data_details_cfg.sub021.run1.suffix = '00e480a9';
+    data_details_cfg.sub021.run2.suffix = '00e480a9';
+
+    data_details_cfg.sub022.run1.suffix = '00e480a9';
+    data_details_cfg.sub022.run2.suffix = '00e480a9';
+
+    data_details_cfg.sub023.run1.suffix = '00e480a9';
+    data_details_cfg.sub023.run2.suffix = '00e480a9';
+
+    data_details_cfg.sub024.run1.suffix = '7e9fd61a';
+    data_details_cfg.sub024.run2.suffix = '7e9fd61a';
+
+    data_details_cfg.sub025.run1.suffix = '00e480a9';
+    data_details_cfg.sub025.run2.suffix = '00e480a9';
+
+    data_details_cfg.sub026.run1.suffix = '7e9fd61a';
+    data_details_cfg.sub026.run2.suffix = '7e9fd61a';
+
+    data_details_cfg.sub027.run1.suffix = '00e480a9';
+    data_details_cfg.sub027.run2.suffix = '00e480a9';
+
+    data_details_cfg.sub028.run1.suffix = '87a85519';
+    data_details_cfg.sub028.run2.suffix = '87a85519';
+
+    data_details_cfg.sub029.run1.suffix = '00e480a9';
+    data_details_cfg.sub029.run2.suffix = '00e480a9';
+
+    data_details_cfg.sub030.run1.suffix = 'af7905d8';
+    data_details_cfg.sub030.run2.suffix = 'af7905d8';
+end

util/make_mask.m

@@ -0,0 +1,55 @@
+function mask = make_mask(dat, thr)
+% This function generates an opacity-ramping mask with a decent plateau and
+% nice ramp-on/ramp-off. Useful for presenting e.g. t-maps thresholded at
+% some % of maximum.
+%
+% Copyright (C) Eelke Spaak, Donders Institute, Nijmegen, The Netherlands, 2019.
+% 
+% This code is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+% 
+% This code is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+% GNU General Public License for more details.
+% 
+% You should have received a copy of the GNU General Public License
+% along with this code. If not, see <https://www.gnu.org/licenses/>.
+%
+% MODIFIED by Mark Henney May-2024 to reinstate sign of data after masking
+
+if nargin < 2
+  thr = [0.5 0.8];
+end
+
+% if data contains negatives, ensure strongly negative values are given as
+% much opacity as strongly positive ones
+dat_sign = sign(dat); % Save the sign for later
+dat = abs(dat);
+
+mask = zeros(size(dat));
+
+thr = max(dat, [], 'all') .* thr; % Modify args for max() to all dims
+
+% everything above thr(2) is fully opaque
+mask(dat > thr(2)) = 1;
+
+% in between thr(1) and thr(2): ramp up nicely
+inds = dat > thr(1) & dat < thr(2);
+x = dat(inds);
+
+% scale between 0 and 1
+x = (x-min(x)) ./ (max(x)-min(x));
+
+% make sigmoidal
+beta = 2;
+x = 1 ./ (1 + (x./(1-x)).^-beta);
+
+mask(inds) = x;
+
+% Reinstate sign
+mask = mask .* dat_sign;
+
+end