1 function varargout = process_beamformer_test( varargin )
   2 % PROCESS_BEAMFORMER_TEST: 
   3 
   4 % @=============================================================================
   5 % This software is part of the Brainstorm software:
   6 % http://neuroimage.usc.edu/brainstorm
   7 % 
   8 % Copyright (c)2000-2013 Brainstorm by the University of Southern California
   9 % This software is distributed under the terms of the GNU General Public License
  10 % as published by the Free Software Foundation. Further details on the GPL
  11 % license can be found at http://www.gnu.org/copyleft/gpl.html.
  12 % 
  13 % FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
  14 % UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
  15 % WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
  16 % MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
  17 % LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
  18 %
  19 % For more information type "brainstorm license" at command prompt.
  20 % =============================================================================@
  21 %
  22 % Authors: 
  23 
  24 eval(macro_method);
  25 end
  26 
  27 
  28 %% ===== GET DESCRIPTION =====
  29 function sProcess = GetDescription() %#ok<DEFNU>
  30     % Description the process
  31     sProcess.Comment     = 'Beamformer test #1';
  32     sProcess.FileTag     = '';
  33     sProcess.Category    = 'Custom';
  34     sProcess.SubGroup    = 'Sources';
  35     sProcess.Index       = 1000;
  36     % Definition of the input accepted by this process
  37     sProcess.InputTypes  = {'data'};
  38     sProcess.OutputTypes = {'results'};
  39     sProcess.nInputs     = 1;
  40     sProcess.nMinFiles   = 1;
  41     % Definition of the options
  42     % === ACTIVE STATE
  43     sProcess.options.active_time.Comment = 'Active state time window: ';
  44     sProcess.options.active_time.Type    = 'timewindow';
  45     sProcess.options.active_time.Value   = [];
  46     % === ACTIVE STATE
  47     sProcess.options.minvar_time.Comment = 'Minimum variance time window: ';
  48     sProcess.options.minvar_time.Type    = 'timewindow';
  49     sProcess.options.minvar_time.Value   = [];
  50     % === REGULARIZATION
  51     sProcess.options.reg.Comment = 'Regularization parameter: ';
  52     sProcess.options.reg.Type    = 'value';
  53     sProcess.options.reg.Value   = {0.03, ' ', 4};
  54     % === Sensor types
  55     sProcess.options.sensortypes.Comment = 'Sensor types or names (empty=all): ';
  56     sProcess.options.sensortypes.Type    = 'text';
  57     sProcess.options.sensortypes.Value   = 'MEG, EEG';
  58 end
  59 
  60 
  61 %% ===== FORMAT COMMENT =====
  62 function Comment = FormatComment(sProcess) %#ok<DEFNU>
  63     Comment = sProcess.Comment;
  64 end
  65 
  66 
  67 %% ===== RUN =====
  68 function OutputFiles = Run(sProcess, sInputs) %#ok<DEFNU>
  69     % Initialize returned list of files
  70     OutputFiles = {};
  71     % Get option values
  72     ActiveTime  = sProcess.options.active_time.Value{1};
  73     MinVarTime  = sProcess.options.minvar_time.Value{1};
  74     Reg         = sProcess.options.reg.Value{1};
  75     SensorTypes = sProcess.options.sensortypes.Value;
  76     
  77     % ===== LOAD CHANNEL FILE =====
  78     % Load channel file
  79     ChannelMat = in_bst_channel(sInputs(1).ChannelFile);
  80     % Find the MEG channels
  81 %     iMEG = good_channel(ChannelMat.Channel, [], 'MEG');
  82 %     iEEG = good_channel(ChannelMat.Channel, [], 'EEG');
  83 %     iSEEG = good_channel(ChannelMat.Channel, [], 'SEEG');
  84 %     iECOG = good_channel(ChannelMat.Channel, [], 'ECOG');
  85     iChannels = channel_find(ChannelMat.Channel, SensorTypes);
  86     
  87     % ===== LOAD HEAD MODEL =====
  88     % Get channel study
  89     [sChannelStudy, iChannelStudy] = bst_get('ChannelFile', sInputs(1).ChannelFile);
  90     % Load the default head model
  91     HeadModelFile = sChannelStudy.HeadModel(sChannelStudy.iHeadModel).FileName;
  92     sHeadModel = load(file_fullpath(HeadModelFile));
  93     % Get number of sources
  94     nSources = length(sHeadModel.GridLoc);
  95     
  96     % ===== LOAD THE DATA =====
  97     % Read the first file in the list, to initialize the loop
  98     DataMat = in_bst(sInputs(1).FileName, [], 0);
  99     nChannels = size(DataMat.F,1);
 100     nTime     = size(DataMat.F,2);
 101     Time = DataMat.Time;
 102     % Find the indices for covariance calculation
 103     iActiveTime = panel_time('GetTimeIndices', Time, ActiveTime);
 104     iMinVarTime = panel_time('GetTimeIndices', Time, MinVarTime);
 105     % Initialize the covariance matrices
 106     ActiveCov = zeros(nChannels, nChannels);
 107     MinVarCov = zeros(nChannels, nChannels);
 108     nTotalActive = zeros(nChannels, nChannels);
 109     nTotalMinVar = zeros(nChannels, nChannels);
 110     % Reading all the input files in a big matrix
 111     for i = 1:length(sInputs)
 112         % Read the file #i
 113         DataMat = in_bst(sInputs(i).FileName, [], 0);
 114         % Check the dimensions of the recordings matrix in this file
 115         if (size(DataMat.F,1) ~= nChannels) || (size(DataMat.F,2) ~= nTime)
 116             % Add an error message to the report
 117             bst_report('Error', sProcess, sInputs, 'One file has a different number of channels or a different number of time samples.');
 118             % Stop the process
 119             return;
 120         end
 121         % Get good channels
 122         iGoodChan = find(DataMat.ChannelFlag == 1);
 123         % Average baseline values
 124         FavgActive = mean(DataMat.F(iGoodChan,iActiveTime), 2);
 125         FavgMinVar = mean(DataMat.F(iGoodChan,iMinVarTime), 2);
 126         % Remove average
 127         DataActive = bst_bsxfun(@minus, DataMat.F(iGoodChan,iActiveTime), FavgActive);
 128         DataMinVar = bst_bsxfun(@minus, DataMat.F(iGoodChan,iMinVarTime), FavgMinVar);
 129         % Compute covariance for this file
 130         fileActiveCov = DataMat.nAvg .* (DataActive * DataActive');
 131         fileMinVarCov = DataMat.nAvg .* (DataMinVar * DataMinVar');
 132         % Add file covariance to accumulator
 133         ActiveCov(iGoodChan,iGoodChan) = ActiveCov(iGoodChan,iGoodChan) + fileActiveCov;
 134         MinVarCov(iGoodChan,iGoodChan) = MinVarCov(iGoodChan,iGoodChan) + fileMinVarCov;
 135         nTotalActive(iGoodChan,iGoodChan) = nTotalActive(iGoodChan,iGoodChan) + length(iActiveTime);
 136         nTotalMinVar(iGoodChan,iGoodChan) = nTotalMinVar(iGoodChan,iGoodChan) + length(iMinVarTime);
 137     end
 138     % Remove zeros from N matrix
 139     nTotalActive(nTotalActive <= 1) = 2;
 140     nTotalMinVar(nTotalMinVar <= 1) = 2;
 141     % Divide final matrix by number of samples
 142     ActiveCov = ActiveCov ./ (nTotalActive - 1);
 143     MinVarCov = MinVarCov ./ (nTotalMinVar - 1);
 144 
 145     
 146     % ===== PROCESS =====
 147     % Processing iChannels
 148     %%%% TO EDIT %%%%%
 149     ImageGridAmp = rand(nSources, nTime);
 150     GridLoc = rand(nSources, 3);
 151     
 152     % ===== SAVE THE RESULTS =====
 153     % Create a new data file structure
 154     ResultsMat = db_template('resultsmat');
 155     ResultsMat.ImagingKernel = [];
 156     ResultsMat.ImageGridAmp  = ImageGridAmp;
 157     ResultsMat.nComponents   = 1;   % 1 or 3
 158     ResultsMat.Comment       = 'TEST';
 159     ResultsMat.Function      = 'NameMethod';
 160     ResultsMat.Time          = Time;           % Leave it empty if using ImagingKernel
 161     ResultsMat.DataFile      = [];
 162     ResultsMat.HeadModelFile = HeadModelFile;
 163     ResultsMat.HeadModelType = sHeadModel.HeadModelType;
 164     ResultsMat.ChannelFlag   = [];
 165     ResultsMat.GoodChannel   = iChannels;
 166     ResultsMat.SurfaceFile   = sHeadModel.SurfaceFile;
 167     ResultsMat.GridLoc       = GridLoc;
 168 
 169     % === NOT SHARED ===
 170     % Get the output study (pick the one from the first file)
 171     iStudy = sInputs(1).iStudy;
 172     % Create a default output filename 
 173     OutputFiles{1} = bst_process('GetNewFilename', fileparts(sInputs(1).FileName), 'results_');
 174 %     % === SHARED ===
 175 %     % Get the output study (pick the one from the first file)
 176 %     iStudy = iChannelStudy;
 177 %     % Create a default output filename 
 178 %     OutputFiles{1} = bst_process('GetNewFilename', fileparts(sInputs(1).ChannelFile), 'results_KERNEL_');
 179 
 180 %     OutputFiles{2} = ...
 181 
 182     % Save on disk
 183     save(OutputFiles{1}, '-struct', 'ResultsMat');
 184     % Register in database
 185     db_add_data(iStudy, OutputFiles{1}, ResultsMat);
 186 end