[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
This function is called by

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 Oct 12 12:05:14 2004
Cross-Directory links are: ON