Diff for "Maxfilter_V2.2" - Meg Wiki
location: Diff for "Maxfilter_V2.2"
Differences between revisions 21 and 23 (spanning 2 versions)
Revision 21 as of 2014-07-24 11:12:57
Size: 21770
Editor: OlafHauk
Comment:
Revision 23 as of 2019-03-25 17:56:54
Size: 20336
Editor: OlafHauk
Comment:
Deletions are marked like this. Additions are marked like this.
Line 4: Line 4:
[[attachment:Maxfilter_Manual_v2pt2.pdf|Maxfilter 2.2 manual.]]
Line 12: Line 14:
For more information on these and other options, see the [[attachment:Maxfilter_Manual_v2pt2.pdf|Maxfilter 2.2 Manual]]. For more information on these and other options, see the Maxfilter 2.2 Manual.
Line 22: Line 24:
Line 37: Line 38:
% Note that have turned off "movecomp" option via "mvcomp_fail" matrix, because  % Note that have turned off "movecomp" option via "mvcomp_fail" matrix, because
Line 60: Line 61:
cbu_codes{1} = {     cbu_codes{1} = {
Line 71: Line 72:
    'meg14_0028' 
    'meg14_108' 
    'meg14_0028'
    'meg14_108'
Line 75: Line 76:
 'meg14_0117'         'meg14_0117'
Line 95: Line 96:
runs{3} = {''};  runs{3} = {''};
Line 121: Line 122:
% user_bad{1}{3} = [1218 1278];   % user_bad{1}{3} = [1218 1278];
Line 125: Line 126:
basestr = [basestr ' -force'];   basestr = [basestr ' -force'];
Line 133: Line 134:
%mvcomp_fail(1,1,1,2) = 1; % Pilot1 object second run 
%mvcomp_fail(1,1,2,2) = 1; % Pilot1 scene second run 
%mvcomp_fail(1,1,1,2) = 1; % Pilot1 object second run
%mvcomp_fail(1,1,2,2) = 1; % Pilot1 scene second run
Line 158: Line 159:
TransDefaultFlag = 1;  TransDefaultFlag = 1;
Line 168: Line 169:
            
Line 170: Line 171:
            for e = 1:length(expts)  
                
            for e = 1:length(expts)
Line 173: Line 174:
                
Line 177: Line 178:
                
Line 179: Line 180:
                
                for r = 1:length(runs{e})                        
                    

                for r = 1:length(runs{e})
Line 183: Line 184:
                    
Line 190: Line 191:
                    
Line 198: Line 199:
                        
Line 208: Line 209:
                        
Line 210: Line 211:
                          
Line 213: Line 214:
                        
Line 216: Line 217:
                        
Line 220: Line 221:
                        
Line 224: Line 225:
                        
Line 228: Line 229:
                        
Line 232: Line 233:
                        
Line 235: Line 236:
                        
Line 243: Line 244:
                                                 
Line 246: Line 247:
                        
Line 248: Line 249:
                        
Line 254: Line 255:
                        
Line 262: Line 263:
                        
Line 264: Line 265:
                        
Line 268: Line 269:
                        
Line 271: Line 272:
                                                 
Line 274: Line 275:
                        
Line 280: Line 281:
                            rik_eval(finstr);                                                        rik_eval(finstr);
Line 284: Line 285:
                        
Line 289: Line 290:
        
Line 291: Line 292:
        
Line 295: Line 296:
            
Line 297: Line 298:
                
Line 299: Line 300:
                
Line 303: Line 304:
                
Line 305: Line 306:
                
Line 307: Line 308:
                    
Line 309: Line 310:
                    
Line 316: Line 317:
                    
Line 323: Line 324:
                    else                                            else
Line 333: Line 334:
                        
Line 335: Line 336:
                                                 
Line 338: Line 339:
                        
Line 341: Line 342:
                        
Line 345: Line 346:
                        
Line 353: Line 354:
                        
Line 357: Line 358:
                        
Line 361: Line 362:
                        
Line 364: Line 365:
                        
Line 372: Line 373:
                                                 
Line 375: Line 376:
                        
Line 377: Line 378:
                        
Line 383: Line 384:
                        
Line 391: Line 392:
                         
Line 393: Line 394:
                        
Line 401: Line 402:
                        
Line 404: Line 405:
                                                 
Line 407: Line 408:
                        
Line 417: Line 418:
                            end                                                        end
Line 420: Line 421:
                        end                                                 end

[ATTACH]

Maxfilter Version 2.2

Maxfilter 2.2 manual.

Example

Version 2.2. is currently the default on our 64-bit machines (see below for more details).

This command will apply Maxfilter including Signal Space Separation (SSS), its temporal extension (ST), and movement compensation.

maxfilter  -f input_file.fif  -o output_file.fif  -st  -origin 0 0 45  -frame head  -autobad on  -movecomp  -format short  -force

For more information on these and other options, see the Maxfilter 2.2 Manual.

Known problems and bugs

Maxfilter 2.2. cannot do downsampling and filtering at the same time as movement compensation. You will therefore have to run those in a separate step, e.g.

maxfilter  -f input_file.fif  -o output_file.fif  -ds 4  -lpfilt 40

Diagnostics

You can get information about movement parameters using these scripts to analyse maxfilter output.

Example Matlab Script

Courtesy of Rik Henson

% Maxfilter 2.2 Matlab script for AD project (R Henson Jan 2013)
%
% This version has three options for running on our Compute machines, based
% on setting of ParType variable below: 0 = run on Login in serial; 1 = run
% one Compute core (using spmd); 2 = run on multiple Compute cores
% (using parfor) - with one core per subject (though could switch parfor
% loop across experiments if running one subject (NB: cannot parfor across
% runs, because they depend on same origin fitting)
%
% Note that have turned off "movecomp" option via "mvcomp_fail" matrix, because
% this caused apparent random MF2.2 crashes for the pilot subject. Could
% set to 0 for next subject, in case specific to pilot.
%
% This version does 2-3 steps, with 1-2 outputs:
%
%   1. SSS with autobad on, to detect bad channels and write out move parameters
%   2. tSSS, downsample by factor 4 (250Hz), and trans to first file if multiple runs
%   3. trans to default helmet position (if flag below set) (note trans'ing
%   direct to default space without realigning runs within subject does not
%   work so well)

clear  % May not want/need to!

dat_wd = '/megdata/cbu/multi_oc';
bas_wd = '/imaging/rh01/Projects/AD_MCI/MEG'; % <-- change to yours

groups = {
    'Pilots'
    'OC1'    % Old controls who did scene-object location task
    'OC2'    % Old controls who did object-quadrant task
    };

cbu_codes{1} = {
    'meg13_0515'  % Pilot1 = Alex; some session names below may not apply)
    'meg14_0003'  % a young control
    };

cbu_codes{2} = {
    'meg14_0008'
    'meg14_0010'
    'meg14_0014'
    'meg14_0018'
%    'meg14_0024' % Two runs of this subject don't maxfilter because of HPI problems (so handled by "badrun" variable below)
    'meg14_0028'
    'meg14_108'
    'meg14_0111'
    'meg14_0113'
        'meg14_0117'
    };

 cbu_codes{3} = { % Object-quadrant task
    'meg14_0031'  % Becki wonders about this one?
    'meg14_0032'
    'meg14_0037'
    'meg14_0043'
    'meg14_0045'
   };


% Experiments
expts = {'object';'scene';'restclosed'; % 'restopen'
    };

% Run labels per experiment (not actually needed with new filenames from Pilot2 onwards)
runs{1} = {'_1'; '_2'};
%runs{2} =  {'_3' '_4'}; % This was only Pilot1
runs{2} = {'_1'; '_2'};
runs{3} = {''};

%trans_run = [2 2 1];  % Which run WITHIN an experiment to trans too (1 for rest if only one rest run)
%trans_run = [0 0 0];   % If want to skip this step, eg trans default only

%% Bad runs:
badrun = zeros(5,50,3,2);  % Assume all runs ok, unless indicated next

subid = find(strcmp(cbu_codes{2},'meg14_0024'));
badrun(2,subid,1,2) = 1; % meg14_0024, object, run 2 cannot be maxfiltered "The origin is only 4 cm from nearest coil!"
badrun(2,subid,2,1) = 1; % meg14_0024, scene,  run 1 cannot be maxfiltered "The origin is only 4 cm from nearest coil!"


% Set up directory structures (only needs to be done once)
for e=1:length(expts)
    eval(sprintf('!mkdir %s',fullfile(bas_wd,expts{e})));
    for g=1:length(groups)
        eval(sprintf('!mkdir %s',fullfile(bas_wd,expts{e},groups{g})));
        for s = 1:length(cbu_codes{g})
%            eval(sprintf('!mkdir %s',fullfile(bas_wd,expts{e},groups{g},sprintf('Sub%02d',s))));
            eval(sprintf('!mkdir %s',fullfile(bas_wd,expts{e},groups{g},cbu_codes{g}{s})));
        end
    end
end

% Any use bad channels? (This option not implemented yet)
% user_bad{1}{3} = [1218 1278];

basestr = ' -ctc /neuro/databases/ctc/ct_sparse.fif -cal /neuro/databases/sss/sss_cal.dat';
basestr = [basestr ' -linefreq 50 -hpisubt amp'];
basestr = [basestr ' -force'];
maxfstr = '!/neuro/bin/util/x86_64-pc-linux-gnu/maxfilter-2.2 '

addpath /imaging/local/meg_misc
addpath /neuro/meg_pd_1.2/

% files where MF2.2 fails when mvcomp included (for some reason) - so turned off below
%mvcomp_fail = zeros(5,50,3,2);
%mvcomp_fail(1,1,1,2) = 1; % Pilot1 object second run
%mvcomp_fail(1,1,2,2) = 1; % Pilot1 scene second run
mvcomp_fail = ones(5,50,3,2);  % Turn off all mvcomp, since seems to fail randomly!

movfile = 'trans_move.txt'; % This file will record translations between runs

%ParType = 0;  % Fun on Login machines (not generally advised!)
%ParType = 1;   % Run maxfilter call on Compute machines using spmd (faster)
ParType = 2;   % Run on multiple Compute machines using parfar (best, but less feedback if crashes)

%% open matlabpool if required
% matlabpool close force CBU_Cluster
if ParType
    if matlabpool('size')==0;
        MaxNsubs = 1;
        if ParType == 2
            for g=1:length(cbu_codes)
                MaxNsubs = max([MaxNsubs length(cbu_codes{g})]);
            end
        end
        P = cbupool(MaxNsubs);
        matlabpool(P);
    end
end

TransDefaultFlag = 1;


%% Main loop (can comment/uncomment lines below if want to parallelise over expts rather than subjects)
for g = 1:length(groups)
    fprintf('\n\n%s\n\n',groups{g})
    if ParType == 2 % parfor loop
        parfor s = 1:length(cbu_codes{g})
            raw_wd    = dir(fullfile(dat_wd,cbu_codes{g}{s},'1*'));  % Just to get date directory (assuming between 20*1*0 and 20*1*9Q!)
            raw_wd    = raw_wd.name;

%        parfor e = 1:length(expts)  % If doing a single subject (note: cannot embed parfor loops unfortunately)
            for e = 1:length(expts)

                transfstfile = ''; orig = []; rad=[]; fit=[]; transtr = {}; raw_stem = {};

                % Output directory
                sub_wd = fullfile(bas_wd,expts{e},groups{g},cbu_codes{g}{s}), cd(sub_wd)
                %try eval(sprintf('!mkdir %s',sub_wd)); end  % Try not allowed n parfor, so make directories in advance above

                rik_eval(sprintf('!touch %s',movfile));

                for r = 1:length(runs{e})

                    raw_file = dir(fullfile(dat_wd,cbu_codes{g}{s},raw_wd,sprintf('%s*%s_raw*',expts{e},runs{e}{r})));  % Get raw FIF file

                    if isempty(raw_file)
                        error('Could not find run %d for grp %s, sub %s, exp %s',r,groups{g},cbu_codes{g}{s},expts{e})
                    else
                        raw_stem = raw_file.name(1:(end-4));
                        raw_file = fullfile(dat_wd,cbu_codes{g}{s},raw_wd,raw_file.name);
                    end

                    if badrun(g,s,e,r)  % Care: need to note that this run not trans'ed (despite name), not sss, etc
                        outfile = fullfile(sub_wd,sprintf('%s_NoSSS',raw_stem));
                        filestr = sprintf(' -f %s -o %s.fif',raw_file,outfile);
                        finstr = [maxfstr filestr ' -nosss -ds 4' sprintf(' -v | tee %s.log',outfile)];
                        rik_eval(finstr);
                        %rik_eval(sprintf('!cp %s %s',raw_file,outfile));
                    else

                        %% Fit sphere (since better than MaxFilter does)
                        if r == 1  % fit sphere doesn't change with run!
                            incEEG = 0;
                            if exist(fullfile(sub_wd,'fittmp.txt')); delete(fullfile(sub_wd,'fittmp.txt')); end
                            if exist(fullfile(sub_wd,sprintf('run_%02d_hpi.txt',r))); delete(fullfile(sub_wd,sprintf('run_%02d_hpi.txt',r)));  end
                            [orig,rad,fit] = meg_fit_sphere(raw_file,sub_wd,sprintf('%s_hpi.txt',raw_stem),incEEG);
                            delete(fullfile(sub_wd,'fittmp.txt'));
                        end
                        origstr = sprintf(' -origin %d %d %d -frame head',orig(1),orig(2),orig(3))

                        badstr  = sprintf(' -autobad %d -badlimit %d',900,7); % 900s is 15mins - ie enough for whole recording!


                        %% 1. Bad channel detection (this email says important if doing tSSS later https://www.jiscmail.ac.uk/cgi-bin/webadmin?A2=NEUROMEG;d3f363f3.1205)

                        outfile = fullfile(sub_wd,sprintf('%s_bad',raw_stem));
                        filestr = sprintf(' -f %s -o %s.fif',raw_file,outfile);

                        % Write out movements too...
                        posfile = fullfile(sub_wd,sprintf('%s_headpos.txt',raw_stem));
                        compstr = sprintf(' -headpos -hpistep 10 -hp %s',posfile);

                        finstr = [maxfstr filestr origstr basestr badstr compstr sprintf(' -v | tee %s.log',outfile)]
                        rik_eval(finstr);
                        delete(sprintf('%s.fif',outfile));

                        % Pull out bad channels from logfile:
                        badfile = sprintf('%s.txt',outfile); delete(badfile);
                        rik_eval(sprintf('!cat %s.log | sed -n -e ''/Detected/p'' -e ''/Static/p'' | cut -f 5- -d '' '' > %s',outfile,badfile));

                        tmp=dlmread(badfile,' '); Nbuf = size(tmp,1);
                        tmp=reshape(tmp,1,prod(size(tmp)));
                        tmp=tmp(tmp>0); % Omit zeros (padded by dlmread):

                        % Get frequencies (number of buffers in which chan was bad):
                        [frq,allbad] = hist(tmp,unique(tmp));

                        % Mark bad based on threshold (currently ~5% of buffers (assuming 500 buffers)):
                        badchans = allbad(frq>0.05*Nbuf);
                        if isempty(badchans)
                            badstr = '';
                        else
                            badstr = sprintf(' -bad %s',num2str(badchans))
                        end


                        %% 2. tSSS and trans to first file (ie, align within subject if multiple runs)

                        tSSSstr = ' -st 10 -corr 0.98'; %'tSSSstr = '';

                        if mvcomp_fail(g,s,e,r) == 1
                            compstr = '';
                        else
                            compstr = sprintf(' -movecomp inter');
                        end

                        outfile = fullfile(sub_wd,sprintf('%s_trans1st',raw_stem))
                        if length(runs{e})>1 & r>1
                            transtr = sprintf(' -trans %s ',transfstfile)
                        else
                            transfstfile = [outfile '.fif'];
                            transtr = '';
                        end

                        dsstr = ' -ds 4';   % downsample to 250Hz

                        filestr = sprintf(' -f %s -o %s.fif',raw_file,outfile);
                        finstr = [maxfstr filestr basestr badstr tSSSstr compstr origstr transtr dsstr sprintf(' -v | tee %s.log',outfile)]
                        rik_eval(finstr);

                        rik_eval(sprintf('!echo ''Trans 1st...'' >> %s',movfile));
                        rik_eval(sprintf('!cat %s.log | sed -n ''/Position change/p'' | cut -f 7- -d '' '' >> %s',outfile,movfile));


                        %% 3. trans to default helmet space (align across subjects)

                        if TransDefaultFlag
                            transdeffile = fullfile(sub_wd,sprintf('%s_trans1stdef',raw_stem))
                            transtr = sprintf(' -trans default -origin %d %d %d -frame head -force',orig+[0 -13 6])
                            filestr = sprintf(' -f %s.fif -o %s.fif',outfile,transdeffile);
                            finstr = [maxfstr filestr transtr sprintf(' -v | tee %s.log',outfile)]
                            rik_eval(finstr);
                            rik_eval(sprintf('!echo ''Trans def...'' >> %s',movfile));
                            rik_eval(sprintf('!cat %s.log | sed -n ''/Position change/p'' | cut -f 7- -d '' '' >> %s',outfile,movfile));
                        end

                    end
                end
            end
        end

    else  %% non parfor...

        for s = 1:length(cbu_codes{g})
            raw_wd    = dir(fullfile(dat_wd,cbu_codes{g}{s},'1*'));  % Just to get date directory (assuming between 20*1*0 and 20*1*9Q!)
            raw_wd    = raw_wd.name;

            for e = 1:length(expts)

                transfstfile = ''; orig = []; rad=[]; fit=[]; transtr = {}; raw_stem = {};

                % Output directory
                sub_wd = fullfile(bas_wd,expts{e},groups{g},cbu_codes{g}{s}), cd(sub_wd)
                %try eval(sprintf('!mkdir %s',sub_wd)); end  % Try not allowed n parfor, so make directories in advance above

                rik_eval(sprintf('!touch %s',movfile));

                for r = 1:length(runs{e})

                    raw_file = dir(fullfile(dat_wd,cbu_codes{g}{s},raw_wd,sprintf('%s*%s_raw*',expts{e},runs{e}{r})));  % Get raw FIF file

                    if isempty(raw_file)
                        error('Could not find run %d for grp %s, sub %s, exp %s',r,groups{g},cbu_codes{g}{s},expts{e})
                    else
                        raw_stem = raw_file.name(1:(end-4));
                        raw_file = fullfile(dat_wd,cbu_codes{g}{s},raw_wd,raw_file.name);
                    end

                    if badrun(g,s,e,r)  % Care: need to note that this run not trans'ed (despite name), not sss, etc
                        outfile = fullfile(sub_wd,sprintf('%s_NoSSS',raw_stem));
                        filestr = sprintf(' -f %s -o %s.fif',raw_file,outfile);
                        finstr = [maxfstr filestr ' -nosss -ds 4' sprintf(' -v | tee %s.log',outfile)];
                        rik_eval(finstr);
                        %rik_eval(sprintf('!cp %s %s',raw_file,outfile));
                    else
                        %% Fit sphere (since better than MaxFilter does)
                        if r == 1  % fit sphere doesn't change with run!
                            incEEG = 0;
                            if exist(fullfile(sub_wd,'fittmp.txt')); delete(fullfile(sub_wd,'fittmp.txt')); end
                            if exist(fullfile(sub_wd,sprintf('run_%02d_hpi.txt',r))); delete(fullfile(sub_wd,sprintf('run_%02d_hpi.txt',r)));  end
                            [orig,rad,fit] = meg_fit_sphere(raw_file,sub_wd,sprintf('%s_hpi.txt',raw_stem),incEEG);
                            delete(fullfile(sub_wd,'fittmp.txt'));
                        end
                        origstr = sprintf(' -origin %d %d %d -frame head',orig(1),orig(2),orig(3))

                        badstr  = sprintf(' -autobad %d -badlimit %d',900,7); % 900s is 15mins - ie enough for whole recording!


                        %% 1. Bad channel detection (this email says important if doing tSSS later https://www.jiscmail.ac.uk/cgi-bin/webadmin?A2=NEUROMEG;d3f363f3.1205)

                        outfile = fullfile(sub_wd,sprintf('%s_bad',raw_stem));
                        filestr = sprintf(' -f %s -o %s.fif',raw_file,outfile);

                        % Write out movements too...
                        posfile = fullfile(sub_wd,sprintf('%s_headpos.txt',raw_stem));
                        compstr = sprintf(' -headpos -hpistep 10 -hp %s',posfile);

                        finstr = [maxfstr filestr origstr basestr badstr compstr sprintf(' -v | tee %s.log',outfile)]
                        if ParType==1
                            spmd; rik_eval(finstr); end
                        else
                            rik_eval(finstr);
                        end
                        delete(sprintf('%s.fif',outfile));

                        % Pull out bad channels from logfile:
                        badfile = sprintf('%s.txt',outfile); delete(badfile);
                        rik_eval(sprintf('!cat %s.log | sed -n -e ''/Detected/p'' -e ''/Static/p'' | cut -f 5- -d '' '' > %s',outfile,badfile));

                        tmp=dlmread(badfile,' '); Nbuf = size(tmp,1);
                        tmp=reshape(tmp,1,prod(size(tmp)));
                        tmp=tmp(tmp>0); % Omit zeros (padded by dlmread):

                        % Get frequencies (number of buffers in which chan was bad):
                        [frq,allbad] = hist(tmp,unique(tmp));

                        % Mark bad based on threshold (currently ~5% of buffers (assuming 500 buffers)):
                        badchans = allbad(frq>0.05*Nbuf);
                        if isempty(badchans)
                            badstr = '';
                        else
                            badstr = sprintf(' -bad %s',num2str(badchans))
                        end


                        %% 2. tSSS and trans to first file (ie, align within subject if multiple runs)

                        tSSSstr = ' -st 10 -corr 0.98'; %'tSSSstr = '';

                        if mvcomp_fail(g,s,e,r) == 1
                            compstr = '';
                        else
                            compstr = sprintf(' -movecomp inter');
                        end

                        outfile = fullfile(sub_wd,sprintf('%s_trans1st',raw_stem))
                        if length(runs{e})>1 & r>1
                            transtr = sprintf(' -trans %s ',transfstfile)
                        else
                            transfstfile = [outfile '.fif'];
                            transtr = '';
                        end

                        dsstr = ' -ds 4';   % downsample to 250Hz

                        filestr = sprintf(' -f %s -o %s.fif',raw_file,outfile);
                        finstr = [maxfstr filestr basestr badstr tSSSstr compstr origstr transtr dsstr sprintf(' -v | tee %s.log',outfile)]
                        if ParType==1
                            spmd; rik_eval(finstr); end
                        else
                            rik_eval(finstr);
                        end

                        rik_eval(sprintf('!echo ''Trans 1st...'' >> %s',movfile));
                        rik_eval(sprintf('!cat %s.log | sed -n ''/Position change/p'' | cut -f 7- -d '' '' >> %s',outfile,movfile));


                        %% 3. trans to default helmet space (align across subjects)

                        if TransDefaultFlag
                            transdeffile = fullfile(sub_wd,sprintf('%s_trans1stdef',raw_stem))
                            transtr = sprintf(' -trans default -origin %d %d %d -frame head -force',orig+[0 -13 6])
                            filestr = sprintf(' -f %s.fif -o %s.fif',outfile,transdeffile);
                            finstr = [maxfstr filestr transtr sprintf(' -v | tee %s.log',outfile)]
                            if ParType==1
                                spmd; rik_eval(finstr); end
                            else
                                rik_eval(finstr);
                            end
                            rik_eval(sprintf('!echo ''Trans def...'' >> %s',movfile));
                            rik_eval(sprintf('!cat %s.log | sed -n ''/Position change/p'' | cut -f 7- -d '' '' >> %s',outfile,movfile));
                        end
                    end
                end
            end
        end
    end
end

if ParType
    matlabpool close force CBU_Cluster
end

CbuMeg: Maxfilter_V2.2 (last edited 2023-03-03 10:53:29 by OlafHauk)