[Master Index]
[Index for Toolbox]
bst_sourceimaging
(Toolbox/bst_sourceimaging.m in BrainStorm 2.0 (Alpha))
Function Synopsis
[Results,OPTIONS] = bst_sourceimaging(varargin);
Help Text
BST_SOURCEIMAGING - Command line call to source imaging routines
function [Results,OPTIONS] = bst_sourceimaging(varargin);
function [Results] = bst_sourceimaging(varargin);
[Results] = bst_sourceimaging(OPTIONS)
[OPTIONS] = bst_sourceimaging % Returns the default OPTIONS parameters
[Results,OPTIONS] = bst_sourceimaging(OPTIONS) % Returns also the OPTIONS parameters supplemented by the default settings applied to empty fields in original OPTIONS
INPUTS
OPTIONS - a structure of parameters
OUTPUTS
Results - a regular BrainStorm Results structure (see documentation for details)
Details of the structure for OPTIONS
Data Information:
.DataFile : A string or Cell array of strings of actual data files requested for localization. If isempty, look for data array(s) in .Data
.TimeSegment: Vector of min and max time samples (in seconds) giving first and last latencies specifying
the time window for analysis.
(Default is empty : i.e. take all time samples)
.BaselineTimeSegment: Vector of min and max time samples (in seconds) giving first and last latencies specifying
the time window for baseline. A specific Results file will be created or .BaselineImageGridAmp and .BaselineImageGridTime fields
will be added to the Results output structure
(Default is empty : i.e. no baseline is used).
.BaselineSegment: Again, Vector of min and max time samples (in seconds) giving first and last latencies specifying
the time window for baseline. The baseline will be
used to noise normalize the cortical maps
.ChannelFlag: a vector of flags (1 for good, -1 for bad) of channels in the data file to process
(Default is ChannelFlag from data file or all good channels from data array )
.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.
.Channel : Channel structure 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
NOT RELEVANT : - 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).
NOT IMPLEMENTED YET :
.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)
.FFNormalization : Apply forward field normalization of the gain matrix
(Default is 1)
.Rank : Specify a rank for signal subspace following SVD of the spatio-temporal data matrix
(Default is full rank)
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))
.ResultFileSaved : Cell array of strings that keeps a log of all result file names that were create during call to this function
(Default is left empty, [])
.ComputeKernel : If 1, compute MN kernel to be applied subsequently to data instead of full ImageGridAmp array
CAVEAT: Relevant to Minimum-Norm source mapping only so far
(Default is 0: compute full ImageGridAmp array)
.SaveModelData : 1 :Save model part of data (Fsynth) into ResultFile (if any)
0 : do not save (disk space saving option)
Default is 0 (do not save)
.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).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%-------------------------------------------%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Cross-Reference Information
This function calls
This function is called by
Listing of function C:\BrainStorm_2001\Toolbox\bst_sourceimaging.m
function [Results,OPTIONS] = bst_sourceimaging(varargin);
%BST_SOURCEIMAGING - Command line call to source imaging routines
% function [Results,OPTIONS] = bst_sourceimaging(varargin);
% function [Results] = bst_sourceimaging(varargin);
% [Results] = bst_sourceimaging(OPTIONS)
% [OPTIONS] = bst_sourceimaging % Returns the default OPTIONS parameters
% [Results,OPTIONS] = bst_sourceimaging(OPTIONS) % Returns also the OPTIONS parameters supplemented by the default settings applied to empty fields in original OPTIONS
%
% INPUTS
% OPTIONS - a structure of parameters
% OUTPUTS
% Results - a regular BrainStorm Results structure (see documentation for details)
%
% Details of the structure for OPTIONS
%
% Data Information:
% .DataFile : A string or Cell array of strings of actual data files requested for localization. If isempty, look for data array(s) in .Data
% .TimeSegment: Vector of min and max time samples (in seconds) giving first and last latencies specifying
% the time window for analysis.
% (Default is empty : i.e. take all time samples)
% .BaselineTimeSegment: Vector of min and max time samples (in seconds) giving first and last latencies specifying
% the time window for baseline. A specific Results file will be created or .BaselineImageGridAmp and .BaselineImageGridTime fields
% will be added to the Results output structure
% (Default is empty : i.e. no baseline is used).
% .BaselineSegment: Again, Vector of min and max time samples (in seconds) giving first and last latencies specifying
% the time window for baseline. The baseline will be
% used to noise normalize the cortical maps
% .ChannelFlag: a vector of flags (1 for good, -1 for bad) of channels in the data file to process
% (Default is ChannelFlag from data file or all good channels from data array )
% .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.
% .Channel : Channel structure 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
% NOT RELEVANT : - 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).
% NOT IMPLEMENTED YET :
% .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)
% .FFNormalization : Apply forward field normalization of the gain matrix
% (Default is 1)
% .Rank : Specify a rank for signal subspace following SVD of the spatio-temporal data matrix
% (Default is full rank)
%
% 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))
% .ResultFileSaved : Cell array of strings that keeps a log of all result file names that were create during call to this function
% (Default is left empty, [])
% .ComputeKernel : If 1, compute MN kernel to be applied subsequently to data instead of full ImageGridAmp array
% CAVEAT: Relevant to Minimum-Norm source mapping only so far
% (Default is 0: compute full ImageGridAmp array)
% .SaveModelData : 1 :Save model part of data (Fsynth) into ResultFile (if any)
% 0 : do not save (disk space saving option)
% Default is 0 (do not save)
% .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).
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%-------------------------------------------%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%<autobegin> ---------------------- 12-Oct-2004 12:01:01 -----------------------
% --------- Automatically Generated Comments Block Using AUTO_COMMENTS ---------
%
% CATEGORY: Inverse Modeling
%
% Alphabetical list of external functions (non-Matlab):
% toolbox\bst_message_window.m
% toolbox\bst_sourceimaging.m NOTE: Routine calls itself explicitly
% toolbox\colnorm.m
% toolbox\data_selector_cb.m
% toolbox\find_brainstorm_structure.m
% toolbox\findclosest.m
% toolbox\get_channel.m
% toolbox\get_user_directory.m
% toolbox\good_channel.m
% toolbox\lcmvbf.m
% toolbox\load_brainstorm_file.m
% toolbox\meg4read.m
% toolbox\minnorm.m
% toolbox\norcol.m
% toolbox\read_gain.m
% toolbox\save_fieldnames.m
%
% At Check-in: $Author: Mosher $ $Revision: 29 $ $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:01 -----------------------
% /---Script Author--------------------------------------\
% | |
% | *** 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 ---------------------------------
% SB 23-Dec-2002 Added the .ComputeKernel option
% .Data option is now working
% SB 21-Jan-2003 .ComputeKernel is now working and was tested
% Added .SaveModelData option
% SB 08-Oct-2003 .ResultFile field is properly returned
% EK 4-Oct-2004 Fusion choise in the GUI and related sections in the code
% are removed
% EK 12-Oct-2004 Minor changes
% ----------------------------- Script History ---------------------------------
% Basic parameters
DataType = {'MEG','EEG'}; % Data type labels
Methods = {...
'Minimum-Norm Imaging',...
'Recursive Min. Norm.',...
'ST-MAP',...
'SPTF',...
'RAP-MUSIC'...
}; % Methods Labels
% Default options settings------------------------------------------------------
Def_OPTIONS = struct(...
'BaselineTimeSegment',[],...
'ChannelFlag',[],...
'ComputeKernel',0,...
'Data',[],...
'DataFile','',...
'DataTypes',[],...
'FFNormalization', 1,...
'Gain',[],...
'GainFile','',...
'GridLoc','',...
'GUIHandles',[],...
'iGrid',1,...
'Method','Minimum-Norm Imaging',...
'Rank',[],...
'ResultFile','default',...
'ResultFileSaved',[],...
'SaveModelData',1,...
'Tikhonov',10,...
'TimeSegment',[],...
'Verbose',1 ...
);
% Check arguments passed to and required from function call
if nargin == 0
if nargout > 1
varargout{1} = Def_OPTIONS;
varargout{2} = Def_OPTIONS;
else
varargout{1} = Def_OPTIONS;
end
return
elseif nargin == 1
if ~isstruct(varargin{1})
errordlg('Uncorrect input parameter type, please check the help file'), varargout = cell(nargout,1);
return,
end
OPTIONS = varargin{1};
elseif nargin == 2 % too many input arguments
errordlg('Wrong number of arguments when calling source imaging tool')
varargout = cell(nargout,1);
return
end
% Check field names of passed OPTIONS and fill missing ones with default values
DefFieldNames = fieldnames(Def_OPTIONS);
for k = 1:length(DefFieldNames)
if ~isfield(OPTIONS,DefFieldNames{k})
OPTIONS = setfield(OPTIONS,DefFieldNames{k},getfield(Def_OPTIONS,DefFieldNames{k}));
end
end
clear Def_OPTIONS
% Default options settings : DONE---------------------------------------------------------------------------------------
% Load Data Information-------------------------------------------------------------------------------------------------
UseBaseline = ~isempty(OPTIONS.BaselineTimeSegment); % Was a baseline defined for the data sets ? (1 is yes)
if ~iscell(OPTIONS.HeadModelFile) % A single headmodel file was specified in a string - transform to generic case : a cell
OPTIONS.HeadModelFile = {OPTIONS.HeadModelFile};
end
% Loop operating on every dataset requested for processing
if ~isempty(OPTIONS.DataFile) % Data is stored in regular BsT datafile(s)
DATA = OPTIONS.DataFile;
DataArray = 0; % Flag meaning data is not stored in an array but in a file
elseif ~isempty(OPTIONS.Data) % Direct access to data array(s)
DATA = [];
DataArray = 1; % Flag meaning data is stored in an array, bot in a file
else % No data was given : return
if OPTIONS.Verbose
bst_message_window('No data was passed to BST_SOURCEIMAGING for processing: abort')
end
return
end
if ~iscell(DATA) % Transform DATA in cell array
DATA = {DATA};
end
for kDATA = 1:length(DATA) % Loop on every data set
if ~isempty(OPTIONS.GUIHandles) & OPTIONS.Verbose & ~DataArray
bst_message_window({' ',DATA{kDATA}})
elseif OPTIONS.Verbose & ~DataArray
disp(char( {' ',DATA{kDATA} } ))
elseif OPTIONS.Verbose & DataArray
bst_message_window('Running inverse process on data array. . .')
end
if kDATA == 1 % Time segment and ChannelFlag assumed to be the same for all data sets passed to function
if ~isempty(OPTIONS.TimeSegment) | isempty(OPTIONS.SegmentIndices) % User provided time indices rather than absolute time samples in seconds.
% Default Time Segment: the whole set of time samples
if ~DataArray
load(DATA{kDATA},'Time') % Load available time samples (in seconds)
TimeInd = findclosest([OPTIONS.TimeSegment(1) OPTIONS.TimeSegment(2)],Time'); %get the indices for start and stop points
OPTIONS.SampleRate = mean(diff(Time(TimeInd))); % (in sec)
OPTIONS.TimeSegment = [Time(TimeInd(1)) Time(TimeInd(end))]; % Time segment (min max only, in sec)
OPTIONS.SegmentIndices = TimeInd(1):TimeInd(end); % Full time segment (in sample indices)
if OPTIONS.NNormalization
BaselineInd = findclosest([OPTIONS.BaselineSegment(1) OPTIONS.BaselineSegment(2)],Time'); %get the indices for start and stop points
OPTIONS.SampleNoiseRate = mean(diff(Time(BaselineInd))); % (in sec)
OPTIONS.BaselineSegment = [Time(BaselineInd(1)) Time(BaselineInd(end))]; % Time segment (min max only, in sec)
OPTIONS.NoiseSegmentIndices = BaselineInd(1):BaselineInd(end); % Full Baseline segment (in sample in
if length(OPTIONS.NoiseSegmentIndices)<3
bst_message_window('wrap', {'Baseline region doesnot have enough time slices. Select a different noise region' })
Results = []; % If selected noise region is too small give a msg to change it
return
end
end
clear Time
else
if ~isempty(OPTIONS.Data.Time)
OPTIONS.TimeSegment = OPTIONS.Data.Time;
else
OPTIONS.TimeSegment = [1:size(OPTIONS.Data.F,2)];
end
OPTIONS.SegmentIndices = [1:size(OPTIONS.Data.F,2)];
end
end
if ~DataArray & isempty(OPTIONS.ChannelFlag) % Data file is specified but no ChannelFlag entry in OPTIONS
load(DATA{kDATA},'ChannelFlag')
elseif DataArray & isempty(OPTIONS.ChannelFlag) % Data matrix is passed to function but no ChannelFlag: use all rows.
OPTIONS.ChannelFlag = ones(size(DATA{kDATA},1),1);
end
% Channel information
if ~DataArray
[tmp,channelFileName] = get_channel(DATA{kDATA});
OPTIONS.ChannelFile = fullfile(fileparts(DATA{kDATA}),channelFileName);
else
if isempty(OPTIONS.Data.Channel)
errordlg('A channel structure needs to be specified in OPTIONS.Data.Channel when data is stored in an array')
end
OPTIONS.Channel = OPTIONS.Data.Channel;
end
if isempty(OPTIONS.DataTypes) & ~DataArray % Modality not specified - use all available, therefore load Channel information from channel file
load(OPTIONS.ChannelFile)
% Now find whether we have EEG and MEG available
tmpChannelFlag = ones(length(Channel),1);
MEGChans = good_channel(Channel,ChannelFlag,'MEG');
EEGChans = good_channel(Channel,ChannelFlag,'MEG');
if ~isempty(MEGChans)
OPTIONS.DataTypes = [OPTIONS.DataTypes,1];
end
if ~isempty(EEGChans)
OPTIONS.DataTypes = [OPTIONS.DataTypes,2];
end
elseif isempty(OPTIONS.DataTypes) & DataArray
error('Please define the .DataTypes OPTION when computing inverse solution from a data array')
return
end
end
% Gain / HeadModel information
% Load Headmodel information
if ~isempty(OPTIONS.HeadModelFile{1}) % User passed HeadModel file names
% Now manage situation where user has passed either a single or multiple headmodel files
if length(OPTIONS.HeadModelFile) == 1 % A single headmodel will be used for all data sets
iHMFile = 1; % Still use first headmodel
else
iHMFile = kDATA; % Use headmodel associated to current data set.
end
if exist(OPTIONS.HeadModelFile{iHMFile},'file') % User has passed a valid Headmodelfile name
HeadModelFileWhos = whos('-file',OPTIONS.HeadModelFile{iHMFile});
if isempty(find(strcmp(cellstr(char(HeadModelFileWhos(:).name)),'GridLoc')))
errordlg({...
sprintf('%s is an HeadModel file from an older BrainStorm version.',OPTIONS.HeadModelFile{iHMFile}),...
sprintf('Please compute it again with new Headmodeling tool.')},...
'Headmodel file from older BrainStorm version')
return
end
% Load Image Support Information
load(OPTIONS.HeadModelFile{iHMFile},'GridLoc');
% Load Tessellation information
if length(GridLoc) > 1 %| size(GridLoc{1},1) ~= 3
errordlg({...
sprintf('%s is an HeadModel file from an older BrainStorm version.',OPTIONS.HeadModelFile{iHMFile}),...
sprintf('Please compute it again with new Headmodeling tool.')},...
'Headmodel file from older BrainStorm version')
return
end
if isempty(OPTIONS.GridLoc)
OPTIONS.GridLoc = GridLoc{1};
end
if ischar(OPTIONS.GridLoc) % Source locations are vertices of a tessellated envelope
if ~exist(OPTIONS.GridLoc,'file') % Tessellation file does not exist
if ~isempty(OPTIONS.GUIHandles) % User plays with GUIs
func = 'errordlg'; % Use a error message window
else
func = 'disp';
end
DisplayString = {
sprintf('Tessellation file %s was not found in Matlab''s path.',OPTIONS.GridLoc),...
sprintf('Please update your session''s Matlab path using the GENPATH and ADDPATH commands with the STUDIES and SUBJECTS folders.')...
};
eval([func,'(char(DisplayString))'])
Results = [];
return
end
% Load Tessellation Information
load(OPTIONS.GridLoc);
HeadModelWhos = whos('-file',OPTIONS.HeadModelFile{iHMFile}); % Look for variables in headmodel file
if ~isempty(find(strcmp(cellstr(char(HeadModelWhos(:).name)),'iGrid'))) % See note above
tmp = load(OPTIONS.HeadModelFile{iHMFile},'iGrid');
OPTIONS.iGrid = tmp.iGrid; clear tmp
else
OPTIONS.iGrid = 1;
end
Faces = Faces{OPTIONS.iGrid}; Vertices = Vertices{OPTIONS.iGrid};
% else
% bst_message_window({...
% sprintf('HeadModel File %s was not found',OPTIONS.HeadModelFile{1}),...
% 'Inverse computation aborted.'...
% })
%
end
end
elseif isempty(OPTIONS.Gain)
% No headmodel file, no gain matrix : error
error(...
sprintf('Please specify a HeadModel file or a Gain matrix when calling %s',upper(mfilename))...
)
end % User has passed a BrainStorm HeadModel file
if strcmpi(OPTIONS.ResultFile,'default') % Use default name
ResultFileNameDefault = 1; % Flag meaning default result file name is assumed for every data set
else
ResultFileNameDefault = 0;
end
for tmp = OPTIONS.DataTypes % For each selected modality
OPTIONS.DataType = tmp;
if OPTIONS.Verbose
bst_message_window('wrap', {...
' ',...
OPTIONS.Method,...
sprintf('Processing %s Data. . .',DataType{OPTIONS.DataType})...
})
end
% Name for Result file ------------------------------------------------------------------------------------------------------
if ResultFileNameDefault == 1;%strcmpi(OPTIONS.ResultFile,'default') % Use default name
ResultFileNameDefault = 1; % Flag meaning default result file name is assumed for every data set
switch(OPTIONS.Method)
case 'Minimum-Norm Imaging'
MethodeCode = ['MNE_',DataType{OPTIONS.DataType}];
case 'LCMV Beamformer'
MethodeCode = ['LCMFBF_',DataType{OPTIONS.DataType}];
case 'Recursive Min. Norm.'
MethodeCode = ['RMN_',DataType{OPTIONS.DataType}];;
case 'RAP-MUSIC'
MethodeCode = ['RAP_',DataType{OPTIONS.DataType}];;
case 'SPTF'
MethodeCode = ['SPTF_',DataType{OPTIONS.DataType}];;
otherwise
errordlg('Please choose a method for source localization from the list')
return
end
% now hash together a results name
if ~DataArray
[PATH,NAME,EXT,VER] = fileparts(DATA{kDATA});
NAME = strrep(NAME,'data','results');
else
PATH = pwd;
NAME = 'DataArray';
end
c = clock;
OPTIONS.ResultFile = fullfile(PATH,[NAME sprintf('_%s_%02.0f%02.0f',MethodeCode,c(4:5)) EXT VER]);
i = 0;
while(exist(OPTIONS.ResultFile,'file')),
i = i+1; % subtract another minute
c(5) = mod(c(5) - 1,60);
OPTIONS.ResultFile = fullfile(PATH,[NAME sprintf('_%s_%02.0f%02.0f',MethodeCode,c(4),c(5)) EXT VER]);
end
%OPTIONS.ResultFile = OPTIONS.ResultFile;
if UseBaseline % Associated Baseline Results file
OPTIONSBaseline = OPTIONS; % Fill out a specific OPTIONS structure for inverse processing of Baseline
OPTIONSBaseline.ResultFile = strrep(OPTIONS.ResultFile,'.mat','_baseline.mat');
while(exist(OPTIONSBaseline.ResultFile,'file')),
i = i+1; % subtract another minute
c(5) = mod(c(5) - 1,60);
OPTIONSBaseline.ResultFile = fullfile(PATH,[NAME sprintf('_%s_%02.0f%02.0f',MethodeCode,c(4),c(5)),'_baseline' EXT VER]);
end
%OPTIONSBaseline.Results = OPTIONS.ResultFile;
end
% else % ResultFile name is provided by user and is different than 'default'
% ResultFileNameDefault = 0; %
end
CurrentOPTIONS = OPTIONS; % Duplicate the OPTIONS field before calling imaging routines on current data set
% Alter fields specific to current data set (HeadModels, GridLoc, etc.)
CurrentOPTIONS.HeadModelFile = OPTIONS.HeadModelFile{iHMFile};
if ~DataArray
CurrentOPTIONS.DataFile = DATA{kDATA};
else
CurrentOPTIONS.DataFile = '';
CurrentOPTIONS.Data = OPTIONS.Data;
UseBaseline = 0;
end
if UseBaseline
OPTIONSBaseline.HeadModelFile = OPTIONS.HeadModelFile{iHMFile};
OPTIONSBaseline.DataFile = DATA{kDATA};
OPTIONSBaseline.TimeSegment = OPTIONS.BaselineTimeSegment;
OPTIONSBaseline.SegmentIndices = OPTIONS.BaselineSegmentIndices;
end
% Switchyard to available imaging method-----------------------------------------------------------------------------------
switch( OPTIONS.Method )
case 'Minimum-Norm Imaging'
% Here we go :
Results = minnorm(CurrentOPTIONS);
case 'LCMV Beamformer'
Results = lcmvbf(CurrentOPTIONS);
end
% Compute residuals --------------------------------------------------------------------------------------------------------------
GoodChannel = good_channel(Results.Channel,OPTIONS.ChannelFlag,DataType{OPTIONS.DataType}); %good channels for selected modality
if ~DataArray
load(Results.DataFile,'F'); % Load data
if ischar(F) % Data is in native file format
F = meg4read(F,OPTIONS.Trial,OPTIONS.SegmentIndices,OPTIONS.Verbose);
F = F(GoodChannel,:);
else
% Block of analyzed data
F = F(GoodChannel,OPTIONS.SegmentIndices);
end
else
F = OPTIONS.Data.F;
% Block of analyzed data
F = F(GoodChannel,OPTIONS.SegmentIndices);
end
clear Channel
% ESEN Do it if Fsynth is available--not available for LCMV if
% neural index is calculated
if length(Results.Fsynth)
Results.Residuals = 100*sqrt(norcol(F-Results.Fsynth)./norcol(F)); % norcol return the square of the column norms
rms_res = 100*norm(F-Results.Fsynth,'fro')/norm(F,'fro');
% Compute basic statistics on residuals
mean_res = mean(Results.Residuals);
std_res = std(Results.Residuals);
min_res = min(Results.Residuals);
max_res = max(Results.Residuals);
if OPTIONS.Verbose
DisplayString = {...
'Residuals:',...
sprintf('RMS: %3.1f%%',rms_res),...
sprintf('Std.: %3.1f%%',std_res),...
sprintf('Min.: %3.1f%%',min_res),...
sprintf('Max.: %3.1f%%',max_res)...
};
bst_message_window(DisplayString);
end
end
clear F;
if OPTIONS.SaveModelData == 0
Results.Fsynth = [];
end
if ~isempty(OPTIONS.ResultFile)
% Save Results into file
save_fieldnames(Results,OPTIONS.ResultFile);
bst_message_window('wrap',...
sprintf('Saved in: %s',strrep(OPTIONS.ResultFile,fileparts(OPTIONS.ResultFile),''))...
); % Display only file name without path
end
% -------------------------------------------------------------------------------------------------------------------------------
if UseBaseline % Compute source maps during baseline
if ~isempty(OPTIONS.GUIHandles)
set(OPTIONS.GUIHandles.Tstatustxt,'String','Now computing source maps on baseline segment...'), drawnow
end
if OPTIONS.Verbose
bst_message_window({'','Analysing baseline segment. . .'})
end
if ~isempty(OPTIONS.ResultFile) % if result file is requested from user
Results.ZScore.BaselineFile = OPTIONSBaseline.ResultFile; % append baseline information to original Results structure
save_fieldnames(Results,OPTIONS.ResultFile,'-append');
end
% Run estimation of Baseline source map.
ResultsBaseline = minnorm(OPTIONSBaseline);
% Compute residuals on Baseline------------------------------------------------------------------------------------------------------
load(ResultsBaseline.DataFile,'F'); % Load data
GoodChannel = good_channel(ResultsBaseline.Channel,OPTIONSBaseline.ChannelFlag,DataType{OPTIONSBaseline.DataType}); %good channels for selected modality
% Block of analyzed data
F = F(GoodChannel,[OPTIONSBaseline.SegmentIndices]);
clear Channel
ResultsBaseline.Residuals = 100*sqrt(norcol(F-ResultsBaseline.Fsynth)./norcol(F));
Baselinerms_res = 100*norm(F-ResultsBaseline.Fsynth,'fro')/norm(F,'fro');
clear F;
% Compute basic statistics on residuals
Baselinemean_res = mean(ResultsBaseline.Residuals);
Baselinestd_res = std(ResultsBaseline.Residuals);
Baselinemin_res = min(ResultsBaseline.Residuals);
Baselinemax_res = max(ResultsBaseline.Residuals);
if OPTIONS.SaveModelData == 0
ResultsBaseline.Fsynth = [];
end
if ~isempty(OPTIONSBaseline.ResultFile)
% Save Results into file
save_fieldnames(ResultsBaseline,OPTIONSBaseline.ResultFile);
bst_message_window('wrap',...
sprintf('Baseline Results saved in: %s',strrep(OPTIONSBaseline.ResultFile,fileparts(OPTIONSBaseline.ResultFile),''))...
); % Display only file name without path
end
if OPTIONS.Verbose & UseBaseline
DisplayString{end+1} = 'Residuals (baseline)';
DisplayString{end+1} = sprintf('RMS: %4.2f%% - Std.: %4.2f%%',Baselinerms_res,Baselinestd_res);
DisplayString{end+1} = sprintf('Min.: %4.2f%% - Max.: %4.2f%%',Baselinemin_res,Baselinemax_res);
DisplayString{end+1} = '---------------------------------------------------';
end
if OPTIONS.Verbose
bst_message_window('wrap',DisplayString)
end
end
% Display Residuals performances on analysis time window
% ESEN We don't have Tstatustxt anymore
% if ~isempty(OPTIONS.GUIHandles)
% set(OPTIONS.GUIHandles.Tstatustxt,'String',DisplayString)
% end
%ESEN
%---------------------------------------------------------------------
end % For each selected modality
if ResultFileNameDefault == 1
Users = get_user_directory;
OPTIONS.ResultFileSaved{end+1} = strrep(OPTIONS.ResultFile,[Users.STUDIES,filesep],''); % Keep log of result files
OPTIONS.ResultFile = 'default'; %for next data set
end
end % For each dataset
if OPTIONS.Verbose
bst_message_window({...
' ',...
'---------> Done with Source Imaging'...
})
end
switch( 0) % Older stuff
case 'SPTF'
%Get Study File
%Filename = (Users.CurrentData.StudyFile);
StudySubject = find_brainstorm_structure(Users.CurrentData.StudyFile);
% Now load in the necessary data
cd(Users.SUBJECTS)
if(~isempty(StudySubject.Subject)),
try
Subject = load(fullfile(Users.SUBJECTS,StudySubject.Subject)); %load all of the information on the subject
catch
cd(Users.SUBJECTS)
Subject = load(fullfile(Users.SUBJECTS,StudySubject.Subject)); %load all of the information on the subject
end
cd(Users.SUBJECTS)
load(Subject.Tesselation,'Vertices');
nv = size(Vertices{OPTIONS.iGrid}',1); % number of grid points
else
Subject = [];
end
Study = load_brainstorm_file(StudySubject.Study); % and the study information
% load up the information needed for the gain matrix
cd(Users.STUDIES)
try
if strcmp(OPTIONS.Column_norm,'n')
HeadModel = load(StudySubject.HeadModel,'Function','Param',...
'ImageGridLoc','ImageGain','ImageGainCovar'); % expensive load
else
HeadModel = load(StudySubject.HeadModel,'Function','Param',...
'ImageGridLoc','ImageGain','ImageGainCovar','ImageGainCovar_ColNorm'); % expensive load
end
catch
errordlg([StudySubject.HeadModel, ' be incomplete. Please make sure the forward problem (head model) has been solved for this study'])
return
end
try
if(isempty(HeadModel.ImageGain{OPTIONS.iGrid})),
msgbox('Your Grid Gain Matrix is empty, please run Head Modeler to build it for the selected cortical surface or try another cortex tessellation available','Error','modal');
return
end
% if(isempty(HeadModel.ImageGainCovar{OPTIONS.iGrid})),
% disp(' ')
% disp('Your Grid Gain Correlation matrix is empty, consider running Head Modeler to build');
% disp(' ')
% end
catch
errordlg([StudySubject.HeadModel, ' may be incomplete. Please make sure the forward problem (head model) has been solved for this study'])
return
end
if(~isempty(findstr('_results',deblank(lower(OPTIONS.DataName))))),
% 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.DataName);
% 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
% assume it is a standard data file
Data = load(OPTIONS.DataName,'F','ChannelFlag','Time','Projector','Comment');
end
Channel = load(StudySubject.Channel);
Channel = Channel.Channel;
% now alter the data according to the bad channels
GoodChannel = good_channel(Channel,OPTIONS.ChannelFlag,DataType{OPTIONS.DataType}); %good channels for selected modality
EEGChannel = good_channel(Channel,ones(size(Channel)),'EEG');
MEGChannel = good_channel(Channel,ones(size(Channel)),'MEG');
Data.F = Data.F(GoodChannel,:); % good channels only
%nv = size(Vertices{OPTIONS.iGrid}',1); % number of grid points
GainFile = fullfile(fileparts(StudySubject.HeadModel),HeadModel.ImageGain{OPTIONS.iGrid}{1});
switch(OPTIONS.DataType)
case 1
HeadModel = load(StudySubject.HeadModel,'Function','Param');
G = feval(HeadModel.Function{MEGChannel(1)},Vertices{OPTIONS.iGrid},...
Channel(MEGChannel),HeadModel.Param(MEGChannel),-1);
%% cd(Users.STUDIES)
%% G = read_gain(GainFile,MEGChannel,1:nv); %% Direction constraint
case 2
G = read_gain(GainFile,EEGChannel,1:nv);
case 3
errordlg('Not implemented yet')
end
% Here we go :
set(handles.Residuals,'String','...')
set(handles.Tstatustxt,'String',sprintf('Status:\nQuite busy at the moment...'))
clear ImageGridAmp
%% ImageGridAmp = sptfd(G, Data.F, [OPTIONS.Segment(1)], [OPTIONS.Segment(end)], 10, 1, 6)' ; %% Direction constraint
ImageGridAmp = sptfd(G, Data.F, [OPTIONS.Segment(1)], [OPTIONS.Segment(end)], 10, 1, 6)' ; %% without Direction constraint
Comment = sprintf('SPTF, rank %.0f',OPTIONS.Rank);
Date = datestr(datenum(now),0);
Subject = Subject;
Study = Study;
Time = [OPTIONS.Segment(1):1:OPTIONS.Segment(end)-10];
NoiseCov = [];
SourceCov = [];
Projector = Visu.Data.Projector;
SourceLoc = [];
SourceOrder = [];
SourceOrientation = cell(1);
ModelGain = [];
IndepTopo = [];
TimeSeries = [];
PatchNdx = [];
PatchAmp = [];
ImageGridTime = Visu.Data.Time(Time);
Function = 'sptfa'; % name of this calling routine
Fsynth = [];%G*ImageGridAmp;
Results = struct('Comment',Comment,'Function',mfilename,'StudySubject',StudySubject,...
'GUI',GUI,'Date','','Subject',Subject,'Study',Study,'Time',Time,...
'ChannelFlag',Visu.Data.ChannelFlag,'NoiseCov',NoiseCov,'SourceCov',SourceCov,...
'Projector',Projector,'SourceLoc',SourceLoc,'SourceOrder',SourceOrder,...
'SourceOrientation',[],'ModelGain',ModelGain,...
'IndepTopo',IndepTopo,'TimeSeries',TimeSeries,'PatchNdx',PatchNdx,...
'PatchAmp',PatchAmp,'ImageGridAmp',ImageGridAmp,'ImageGridTime',ImageGridTime,'Fsynth',Fsynth);
Results.SourceOrientation = SourceOrientation; % struct command could kron out
ChannelFlag = Visu.Data.ChannelFlag;
if(isfield(GUI,'Results')), %user provided such a string
if(~isempty(OPTIONS.ResultFile)), % and it's not empty
cd(Users.STUDIES)
save(OPTIONS.ResultFile,'Comment','Function','StudySubject','GUI','Date','Subject','Study','Time',...
'ChannelFlag','NoiseCov','SourceCov','Projector','SourceLoc','SourceOrder',...
'SourceOrientation','ModelGain','IndepTopo','TimeSeries','PatchNdx',...
'PatchAmp','ImageGridAmp','ImageGridTime','Fsynth');
msgbox(sprintf('Wrote results to %s',OPTIONS.ResultFile),'Notice','modal');
end
end
set(handles.Tstatustxt,'String',sprintf('Status:\n\nReady'))
%---------------------------------------------------------------------
case 'Recursive Min. Norm.'
OPTIONS.Lamda = str2num(get(handles.Lamda,'String'));
errpot = .1;
% Load Gain Matrix
StudySubject = find_brainstorm_structure(Users.CurrentData.StudyFile);
[path,file] = fileparts(StudySubject.Study);
cd(fullfile(Users.STUDIES,path));
if(~isempty(StudySubject.Subject)),
Subject = load(fullfile(Users.SUBJECTS,StudySubject.Subject)); %load all of the information on the subject
else
Subject = [];
end
Study = load_brainstorm_file(StudySubject.Study); % and the study information
load(StudySubject.HeadModel,'ImageGain');
iCortex = get(handles.LHeadModels,'Userdata');
G = read_gain(ImageGain{iCortex(get(handles.LHeadModels,'Value'))}{1},GoodChannel,[]);
[ImageGridAmp,lamda,res] = bst_mnef(G,Visu.Data.F(GoodChannel,[OPTIONS.Segment(1):OPTIONS.Segment(end)]),OPTIONS.Lamda,errpot);
% now map to results
Comment = sprintf('FOCAL MIN NORM, rank %.0f',OPTIONS.Rank);
Date = datestr(datenum(now),0);
Subject = Subject;
Study = Study;
Time = [OPTIONS.Segment(1):OPTIONS.Segment(end)];
NoiseCov = [];
SourceCov = [];
Projector = Visu.Data.Projector;
SourceLoc = [];
SourceOrder = [];
SourceOrientation = cell(1);
ModelGain = [];
IndepTopo = [];
TimeSeries = [];
PatchNdx = [];
PatchAmp = [];
ImageGridTime = Visu.Data.Time(Time);
Function = 'focuss'; % name of this calling routine
Fsynth = G*ImageGridAmp;
Results = struct('Comment',Comment,'Function',mfilename,'StudySubject',StudySubject,...
'GUI',GUI,'Date','','Subject',Subject,'Study',Study,'Time',Time,...
'ChannelFlag',Visu.Data.ChannelFlag,'NoiseCov',NoiseCov,'SourceCov',SourceCov,...
'Projector',Projector,'SourceLoc',SourceLoc,'SourceOrder',SourceOrder,...
'SourceOrientation',[],'ModelGain',ModelGain,...
'IndepTopo',IndepTopo,'TimeSeries',TimeSeries,'PatchNdx',PatchNdx,...
'PatchAmp',PatchAmp,'ImageGridAmp',ImageGridAmp,'ImageGridTime',ImageGridTime,'Fsynth',Fsynth);
Results.SourceOrientation = SourceOrientation; % struct command could kron out
ChannelFlag = Visu.Data.ChannelFlag;
if(isfield(GUI,'Results')), %user provided such a string
if(~isempty(OPTIONS.ResultFile)), % and it's not empty
cd(Users.STUDIES)
save(OPTIONS.ResultFile,'Comment','Function','StudySubject','GUI','Date','Subject','Study','Time',...
'ChannelFlag','NoiseCov','SourceCov','Projector','SourceLoc','SourceOrder',...
'SourceOrientation','ModelGain','IndepTopo','TimeSeries','PatchNdx',...
'PatchAmp','ImageGridAmp','ImageGridTime','Fsynth');
msgbox(sprintf('Wrote results to %s',OPTIONS.ResultFile),'Notice','modal');
end
end
%-----------------------------------------------------------------------------------------------------
case 'ST-MAP'
OPTIONS.LamdaTime = str2num(get(handles.LamdaTime,'String'));
if isempty(OPTIONS.LamdaTime)
OPTIONS.LamdaTime = 0;
end
OPTIONS.Lamda = str2num(get(handles.Lamda,'String'));
errpot = .1;
% Load Gain Matrix
StudySubject = find_brainstorm_structure(Users.CurrentData.StudyFile);
[path,file] = fileparts(StudySubject.Study);
cd(fullfile(Users.STUDIES,path));
if(~isempty(StudySubject.Subject)),
Subject = load(fullfile(Users.SUBJECTS,StudySubject.Subject)); %load all of the information on the subject
else
Subject = [];
end
Study = load_brainstorm_file(StudySubject.Study); % and the study information
load(StudySubject.HeadModel,'ImageGain');
iCortex = get(handles.LHeadModels,'Userdata');
G = read_gain(ImageGain{iCortex(get(handles.LHeadModels,'Value'))}{1},GoodChannel,[]);
%Neighborhood System File
cd(Users.SUBJECTS)
NSFile = strrep(Visu.Tesselation,'.mat','_NeighSys.mat');
if exist(NSFile,'file')
load(NSFile)
else
bst_sourceimaging('CorticalNeighborhoodSystem')
end
load(Visu.Tesselation,'Vertices');
Vertices = Vertices{iCortex}';
method = 'mn' ; % 'mn' or 'mg' minnnom or minnorm gradient
[ImageGridAmp, lamda, residuals] = ...
bst_STMAP(G,Visu.Data.F(GoodChannel,[OPTIONS.Segment(1):OPTIONS.Segment(end)]),Vertices,NSFile,OPTIONS.Lamda,OPTIONS.LamdaTime,method,errpot);
case 'RAP-MUSIC'
bst_sourceimaging SubspaceSelection;
data_selector_cb('rap-music')
end
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