Please read this about paths and files before you proceed!

Since the workshop was not intended to cover the preprocessing steps, the raw files were not actually provided.
Therefore, the preprocessing tutorials rather serve the purpose of showing what was done to the preprocessed data.
They of course also serve as examples for you to model your own analysis if deemed useful

Set up general paths

Change these to appropriate paths for your operating system and setup

%% paths

restoredefaultpath
addpath /home/lau/matlab/fieldtrip-20170613_workshop_natmeg/
ft_defaults

raw_meg_path = '/archive/20067_workshop_source_reconstruction/MEG/';
meg_path = '/home/share/workshop_source_reconstruction/data/MEG/';
mri_path = '/home/share/workshop_source_reconstruction/data/MRI/';

Set up subject specific paths

Make subject and recording specific paths (the cell array “subjects_and_dates” can be expanded)

%% subjects_and_dates

subjects_and_dates = ...
                    {
                        'NatMEG_0177/170424/'
                    };
                
subjects = ...
            {
                'NatMEG_0177'
            };
                
filenames = {...
        'tactile_stim_raw_tsss_mc.fif'
        'tactile_stim_raw_tsss_mc-1.fif'
        'tactile_stim_raw_tsss_mc-2.fif'
            };
        
dicom_name = '00000001.dcm';

output_path = fullfile(meg_path, subjects_and_dates{1});

Read in dicom files

%% read in mri

dicom_path = fullfile(mri_path, 'dicoms', subjects{1}, dicom_name);

mri = ft_read_mri(dicom_path);

save([output_path 'mri'], 'mri');

Indicate coordinate system

%% determine coordinate system

mri_coordsys = ft_determine_coordsys(mri);

save([output_path 'mri_coordsys'], 'mri_coordsys');

Co-register mri to neuromag based on fiducials (nasion and left and right pre-auricular points)

%% co-register fieldtrip and neuromag

cfg = [];
cfg.method = 'interactive';
cfg.coordsys = 'neuromag';

mri_realigned_1 = ft_volumerealign(cfg, mri_coordsys);

save([output_path 'mri_realigned_1'], 'mri_realigned_1');

Co-register with the extra digitization points (along the scalp and the face)

An iterative closest point (icp) procedure is used

%% co-register to the remaining points

raw_filepath = fullfile(raw_meg_path, subjects_and_dates{1}, filenames{1});
input_path  = fullfile(meg_path, subjects_and_dates{1}, filenames{1});

% headshape = ft_read_headshape(raw_filepath);
headshape = ft_read_headshape(input_path);
cfg = [];
cfg.method = 'headshape';
cfg.headshape.headshape = headshape;
cfg.headshape.icp = 'yes';
cfg.coordsys = 'neuromag';

mri_realigned_2 = ft_volumerealign(cfg, mri_realigned_1);

save([output_path 'mri_realigned_2'], 'mri_realigned_2');

Check co-registration

This step is really just used to check whether the coregistration went okay

%% checking co-registration

cfg = [];
cfg.method = 'headshape';
cfg.headshape.headshape = headshape;
cfg.coordsys = 'neuromag';
cfg.headshape.icp = 'no';

mri_realigned_3 = ft_volumerealign(cfg, mri_realigned_2);

save([output_path 'mri_realigned_3'], 'mri_realigned_3');

Reslice MRI data

%% reslice mri

cfg = [];
cfg.resolution = 1;

mri_resliced = ft_volumereslice(cfg, mri_realigned_3);

save([output_path 'mri_resliced'], 'mri_resliced');

mri_resliced_cm = ft_convert_units(mri_resliced, 'cm');

save([output_path 'mri_resliced_cm'], 'mri_resliced_cm');

Segment data into brain, skull and scalp

%% segment data

cfg = [];
cfg.output = {'brain' 'skull' 'scalp'};

mri_segmented = ft_volumesegment(cfg, mri_resliced);

save([output_path 'mri_segmented'], 'mri_segmented');

Correct segmentation

%% correct segmentation see http://www.fieldtriptoolbox.org/tutorial/natmeg/dipolefitting for full details

binary_brain = mri_segmented.brain;
binary_skull = mri_segmented.skull | binary_brain;
binary_scalp = mri_segmented.scalp | binary_brain | binary_skull;
binary_scalp = mri_segmented.scalp + binary_skull;

% use boolean logic together with IMERODE
binary_skull = binary_skull & imerode(binary_scalp, strel_bol(2)); % fully contained inside eroded scalp
binary_brain = binary_brain & imerode(binary_skull, strel_bol(2)); % fully contained inside eroded skull

Create a new segmented mri

%% insert mri_segmented 2

mri_segmented_2 = mri_segmented;
% insert the updated binary volumes, taking out the center part for skull and scalp
mri_segmented_2.brain    = binary_brain;
mri_segmented_2.skull    = binary_skull & ~binary_brain;
mri_segmented_2.scalp    = binary_scalp & ~binary_brain & ~binary_skull;

Constuct meshes for each segment

cfg = [];
cfg.method = 'projectmesh';
cfg.tissue = 'brain';
cfg.numvertices = 3000;

mesh_brain = ft_prepare_mesh(cfg, mri_segmented_2);

cfg.tissue = 'skull';
cfg.numvertices = 2000;

mesh_skull = ft_prepare_mesh(cfg, mri_segmented_2);

cfg.tissue = 'scalp';
cfg.numvertices = 1000;

mesh_scalp = ft_prepare_mesh(cfg, mri_segmented_2);

meshes = [mesh_brain mesh_skull mesh_scalp];

save([output_path 'meshes'], 'meshes');

Construct separate head models for EEG and MEG

Head models define how electric potentials and magnetic fields spread throughout the conductor (the head)
Note the simplicity of the MEG model compared to the EEG model

cfg = [];
cfg.method = 'bemcp';
cfg.conductivity = [1 1/20 1] .* (1/3);

headmodel_eeg = ft_prepare_headmodel(cfg, meshes);

cfg = [];
cfg.method = 'singleshell';

headmodel_meg = ft_prepare_headmodel(cfg, mesh_brain);

save([output_path 'headmodel_eeg'], 'headmodel_eeg');
save([output_path 'headmodel_meg'], 'headmodel_meg');