[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
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