Compute average sensor array
The following Matlab code will produce the average sensor array from a list of input files, and select the input file which is closest to the average (which I would recommend to use for interpolation). As input files, you can for example specify the on-line averages for your subjects.
% From a list of fiff-files (e.g. on-line average files), compute the one
% with average sensor transformation, and find the input file which is
% closest to the average (according to both translation and rotation parameters)
% file list should be specified in cell array "files{}"
% output file should be specified in string "file_out"
% both average file "_avg.fif" and file closest to it will be written
% if not output file specified, output will be written to current working directory
% OH, March 2009
% Output file for average sensor array
file_out = '/fullpath/OutputFile.fif';
% you can specify a list of files here (e.g. on-line averages for all subjects in your study)
fiff_files = {'/fullpath/file4subj1.fif', ...
'/fullpath/file4subj2.fif', ...
'/fullpath/file4subj3.fif'};
fid = fopen(file_out, 'a');
if fid==-1,
fprintf(1, 'Cannot access output file %s\n', file_out);
return;
end;
fclose(fid);
nr_files = length(fiff_files);
epoch_gm = [0];
trans_gm = [0];
for f = 1:nr_files, % read all files and average all that's relevant
fprintf(1, '%s\n', fiff_files{f});
data = fiff_read_evoked(fiff_files{f});
epoch_gm = epoch_gm + data.evoked.epochs; % average MEG data epoch
trans_gm = trans_gm + data.info.dev_head_t.trans; % average sensor transformation matrix
all_trans(:,:,f) = data.info.dev_head_t.trans; % remember individual transformations
end;
epoch_gm = epoch_gm/nr_files;
trans_gm = trans_gm/nr_files;
data_new = data;
data_new.evoked.epoch = epoch_gm;
data_new.info.dev_head_t.trans = trans_gm;
[thispath, thisfile,thisext,thisversn] = fileparts(file_out);
file_out_avg = fullfile(thispath, [thisfile '_avg' thisext]);
fiff_write_evoked(file_out_avg, data_new); % Write average coordinate transformation
% find input file that is closest to average...
for f = 1:nr_files,
diff_trans_t = all_trans(1:3,4,f) - trans_gm(1:3,4); % translation difference
diff_trans_r = all_trans(1:3,1:3,f) - trans_gm(1:3,1:3); % rotation difference
norm_t(f) = norm( diff_trans_t );
norm_r(f) = norm( diff_trans_r );
end;
[y_t, i_t] = sort( norm_t ); % sort according to translation difference
[y_r, i_r] = sort( norm_r ); % sort according to rotation difference
for f = 1:nr_files, % combine the two sorted lists
rank_t(f) = find( (i_t-f)==0 ); % rank for translation
rank_r(f) = find( (i_r-f)==0 ); % rank for rotation
avg_rank(f) = (rank_t(f) + rank_r(f))/2;
[min_val, min_file] = min( avg_rank ); % best average rank
end;
min_file = fiff_files{min_file}; % filename for file with best average rank
fprintf('\nThe winner is... %s\n', min_file);
copyfile(min_file, file_out); % copy file with best average rank to output file
% The end of this