[Master Index] [Index for Toolbox]

# lcmvbf

## (Toolbox/lcmvbf.m in BrainStorm 2.0 (Alpha))

### Function Synopsis

    Results = lcmvbf(OPTIONS);

### Help Text

    LCMVBF - Spatial filtering solution as called from command line or BST_SOURCEIMAGING.M
     function Results = lcmvbf(OPTIONS);
     -------------------------------------------------------------------------------------------------
     NOTICE:
     This function is not optimized for stand-alone command calls.
     Please use the generic BST_SOURCEIMAGING function to access Min.Norm. source map computations.
     -------------------------------------------------------------------------------------------------

     INPUTS
           OPTIONS - a structure of parameters
     OUTPUTS
           Results - a regular BrainStorm Results structure (see documentation for details)

     Details of the OPTIONS structure

       Data Information:
              .DataFile   : A string actual data file requested for localization. If isempty, look for data array(s) in .Data
              .Data       : Array of nChannel x nTime array(s) of data.
              .SegmentIndices: Vector of all time samples (in sample values)
                               of the latencies specifying the time window for analysis.
              .Trial      : A scalar specifying the trial from which to extract the data segment
                            (not empty only when data in .DataFile is stored in native format)
              .ChannelFlag: a vector of flags (1 for good, -1 for bad) of channels in the data file to process
              .DataTypes  : A scalar or a vector of codes to specidy which is the type of data to use:
                            1 is MEG, 2 is EEG, 3 is Fusion
                            (Default is all data type available)
              .FusionType : Code for the type of fusion to apply between MEG and EEG data
                            1 is linear, 2 is non-linear
                            (Defaut is 1)
              .Data       : Specify some data to be passed independently of any data file
                                .F : Cell array of nChannel x nTime array(s) of data.
                                     If not empty: uses this array as input for inverse process.
                                .ChannelFile : Channel file corresponding to data in .F
                                .Time : Vector specifying the time instants of data samples in .Data.F
                            (Default for .Data is empty)

       Source Support Information:
              .GridLoc    : a string or a 3xN matrix
                            - if a string : specifies the name of the BrainStorm tessellation file where cortical support is stored
                            - if a matrix : contains the coordinates of the cortical locations where to compute the cortical current density
              .iGrid      : when .GridLoc is a file name, specify which imaging grid to use (Default is 1).

       Forward Field Information :
              .HeadModelFile: Cell array of strings for filename of the BrainStorm headmodel file(s) to use for imaging. If is empty, look for.GainFile field.
                              There are only 2 options at this point:
                                    - .HeadModelFile is of length 1, a single headmodel file is used for all data sets passed to the bst_sourceimaging.
                                    - .HeadModelFile is of length the number of data sets passed to bst_sourceimaging (i.e. the length of .DataFile).
              .GainFile   : Cell array of filenames of image gain matrix (.bin) for the selected cortical support. If is empty, look for .Gain field.
                            Remarks about the length of .HeadModelFile hold for .GainFile too.
              .Gain       : Cell array containing the gain matrix of the cortical support (either a single cell or a cell array the length the number of data sets passed to function)
                            Remarks about the length of .HeadModelFile hold for .GainFile too.
                            (Default is empty).

       Processing Parameters :
              .Method     : A string that specifies the imaging method
                            'Minimum-Norm Imaging': MN estimate of current density using MN priors on source amplitudes (Default)
                            'Recursive Min. Norm.': Recursive MN estimate based on adaptation of FOCUSS (NOT AVAILABLE)
                            'SPTF' : Scanning technique constrained on cortical surface (see Denis Schwartz) (UNDER DEVELOPMENT)
                            'ST-MAP' : Non-linear imaging technique based on sparse-focal image models (NOT AVAILABLE)
                            'RAP-MUSIC' : link to Mosher's GUI (UNDER DEVELOPMENT)
                            (Default is 'Minimum-Norm Imaging')
              .Tikhonov   : Hyperparameter used in Tykhonov regularization
                            (Default is 10)
              .Rank       : Specify a rank for signal subspace following SVD of the spatio-temporal data matrix
                            (Default is full rank)
              .ComputeKernel : If 1, compute MN kernel to be applied subsequently to data instead of full ImageGridAmp array
                             (Default is 0: compute full ImageGridAmp array)

       Output Information:
              .ResultFile : optional string for filename where to write results
                            if 'default', filename as used in GUI environment is used
                            (Default is '' (no Result file created))
              .Verbose    : Turn ON/OFF verbose mode
                            (Default is 1 (ON))
              .GUIHandles : A vector of GUI handles to the SOURCEIMAGING GUI
                            when operating in graphic mode: passed automatically to bst_sourceimaging by the source imaging switchyard
                            in command line mode: leave this field empty (default).
              .NNormalization : Apply noise normalization
                            (Default is 0)

