%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%user specific settings

%would you like to use an mni aligned grid?
mni_aligned_grid='yes';

%path to MRI; MRI should not be co-registered yet
mri_path='....fif';

%path to raw data file with sensor position information
raw_data_path='....fif';

%voxel coordinates of fiducials. can be  taken from the Neuromag GUI for MRI-MEG Integration
%but x and y coordinates need to be swapped!
rpa=[y x z];
nas=[y x z];
lpa=[y x z];

%where would you like to save the results
output_fname='...';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%read mri and raw data file header; paths need to be adapted!
mri  = ft_read_mri(mri_path);
hdr=ft_read_header(raw_data_path);

%determine head coordinate system
cfg=[];
%this is taken from the Neuromag GUI for MRI-MEG Integration. x and y coordinates are swapped!
%values need to be adapted!
cfg.fiducial.rpa=rpa;
cfg.fiducial.nas=nas;
cfg.fiducial.lpa=lpa;
cfg.coordsys='neuromag';
cfg.method='fiducial';
mri=ft_volumerealign(cfg,mri);

%segment the brain into gray, white and csf matter 
cfg                = [];
cfg.keepintermediate = 'no';
cfg.coordsys = 'neuromag';
segmentedmri   = ft_volumesegment(cfg, mri);

%check that segmentation fits the mri; press q to go on
test=segmentedmri;
test.anatomy=mri.anatomy;
cfg=[];
cfg.funparameter='gray';
cfg.interactive='yes';
ft_sourceplot(cfg,test);
clear test

%create single shell, realistic headmodel
cfg                = [];
hdm                = ft_prepare_singleshell(cfg,segmentedmri);

if strcmp(mni_aligned_grid,'yes')
    
%get the individual-to-template transformation matrix and calculate inverse
cfg                = [];
cfg.template       = '/data/apps/spm/spm8/templates/T1.nii'; 
cfg.downsample     = 2;
cfg.coordsys       = 'yes';
cfg.nonlinear      = 'no';
cfg.keepintermediate='no';
norm               = ft_volumenormalise(cfg,mri); 
%this is the template-to-individual matrix
N=inv(norm.cfg.final);

load /data/apps/fieldtrip/templates/template_grid1cm
subj_grid               = [];
%warp mni template grid to individual brain 
subj_grid.pos           = warp_apply(N, template_grid.pos, 'homogenous')./10; 
subj_grid.inside        = template_grid.inside;
subj_grid.outside       = template_grid.outside;
subj_grid.dim           = template_grid.dim;

else
    
    subj_grid.xgrid=-20:1:20;
    subj_grid.ygrid=-20:1:20;
    subj_grid.zgrid=-10:1:20;

end

%to see if everything has worked:
%make a figure of the single subject headmodel, sensor positions and subj_grid positions
cfg                = [];
cfg.vol            = hdm;
cfg.grid           = subj_grid;
cfg.grad       = hdr.grad;
figure;
headmodelplot(cfg);

save(output_fname,'hdm','subj_grid','mri');