### Cross-Reference Information

This function calls

    bst_message_window                 C:\BrainStorm_2001\Toolbox\bst_message_window.m
    colnorm                            C:\BrainStorm_2001\Toolbox\colnorm.m
    good_channel                       C:\BrainStorm_2001\Toolbox\good_channel.m
    inorcol                            C:\BrainStorm_2001\Toolbox\inorcol.m
    meg4read                           C:\BrainStorm_2001\Toolbox\meg4read.m
    norcol                             C:\BrainStorm_2001\Toolbox\norcol.m
    read_gain                          C:\BrainStorm_2001\Toolbox\read_gain.m
    regcheck                           C:\BrainStorm_2001\Toolbox\regcheck.m

This function is called by

    bst_sourceimaging                  C:\BrainStorm_2001\Toolbox\bst_sourceimaging.m
    dataplot_cb                        C:\BrainStorm_2001\Toolbox\dataplot_cb.m

### Listing of function C:\BrainStorm_2001\Toolbox\lcmvbf.m

    function Results = lcmvbf(OPTIONS);
    %LCMVBF - Spatial filtering solution as called from command line or BST_SOURCEIMAGING.M
    % function Results = lcmvbf(OPTIONS);
    % -------------------------------------------------------------------------------------------------
    % NOTICE:
    % This function is not optimized for stand-alone command calls.
    % Please use the generic BST_SOURCEIMAGING function to access Min.Norm. source map computations.
    % -------------------------------------------------------------------------------------------------
    %
    %
    % INPUTS
    %       OPTIONS - a structure of parameters
    % OUTPUTS
    %       Results - a regular BrainStorm Results structure (see documentation for details)
    %
    % Details of the OPTIONS structure
    %
    %   Data Information:
    %          .DataFile   : A string actual data file requested for localization. If isempty, look for data array(s) in .Data
    %          .Data       : Array of nChannel x nTime array(s) of data.
    %          .SegmentIndices: Vector of all time samples (in sample values)
    %                           of the latencies specifying the time window for analysis.
    %          .Trial      : A scalar specifying the trial from which to extract the data segment
    %                        (not empty only when data in .DataFile is stored in native format)
    %          .ChannelFlag: a vector of flags (1 for good, -1 for bad) of channels in the data file to process
    %          .DataTypes  : A scalar or a vector of codes to specidy which is the type of data to use:
    %                        1 is MEG, 2 is EEG, 3 is Fusion
    %                        (Default is all data type available)
    %          .FusionType : Code for the type of fusion to apply between MEG and EEG data
    %                        1 is linear, 2 is non-linear
    %                        (Defaut is 1)
    %          .Data       : Specify some data to be passed independently of any data file
    %                            .F : Cell array of nChannel x nTime array(s) of data.
    %                                 If not empty: uses this array as input for inverse process.
    %                            .ChannelFile : Channel file corresponding to data in .F
    %                            .Time : Vector specifying the time instants of data samples in .Data.F
    %                        (Default for .Data is empty)
    %
    %   Source Support Information:
    %          .GridLoc    : a string or a 3xN matrix
    %                        - if a string : specifies the name of the BrainStorm tessellation file where cortical support is stored
    %                        - if a matrix : contains the coordinates of the cortical locations where to compute the cortical current density
    %          .iGrid      : when .GridLoc is a file name, specify which imaging grid to use (Default is 1).
    %
    %   Forward Field Information :
    %          .HeadModelFile: Cell array of strings for filename of the BrainStorm headmodel file(s) to use for imaging. If is empty, look for.GainFile field.
    %                          There are only 2 options at this point:
    %                                - .HeadModelFile is of length 1, a single headmodel file is used for all data sets passed to the bst_sourceimaging.
    %                                - .HeadModelFile is of length the number of data sets passed to bst_sourceimaging (i.e. the length of .DataFile).
    %          .GainFile   : Cell array of filenames of image gain matrix (.bin) for the selected cortical support. If is empty, look for .Gain field.
    %                        Remarks about the length of .HeadModelFile hold for .GainFile too.
    %          .Gain       : Cell array containing the gain matrix of the cortical support (either a single cell or a cell array the length the number of data sets passed to function)
    %                        Remarks about the length of .HeadModelFile hold for .GainFile too.
    %                        (Default is empty).
    %
    %   Processing Parameters :
    %          .Method     : A string that specifies the imaging method
    %                        'Minimum-Norm Imaging': MN estimate of current density using MN priors on source amplitudes (Default)
    %                        'Recursive Min. Norm.': Recursive MN estimate based on adaptation of FOCUSS (NOT AVAILABLE)
    %                        'SPTF' : Scanning technique constrained on cortical surface (see Denis Schwartz) (UNDER DEVELOPMENT)
    %                        'ST-MAP' : Non-linear imaging technique based on sparse-focal image models (NOT AVAILABLE)
    %                        'RAP-MUSIC' : link to Mosher's GUI (UNDER DEVELOPMENT)
    %                        (Default is 'Minimum-Norm Imaging')
    %          .Tikhonov   : Hyperparameter used in Tykhonov regularization
    %                        (Default is 10)
    %          .Rank       : Specify a rank for signal subspace following SVD of the spatio-temporal data matrix
    %                        (Default is full rank)
    %          .ComputeKernel : If 1, compute MN kernel to be applied subsequently to data instead of full ImageGridAmp array
    %                         (Default is 0: compute full ImageGridAmp array)
    %
    %   Output Information:
    %          .ResultFile : optional string for filename where to write results
    %                        if 'default', filename as used in GUI environment is used
    %                        (Default is '' (no Result file created))
    %          .Verbose    : Turn ON/OFF verbose mode
    %                        (Default is 1 (ON))
    %          .GUIHandles : A vector of GUI handles to the SOURCEIMAGING GUI
    %                        when operating in graphic mode: passed automatically to bst_sourceimaging by the source imaging switchyard
    %                        in command line mode: leave this field empty (default).
    %          .NNormalization : Apply noise normalization
    %                        (Default is 0)

    %<autobegin> ---------------------- 12-Oct-2004 12:01:16 -----------------------
    % --------- Automatically Generated Comments Block Using AUTO_COMMENTS ---------
    %
    % CATEGORY: Unknown Category
    %
    % Alphabetical list of external functions (non-Matlab):
    %   toolbox\bst_message_window.m
    %   toolbox\colnorm.m
    %   toolbox\good_channel.m
    %   toolbox\inorcol.m
    %   toolbox\meg4read.m
    %   toolbox\norcol.m
    %   toolbox\read_gain.m
    %   toolbox\regcheck.m
    %
    % Subfunctions in this file, in order of occurrence in file:
    %   varargout = gainmat_covar(varargin)
    %
    % At Check-in: $Author: Mosher $  $Revision: 3 $  $Date: 10/12/04 10:23a $
    %
    % This software is part of BrainStorm Toolbox Version 2.0 (Alpha) 12-Oct-2004
    %
    % Principal Investigators and Developers:
    % ** Richard M. Leahy, PhD, Signal & Image Processing Institute,
    %    University of Southern California, Los Angeles, CA
    % ** John C. Mosher, PhD, Biophysics Group,
    %    Los Alamos National Laboratory, Los Alamos, NM
    % ** Sylvain Baillet, PhD, Cognitive Neuroscience & Brain Imaging Laboratory,
    %    CNRS, Hopital de la Salpetriere, Paris, France
    %
    % See BrainStorm website at http://neuroimage.usc.edu for further information.
    %
    % Copyright (c) 2004 BrainStorm by the University of Southern California
    % This software distributed  under the terms of the GNU General Public License
    % as published by the Free Software Foundation. Further details on the GPL
    % license can be found at http://www.gnu.org/copyleft/gpl.html .
    %
    % FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
    % UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
    % WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
    % MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
    % LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
    %<autoend> ------------------------ 12-Oct-2004 12:01:16 -----------------------

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % /---Script Author--------------------------------------\
    % | *** Esen Kucukaltun-Yildirim                           |
    % |  USC-SIPI                                              |
    % \------------------------------------------------------/

    % Date of creation: October 2004 (adapted from minnorm.m)
    %
    % ----------------------------- Script History ---------------------------------
    % EK  7-Oct-2004  New outputput formats 'Source Power' and Noise Normalized
    %                 Filter Outputs are added
    % JCM 12-Oct-2004 Header comments
    % ----------------------------- Script History ---------------------------------

    % General hard-wired parameter values:--------------------------------------------------------------
    DataType = {'MEG','EEG'}; % String labels corresponding to OPTIONS.DataTypes flags
    OPTIONS.BLOCK = 10000; % process by blocks
    OPTIONS.ExpoiG = .4; % Exponential weighting for column normalization
    %---------------------------------------------------------------------------------------------------

    % Load Headmodel information -----------------------------------------------------------------------

    % From file, if passed to function
    if ~isempty(OPTIONS.HeadModelFile)
        HeadModel = load(OPTIONS.HeadModelFile,'Function','Param',...
            'GridLoc','Gain'); % Load Basics

        if ~isempty(OPTIONS.GridLoc)
            HeadModel.GridLoc = {OPTIONS.GridLoc};
        end
        if ~isempty(OPTIONS.Gain)
            HeadModel.Gain = OPTIONS.Gain;
        end

        HeadModelFileWhos = whos('-file',OPTIONS.HeadModelFile);

        clear tmp

        % Gain matrix full filename
        gainfile = fullfile(fileparts(OPTIONS.HeadModelFile),HeadModel.Gain{1});

    end

    % Load Data Information -----------------------------------------------------------------------

    if(~isempty(findstr('_results',deblank(lower(OPTIONS.DataFile))))),
        % CHEAT
        % The DataName has the string '_results' in it
        % We are trying to treat independent topographies as min norm data. Eventually this
        %  capability will be moved to a new routine.
        % name has the string '_results' in it
        Results = load(OPTIONS.DataFile);
        % now synthesize a Data structure from this
        [ignorm,Anrm] = colnorm(Results.IndepTopo);

        Data = struct('F',Anrm,'ChannelFlag',Results.ChannelFlag,...
            'Time',[1:size(Results.IndepTopo,2)],... % time is now simply an indexer
            'NoiseCov',Results.NoiseCov,'SourceCov',Results.SourceCov,'Projector',Results.Projector,...
            'Comment',sprintf('Synthesized Data Set from Results: %s',Results.Comment));
    else
        if isempty(OPTIONS.Data)
            % assume it is a standard data file
            DataArray = 0; % Correponsing flag
            Data = load(OPTIONS.DataFile,'F','ChannelFlag','Time','Projector','Comment');
        else
            DataArray = 1;
            Data = OPTIONS.Data;
            Data.Projector = [];
            Data.Comment = '';
            Data.ChannelFlag = ones(size(Data.F,1),1);
        end

    end

    % Data & Channel information ------------------------------------------------------------------
    if isfield(OPTIONS,'ChannelFile')
        Channel = load(OPTIONS.ChannelFile);
        OPTIONS.Channel = Channel.Channel; clear Channel
    end

    % now alter the data according to the bad channels
    OPTIONS.GoodChannel = good_channel(OPTIONS.Channel,OPTIONS.ChannelFlag,DataType{OPTIONS.DataType}); %good channels for selected modality
    OPTIONS.EEGChannel = good_channel(OPTIONS.Channel,ones(size(OPTIONS.Channel)),'EEG');
    OPTIONS.MEGChannel = good_channel(OPTIONS.Channel,ones(size(OPTIONS.Channel)),'MEG');

    if ischar(Data.F) % Data is in native file format
        VERBOSE = 1;
        Data.F = meg4read(Data.F,OPTIONS.Trial,OPTIONS.SegmentIndices,OPTIONS.Verbose);
        Data.FBaseline = meg4read(Data.F,OPTIONS.Trial,OPTIONS.NoiseSegmentIndices,OPTIONS.Verbose);
    else
        % Trim the data to the requested segments
        if OPTIONS.NNormalization
            Data.FBaseline = Data.F(:,OPTIONS.NoiseSegmentIndices);
        end
        Data.F = Data.F(:,OPTIONS.SegmentIndices);
    end

    try
        Data.F = Data.F(OPTIONS.GoodChannel,:); % good channels only
        if OPTIONS.NNormalization
            Data.FBaseline = Data.FBaseline(OPTIONS.GoodChannel,:); % good channels for noise segment too
        end

    catch
        errordlg(sprintf('Data array has %d rows while there are %d %s channels. Possible problem is wrong type assignement of reference vs. measurement channels in channel file %s',...
            size(Data.F,1),length(OPTIONS.GoodChannel), DataType{OPTIONS.DataType}, OPTIONS.ChannelFile ),'Inconsistent channel and data array sizes')
        return
    end

    Time = Data.Time(OPTIONS.SegmentIndices);

    % Is there any projector for the data ?
    if(~isempty(Data.Projector)),
        Data.Projector = Data.Projector(OPTIONS.GoodChannel,:);
    end

    % Get HeadModel parmaters for GoodChannels
    HeadModel.Param = HeadModel.Param(OPTIONS.GoodChannel);

    % Update MEG and EEG channel indices by removing bad channels.
    OPTIONS.EEGChannel = intersect(OPTIONS.GoodChannel,OPTIONS.EEGChannel);
    OPTIONS.MEGChannel = intersect(OPTIONS.GoodChannel,OPTIONS.MEGChannel);

    if OPTIONS.DataType == 2 & ...
            length(OPTIONS.GoodChannel) < length(OPTIONS.EEGChannel) % EEG data with bad channels

        % Adjust average reference contribution from covariance matrix
        if isempty(find(mean(Data.F >eps)))
            % Average Reference (don't use the channel.comment field that may not have been altered following manipulations in the data viewer tool)
            % Load gain matrix
            G = read_gain(gainfile,[]);
            if strcmp(OPTIONS.Col_norm,'y');
                G = G*inorcol(G);
            end

            Gm = mean(G,1)';

            I = ones(length(OPTIONS.GoodChannel),1);
            % Reopen the Gain file that was closed during call to read_gain
            fid = fopen(...
                fullfile(fileparts(OPTIONS.HeadModelFile),HeadModel.Gain{1}),...
                'r','ieee-be'); % HeadModel.Gain{1} - CHEAT Dipoles only
        else
            % Do nothing
        end
    end

    OPTIONS = regcheck(OPTIONS); % returns OPTIONS.REG set to the appropriate field

    % Did the user provide an existing control?
    if(isfield(Data,'Projector')),
        %Decompose it:
        if(isempty(Data.Projector)),
            A = [];
            Ua = [];
        else
            A = Data.Projector; % initialize
            Ua = orth(A);
        end
    else % user did not give a projector
        A = [];
        Ua = [];
    end

    if isempty(OPTIONS.Rank)
        OPTIONS.Rank = min(size(Data.F));
    end

    if(OPTIONS.Rank < min(size(Data.F))), % reduced rank

        % User has asked for an rank RANK subspace search.
        % Form the signal subspace

        if(size(Data.F,1) >= size(Data.F,2)), % skinny matrix
            [Uf,Sf,Vf] = svd(Data.F,0); % economy
        else % wide matrix
            [Vf,Sf,Uf] = svd(Data.F',0);
        end

        % form the reduced rank data
        Data.F = Uf(:,1:OPTIONS.Rank)*(Sf(1:OPTIONS.Rank,1:OPTIONS.Rank)*Vf(:,1:OPTIONS.Rank)');

        % project away if needed
        if(~isempty(Ua)),
            Data.F = Data.F - Ua*(Ua'*Data.F);
        end

        clear Ua Sa Va Uf Sf Vf
    end

    %                             --------------------------------------------------------------------------------------

    % Begin LCMV BF estimate --------------------------------------------------------------------------------------

    %                             --------------------------------------------------------------------------------------

    Kernel = read_gain(gainfile,OPTIONS.GoodChannel,[]);
    nSources = size(Kernel,2);
    nChannels = size(Kernel,1);

    Cm = cov(Data.F');
    %invert data covariance
    [U, S, V] = svd(Cm);
    % if user select not to regularize .Tikhonov is 0 and it is a direct
    % inverse
    S = diag(S);
    Lambda = S(1)*OPTIONS.Tikhonov/100 ;
    Si = diag(1../(S+Lambda)); % filtered inverse
    Cm_inv = V*Si*U';   clear U S Si V
    %design beamformer
    source_power_inv = sum( Kernel.*(Cm_inv*Kernel) );
    if get(OPTIONS.GUIHandles.sourcePower, 'Value') ==1 %Calculate source power
        bst_message_window('Calculating the Source Power. . .');

        ImageGridAmp = (1 ./ source_power_inv)' ;
        ImagingKernel = [];
        Fsynth = [];
    elseif get(OPTIONS.GUIHandles.normalizedOut, 'Value') ==1
        bst_message_window('Calculating the spatial filter. . .');
        spatialFilter = Kernel'*Cm_inv;
        spatialFilter= (spatialFilter) ./ repmat(source_power_inv', 1,size(spatialFilter,2));
        ImageGridAmp=[];
        Fsynth = [];

        tmp = (spatialFilter*Data.FBaseline);
        mean_Baseline= mean(tmp, 2);
        std_Baseline = std(tmp, 0, 2);
        tmp = (spatialFilter);
        %     ImagingKernel = (tmp - repmat(mean_Baseline, 1, size(tmp, 2))) ./  repmat(std_Baseline, 1, size(tmp, 2));
        ImagingKernel = tmp  ./  repmat(std_Baseline, 1, size(tmp, 2));
        ImagingKernel = single(ImagingKernel);

    elseif OPTIONS.NNormalization  %Calculate the neural index
        bst_message_window('Calculating the Neural activity index. . .');
        Cn = cov(Data.FBaseline');
        [U, S, V] = svd(Cn);S = diag(S);
        Lambda = S(1)*OPTIONS.Tikhonov/100 ;
        Si = diag(1../(S+Lambda)); % filtered inverse
        Cn_inv = V*Si*U';   clear U S Si V

        noise_power_inv = sum( Kernel.*(Cn_inv*Kernel) );
        ImageGridAmp = (noise_power_inv ./ source_power_inv)' ;
        ImagingKernel = [];
        Fsynth = [];
    elseif OPTIONS.NNormalization==0 %Calculate the current density kernel multiply by the data Data.F to get the reconstruction map
        bst_message_window('Calculating the spatial filter. . .');
        spatialFilter = Kernel'*Cm_inv;
        ImagingKernel= (spatialFilter) ./ repmat(source_power_inv', 1,size(spatialFilter,2));
        ImageGridAmp=[]; clear spatialFilter
        Fsynth = [];
    end

    if OPTIONS.Verbose
        bst_message_window({...
            'Inverse Transformation -> DONE.'...
            })
    end

    %
    % beamformer estimate - DONE --------------------------------------------------------------------------------------------------------------------------

    % Save Results in Result file -----------------------------------------------------------------------------------------------------------------
    %
    % Fields from Results structure
    Channel = OPTIONS.Channel;
    ChannelFlag = Data.ChannelFlag;
    Comment = sprintf('MIN NORM, rank %.0f',OPTIONS.Rank);
    DataFile = OPTIONS.DataFile;
    Date = datestr(datenum(now),0);
    Time = [OPTIONS.SegmentIndices];
    ImageGridTime = Data.Time(Time);
    ImageGridAmp = ImageGridAmp;  % monster matrix when ImagingKernel is not computed
    Function = mfilename; % name of this calling routine
    HeadModelFile = OPTIONS.HeadModelFile;
    ModelGain = HeadModel.Gain{1};
    Projector = Data.Projector;
    SourceLoc = HeadModel.GridLoc{1};
    SourceOrder = -1;

    NoiseCov = [];
    SourceCov = [];
    PatchNdx = [];
    PatchAmp = [];

    % Reduced Results structure for call from bst_sourceimaging.
    Results = struct('Comment',Comment,'Function',mfilename,...
        'OPTIONS',OPTIONS,'Date','','Time',Time,'HeadModelFile',HeadModelFile,...
        'Channel',Channel,'ChannelFlag',Data.ChannelFlag,'DataFile',OPTIONS.DataFile,...
        'NoiseCov',NoiseCov,'SourceCov',SourceCov,...
        'Projector',Projector,'SourceLoc',SourceLoc,'SourceOrder',SourceOrder,...
        'SourceOrientation',[],'TimeSeries',[],'ImagingKernel',ImagingKernel,'ModelGain',ModelGain,...
        'PatchNdx',PatchNdx,'PatchAmp',PatchAmp, ...
        'ImageGridAmp',ImageGridAmp,'ImageGridTime',ImageGridTime,'Fsynth',Fsynth);

    return

    % -----------------------------------------------------------------------------------------------------
    %
    %
    %                                             SubFunctions
    %
    %
    % -----------------------------------------------------------------------------------------------------

    function varargout = gainmat_covar(varargin)

    % GAINMAT_COVAR Computation of gain matrix covariance
    % function [GainCovar] = gainmat_covar(gainfile,OPTIONS)
    % Computation of gain matrix covariance
    %
    % INPUTS
    %
    % gainfile : a string - filename of the binary gain matrix file (imaging grid)
    % OPTIONS  : a structure specifying the options for the computation
    %            .Verbose : turn on/off verbose mode (1/0)
    %            .DataType: scalar specifying the type of data
    %                       1 is MEG
    %                       2 is EEG
    %                       3 is MEG/EEG fusion
    %            .FFNormalization : flag specifying whether forward field normalization is requested
    %                       0 off
    %                       1 on
    %            .GoodChannel : Indices of good channels for current DataType
    %            .MEGChannel  : Indices of all MEG channels
    %            .EEGChannel  : Indcies of all EEG channels
    %
    % OUTPUTS
    % GainCovar : a structure with the proper gain covar matrix outputs
    %
    % /---Script Author--------------------------------------\
    % | *** John C. Mosher, Ph.D.                            |
    % |  Biophysics Group                                    |
    % |                                                      |
    % | *** Sylvain Baillet, Ph.D.                           |
    % | Cognitive Neuroscience & Brain Imaging Laboratory    |
    % | CNRS UPR640 - LENA                                   |
    % | Hopital de la Salpetriere, Paris, France             |
    % | sylvain.baillet@chups.jussieu.fr                     |
    % \------------------------------------------------------/
    %
    % Date of creation: October 2002
    % Script History -----------------------------------------------------------------------------------
    %---------------------------------------------------------------------------------------------------

    % Map entries
    gainfile = varargin{1};
    OPTIONS = varargin{2};

    % Read number of sources and total number of channels
    tmp = read_gain(gainfile,OPTIONS.GoodChannel(1),[]);
    nsrc = length(tmp); clear tmp % Number of sources
    tmp = read_gain(gainfile,[],1);
    nchan = length(tmp); clear tmp % Number of channels

    % if OPTIONS.Verbose
    %     bst_message_window(sprintf('Image support: %.0f grid points . . .',nsrc));
    % end

    AAt = zeros(nchan); %  Initialize GainCovar matrix
    Gain_RowNorm = [];  %  Initialize matrix containing norms of gain matrix rows (useless and therefore empty if DataType ~= 3)

    for i = 1:OPTIONS.BLOCK:nsrc,

        if OPTIONS.Verbose & ~isempty(OPTIONS.GUIHandles)
            waitbar(i/nsrc,OPTIONS.hwaitbar);
        end

        ndx = [0:(OPTIONS.BLOCK-1)]+i;
        if(ndx(end) > nsrc),  % last block too long
            ndx = [ndx(1):nsrc];
        end

        % Normalize EEG and MEG forward fields separately
        % Fusion Case
        if OPTIONS.DataType == 3 % Compute gain matrix row norm on first block only
            if i == 1
                Gain_RowNorm = zeros(length(OPTIONS.Channel));
                if ~isempty(OPTIONS.GUIHandles)
                    hkk = waitbar(0,'Processing Lead-Field Normalization...');
                end
                for kk = OPTIONS.GoodChannel
                    Gain_RowNorm(kk,kk) = 1/norm(read_gain(gainfile,kk,1:nsrc));
                    if ~isempty(OPTIONS.GUIHandles)
                        waitbar((kk-OPTIONS.GoodChannel(1)+1)/length(OPTIONS.GoodChannel))
                    end

                end

                if ~isempty(OPTIONS.GUIHandles)
                    close(hkk)
                end
            end

            temp = Gain_RowNorm([OPTIONS.GoodChannel],[OPTIONS.GoodChannel])*read_gain(gainfile,[OPTIONS.GoodChannel],ndx);
            OPTIONS.iG = sqrt(norcol(temp)'); % NORCOL computes the column norms to the SQUARE
            OPTIONS.iG(OPTIONS.iG == 0)= min(OPTIONS.iG(OPTIONS.iG ~= 0));
            OPTIONS.iG = spdiags((1./OPTIONS.iG).^OPTIONS.ExpoiG, 0, length(OPTIONS.iG), length(OPTIONS.iG)) ;
            AAt(OPTIONS.GoodChannel,OPTIONS.GoodChannel) = AAt(OPTIONS.GoodChannel,OPTIONS.GoodChannel) + temp*(OPTIONS.iG*OPTIONS.iG)*temp';  % next chunk of correlations

        else % MEG | EEG
            temp = read_gain(gainfile,[OPTIONS.MEGChannel],ndx);
            if ~isempty(OPTIONS.MEGChannel) & ~isempty(temp)
                temp(find(isnan(temp))) = 0;
                if OPTIONS.FFNormalization % Forward field normalization is requested
                    OPTIONS.iG = sqrt(norcol(temp)');  % NORCOL computes the column norms to the SQUARE
                    OPTIONS.iG(OPTIONS.iG == 0)= min(OPTIONS.iG(OPTIONS.iG ~= 0));
                    OPTIONS.iG = spdiags((1./OPTIONS.iG).^OPTIONS.ExpoiG, 0, length(OPTIONS.iG), length(OPTIONS.iG)) ;
                    AAt(OPTIONS.MEGChannel,OPTIONS.MEGChannel) = AAt(OPTIONS.MEGChannel,OPTIONS.MEGChannel) + temp*(OPTIONS.iG*OPTIONS.iG)*temp';  % next chunk of correlations
                else % no FF normalization
                    AAt(OPTIONS.MEGChannel,OPTIONS.MEGChannel) = AAt(OPTIONS.MEGChannel,OPTIONS.MEGChannel) + temp*temp';  % next chunk of correlations
                end
            end

            clear temp OPTIONS.iG

            temp = read_gain(gainfile,OPTIONS.EEGChannel,ndx);
            if ~isempty(OPTIONS.EEGChannel) & ~isempty(temp)
                temp(find(isnan(temp))) = 0;
                if OPTIONS.FFNormalization
                    OPTIONS.iG = sqrt(norcol(temp)');  % NORCOL computes the column norms to the SQUARE
                    if min(OPTIONS.iG(:)) == max(OPTIONS.iG(:)) % as it happens when EEG channels stand for ECoG or EMG
                        OPTIONS.iG = eps;
                    end

                    OPTIONS.iG(OPTIONS.iG == 0) = min(OPTIONS.iG(OPTIONS.iG ~= 0));

                    OPTIONS.iG = spdiags((1./OPTIONS.iG).^OPTIONS.ExpoiG, 0, length(OPTIONS.iG), length(OPTIONS.iG)) ;
                    AAt(OPTIONS.EEGChannel,OPTIONS.EEGChannel) = AAt(OPTIONS.EEGChannel,OPTIONS.EEGChannel) + temp*(OPTIONS.iG*OPTIONS.iG)*temp';  % next chunk of correlations
                else
                    AAt(OPTIONS.EEGChannel,OPTIONS.EEGChannel) = AAt(OPTIONS.EEGChannel,OPTIONS.EEGChannel) + temp*temp';  % next chunk of correlations
                end
            end

            clear temp OPTIONS.iG

        end % datatype

    end % for each block of sources

    varargout{1} = AAt;
    varargout{2} = Gain_RowNorm;

    % -----------------------------------------------------------------------------------------------------------------------

Produced by color_mat2html, a customized BrainStorm 2.0 (Alpha) version of mat2html on Tue Nov 2 16:17:59 2004
Cross-Directory links are: ON