[Master Index]
[Index for Toolbox]
sourceimaging
(Toolbox/sourceimaging.m in BrainStorm 2.0 (Alpha))
Function Synopsis
varargout = sourceimaging(varargin)
Help Text
SOURCEIMAGING - Main switch to source imaging routines
function varargout = sourceimaging(varargin)
function sourceimaging(varargin)
GUI switchyard for cortical imaging routines
(was sourceimaging.m before major upgrade on 07-Aug-2002)
Cross-Reference Information
This function calls
This function is called by
Listing of function C:\BrainStorm_2001\Toolbox\sourceimaging.m
function varargout = sourceimaging(varargin)
%SOURCEIMAGING - Main switch to source imaging routines
% function varargout = sourceimaging(varargin)
% function sourceimaging(varargin)
% GUI switchyard for cortical imaging routines
% (was sourceimaging.m before major upgrade on 07-Aug-2002)
%<autobegin> ---------------------- 12-Oct-2004 12:02:50 -----------------------
% --------- Automatically Generated Comments Block Using AUTO_COMMENTS ---------
%
% CATEGORY: Visualization
%
% Alphabetical list of external functions (non-Matlab):
% toolbox\bst_color_scheme.m
% toolbox\bst_layout.m
% toolbox\bst_message_window.m
% toolbox\bst_sourceimaging.m
% toolbox\bst_static_taskbar.m
% toolbox\bst_win_manager.m
% toolbox\dataplot_cb.m
% toolbox\find_brainstorm_structure.m
% toolbox\get_user_directory.m
% toolbox\headmodeler_gui.m
% toolbox\inorcol.m
% toolbox\load_brainstorm_file.m
% toolbox\mutincomp.m
% toolbox\norcol.m
% toolbox\read_gain.m
% toolbox\readmarkerfile_ctf.m
% toolbox\regsubspace.m
% toolbox\sourceimaging.m NOTE: Routine calls itself explicitly
%
% Subfunctions in this file, in order of occurrence in file:
% AAt = bst_ImageGainCovar(ImageGainFile,BLOCK, VERBOSE);
% [ImageGridAmp, Fsynth] = bst_blk_minnorm(F,AAt,ImageGainFile,nv);
% VecInd = findclosest(VecGuess, VecRef);
%
% Group : Preference data and their calls in this file:
% 'BrainStorm' : 'CurrentData'
%
% setpref('BrainStorm','CurrentData',CurrentData);
%
% CurrentData = getpref('BrainStorm','CurrentData');
%
% Application data and their calls in this file:
% 'DataFiles'
% 'DataMarkerFileTable'
% 'ImagingParam'
% 'ResultFiles'
% 'TileType'
%
% setappdata(SrcImagingGUI,'DataFiles',CurrentData.DataFile);
% setappdata(SrcImagingGUI,'DataFiles',DataFiles);
% setappdata(SrcImagingGUI,'DataMarkerFileTable',DSMarker);
% setappdata(SrcImagingGUI,'ImagingParam',ImagingParam);
% setappdata(SrcImagingGUI,'ResultFiles',OPTIONS.ResultFileSaved)
% setappdata(SrcImagingGUI,'TileType','d');
%
% ImagingParam = getappdata(SrcImagingGUI, 'ImagingParam');
% ImagingParam = getappdata(SrcImagingGUI,'ImagingParam');
% RESULTS = getappdata(fig,'RESULTS');
% [DataFile,tmp] = bst_static_taskbar('GET','DATA');%getappdata(SrcImagingGUI,'ResultFiles');
% [ResultFile,tmp] = bst_static_taskbar('GET','RESULTS');%getappdata(SrcImagingGUI,'ResultFiles');
% fname = getappdata(fig,'RESULTS');
%
% Figure Files opened by this function:
% 'bst_MakeVoisins.fig'
%
% Format of strings below: Type:Style:Tag, "String", CallBack Type and Call
% <automatic> callback is <Tag>_Callback by Matlab default
%
% Callbacks by figure bst_MakeVoisins.fig
% Unable to open bst_MakeVoisins.fig
%
% At Check-in: $Author: Mosher $ $Revision: 25 $ $Date: 10/12/04 10:24a $
%
% 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:02:50 -----------------------
% /---Script Authors--------------------------------------\
% | |
% | *** Sylvain Baillet, Ph.D. |
% | Cognitive Neuroscience & Brain Imaging Laboratory |
% | CNRS UPR640 - LENA |
% | Hopital de la Salpetriere, Paris, France |
% | sylvain.baillet@chups.jussieu.fr |
% | |
% | *** John C. Mosher, Ph.D. |
% | Biophysics Group |
% | Los Alamos National Laboratory |
% | Los Alamos, New Mexico, USA |
% | |
% \-------------------------------------------------------/
%
% Date of creation: Sept-Nov 2001
% (adaptation to BrainStorm of routines by SB, 1997-98) + original features
% History --------------------------------------------------------------------
%
% SB 10-May-2002 Make Manage CurrentData.DataFile as a referential filename
% JCM 28-May-2002 Bug between User.CurrentData and Users.CurrentData.
% Now all are Users.CurrentData
% SB 07-Aug-2002 Starting major upgrade for BrainStorm MMII release
% SB 05-Nov-2002 When GUI is created, display only available surface grid headmodels
% SB 21-Jan-2003 Added minor GUI callbacks
% JCM 08-Sep-2003 Commenting
% SB 07-Oct-2003 Adapting to new window layout and facilitate link to new Viewer
% SB 05-Jan-2004 Secured call when no headmodel is available
% EK 2-Oct-2004 Changed .Tikhonov to be defined as the percentage of the
% maximum sigular value.
% EK 4-Oct-2004 Fusion choise in the GUI and related sections in the code
% are removed
% EK 8-Oct-2004 LCMV beamformer related functions are added
% ----------------------------------------------------------------------------
global data % get global information from Viewer
Users = get_user_directory; % Get current database settings + data sets loaded into memory
CurrentData = getpref('BrainStorm','CurrentData');
VERBOSE = 1;
switch(varargin{1})
case 'create' % Create GUI and load all necessary information
SrcImagingGUI = varargin{2};
setappdata(SrcImagingGUI,'TileType','d');
bst_color_scheme(SrcImagingGUI);
bst_layout('align',SrcImagingGUI)
bst_win_manager(SrcImagingGUI,'Source'); % Define this window as the one in charge of source estimation
handles = guihandles(SrcImagingGUI);
% Define default imaging parameters
ImagingParam.Lamda = str2num(get(handles.cAutomaticParameterSetUp, 'String'));
ImagingParam.MNOnAmplitudes = 1; % spatial regularization based on source amplitudes
ImagingParam.LamdaTime = 0;% Do not temporal regularization
ImagingParam.CortNeighborSyst = 0;% Do not Compute a CorticalNeighborhoodSystem on cortical tessellation
% Use forward field normalization (ie gain matrix is column-normalized)
if isfield(handles, 'FFNormalization')
ImagingParam.FFNormalization = get(handles.FFNormalization, 'Value');
else
ImagingParam.FFNormalization = 0;
end
ImagingParam.ComputeKernel = 0;
setappdata(SrcImagingGUI,'ImagingParam',ImagingParam);
if isfield(handles, 'noiseNormalization') % set the ImagingParam.NNormalization
sourceimaging('noiseNormalization', SrcImagingGUI);
ImagingParam = getappdata(SrcImagingGUI, 'ImagingParam');
end
% Store new parameters in sourceimaging GUI appdata
set(SrcImagingGUI,'Visible','on');
OPTIONS.GUIHandles = handles; % Save for later use when calling bst_sourceimaging.
guidata(SrcImagingGUI,OPTIONS)
%ESEN sourceimaging('Display',ImagingParam) % Display of parameters.
% SrcImagingGUI = bst_win_manager;
% SrcImagingGUI = SrcImagingGUI.Source;
% handles = guihandles(SrcImagingGUI);
% OPTIONS = guidata(SrcImagingGUI);
% Display available headmodels for selected study
if isempty(CurrentData.SubjectFile)
errordlg('Please assign a subject to this dataset first by choosing one from the database.')
return
end
cd(Users.STUDIES)
CurrentData.DataFile = bst_static_taskbar('GET','DATA');
if length(CurrentData.DataFile) == 1 %~isfield(Users.CurrentData,'OtherDataFiles') % Only one Current data file selected
DataFiles = strrep(CurrentData.DataFile,'.mat',''); % Trim out file name for more compact display
% Save data file names in an appdata
setappdata(SrcImagingGUI,'DataFiles',CurrentData.DataFile);
%set(handles.DataFile,'String',DataFiles,'Style','edit')
else
DataFiles = CurrentData.DataFile;
%DataFiles = [DataFiles,strrep(Users.CurrentData.OtherDataFiles,[StudyPath,filesep],'')];
% Save data file names in an appdata
setappdata(SrcImagingGUI,'DataFiles',DataFiles);
DataFiles = strrep(DataFiles,'.mat','');
%set(handles.DataFile,'String',DataFiles,'Style','popupmenu','Max',length(DataFiles),'Value',1)
end
% Available time window ----------------------------------------------------------------------------
% Fill out OPTIONS entries for call to bst_sourceimaging
OPTIONS.Time = data.Display.Time;
OPTIONS.SampleRate = data.Display.SamplingRate;
OPTIONS.TimeSegment = data.Display.Time;
OPTIONS.BaselineSegment = [];
[tmp1,tmp2] = dataplot_cb('GET','time');
OPTIONS.SegmentIndices = tmp2; clear tmp1 tmp2
set(handles.data_start,'String', num2str(1000*OPTIONS.Time(1)))%to ms
set(handles.data_stop,'String', num2str(1000*OPTIONS.Time(2)))
OPTIONS.Method = get(SrcImagingGUI, 'Name'); % see the
guidata(SrcImagingGUI,OPTIONS)
% Check whether file is stored in native format - if yes: fetch marker information
if isempty(CurrentData.DataFile)
errordlg('Please first select a data file in taskbar.','No data selected')
close(SrcImagingGUI)
return
end
whoos = whos('-file',CurrentData.DataFile{1});
iF = strmatch('F',char(whoos.name));
if isempty(iF)
errordlg(...
{sprintf('%s',CurrentData.DataFile{1}),...
'has unrecognized file format.'...
},...
'Wrong data file format')
return
end
if strcmp(whoos(iF).class,'char') % Data in native format
load(CurrentData.DataFile,'F')
k = 1; % in prevision of future loop on a number of data sets loaded into the SourceImaging tool
DataDir{k} = F;
% Check for data file type
[PATHSTR,NAME,EXT,VERSN] = fileparts(DataDir{k});
switch(EXT)
case '.ds' % CTF data file
[DSMarker.names,DSMarker.nsamples ,DSMarker.trial_time] = readmarkerfile_ctf(DataDir{k});
DSMarker.Selected = zeros(length(DSMarker.names),1); % Flag for selection for source estimation
DSMarker.SelectedTrials = cell(length(DSMarker.names),1); % Set of selected trials for each marker
if isempty(DSMarker(1).names) % found no MarkerFile
set(handles.MarkerInformation,'String','No markers available')
set(handles.SelectMarkers,'enable','off')
else
set(handles.SelectMarkers,'enable','on')
setappdata(SrcImagingGUI,'DataMarkerFileTable',DSMarker);
set(handles.MarkerInformation,'String',sprintf('%d marker(s) available', length(DSMarker.names)))
end
disp('')
otherwise
errordlg(...
{sprintf('%s',DataDir{k}),...
'has unrecognized file format.'...
},...
'Wrong data file format')
return
end
end
% ------------------------------------------------------------------------------------------------------
setpref('BrainStorm','CurrentData',CurrentData); % Update prefs
% Look for available headmodels
sourceimaging LHeadModels
case 'LHeadModels' % User has selected a pre-computed headmodel in popup menu - display number of vertices / faces
SrcImagingGUI = bst_win_manager;
SrcImagingGUI = SrcImagingGUI.Source;
handles = guihandles(SrcImagingGUI);
OPTIONS = guidata(SrcImagingGUI);
% Load Headmodel information from model selected headmodel in bst_static_taskbar
[filename,tmp] = bst_static_taskbar('GET','HEADMODEL'); clear tmp
if isempty(filename) % No headmodel available
ButtonName=questdlg('No Head Model is available for this data. Would you like to compute one now ?' ,'No headmodel available',...
'Yes','No','Yes');
switch(ButtonName)
case 'Yes'
Users = get_user_directory;
headmodeler_gui('LoadStudy',Users.CurrentData.StudyFile)
allFigs = bst_win_manager; % Get all figure handles
waitfor(allFigs.Head); % go on with imaging only when headmodelling is done
[filename,tmp] = bst_static_taskbar('GET','HEADMODEL'); clear tmp
case 'No' % nope - cancel computation
delete(SrcImagingGUI) % and close GUI window
return
end
end
OPTIONS.HeadModelFile = {fullfile(Users.STUDIES,filename)};
load(OPTIONS.HeadModelFile{1},'GridName','GridLoc','HeadModelName');
% Now look for the image support -
% Important note - CHEAT (well, not a cheat per se, but I'm addding this keyword for fast retrieval in script):
% As per our new MMII conventions, GridLoc should be a cell array of length one containing
% either a 3xN matrix of source locations or a string pointing at the tessellation file containing the
% tessellation cortical image support.
% If there is an ambiguity in the index of the cortical support in the tessellation file
% then there is a iGrid variable in the headmodel file that specifies this index. Therefore, if iGrid does not
% exist in the headmodel file, this means iGrid = 1 by default.
% Load Tessellation information
if length(GridLoc) > 1 | ~ischar(GridLoc{1})
% msgbox({...
% sprintf('%s is an HeadModel file from an older BrainStorm version.',OPTIONS.HeadModelFile{1}),...
% sprintf('Please compute it again with new Headmodeling tool.')},...
% 'Headmodel file from older BrainStorm version')
errordlg({...
sprintf('Please make sure %s',OPTIONS.HeadModelFile{1}),...
sprintf('is a headmodel of type ''SURF''.'),...
'If not, use BrainStorm''s HeadModel tool to create a new one.'},...
'Wrong Headmodel')
close(SrcImagingGUI)
return
end
OPTIONS.GridLoc = fullfile(Users.SUBJECTS,GridLoc{1});
if ~exist(OPTIONS.GridLoc,'file') % Tessellation file does not exist
msgbox({...
sprintf('Tessellation file %s does not exist on this plateform.',OPTIONS.GridLoc),...
sprintf('Please update your database to have it installed in %s',Users.SUBJECTS)},...
'Missing Tessellation File')
end
% Load Tessellation Information
load(OPTIONS.GridLoc);
HeadModelWhos = whos('-file',OPTIONS.HeadModelFile{1}); % Look for variables in headmodel file
if ~isempty(find(strcmp(cellstr(char(HeadModelWhos(:).name)),'iGrid'))) % See note above
load(OPTIONS.HeadModelFile{1},'iGrid');
else
iGrid = 1;
end
try
Faces = Faces{iGrid}; Vertices = Vertices{iGrid};
catch
errordlg({'Could not find cortical support grid corresponding to selected headmodel file.',...
sprintf('Headmodel file needs grid #%d while there is only %d tessellation(s) in %s.',iGrid,length(Vertices),OPTIONS.GridLoc)...
})
return
end
% Is it a MEG or/and EEG headmodel file
% Initialize GUI checkboxes accordingly
if ~isempty(find(strcmp(cellstr(char(HeadModelWhos(:).name)),'EEGMethod'))) % EEG Headmodel
set(handles.EEG,'Enable','on')
set(handles.EEG,'Enable','on','Value',1)
flagEEG = 1;
else
flagEEG = 0;
set(handles.EEG,'Enable','off','Value',0)
end
if ~isempty(find(strcmp(cellstr(char(HeadModelWhos(:).name)),'MEGMethod'))) % MEG Headmodel
set(handles.MEG,'Enable','on','Value',1)
flagMEG = 1;
else
flagMEG = 0;
set(handles.MEG,'Enable','off','Value',0)
end
guidata(SrcImagingGUI, OPTIONS)
% --------------------------------------------------------------------
case 'Go' % Launch source estimation
SrcImagingGUI = bst_win_manager;
SrcImagingGUI = SrcImagingGUI.Source;
handles = guihandles(SrcImagingGUI);
OPTIONS = guidata(SrcImagingGUI);
ImagingParam = getappdata(SrcImagingGUI, 'ImagingParam');
%dimitrios, read noise segment
% OPTIONS.NNormalization = get(handles.noiseNormalization,'value');
if ImagingParam.NNormalization %if noise normalization is requested
noise_start = get(handles.noise_start,'string');
noise_stop = get(handles.noise_stop,'string');
if(isempty(noise_start) | isempty(noise_stop))
bst_message_window('wrap',{'Noise normalization requires the baseline timeperiod',' '});
return;
end
if ~length( OPTIONS.BaselineSegment)
OPTIONS.BaselineSegment = [str2num(noise_start) str2num(noise_stop)];
end
end
set(SrcImagingGUI,'Pointer','Watch') % Make user wait. . .
% Basic parameters
OPTIONS.Verbose = 1;
DataType = {'MEG','EEG'};
% Now here is a loop operating on every dataset previously selected in the DataManager
% CHEAT - need to enhance things here to manage multiple data sets
OPTIONS.DataFile = fullfile(Users.STUDIES, CurrentData.DataFile{1});
[path,file,ext] = fileparts(OPTIONS.DataFile);
if isempty(ext)
OPTIONS.DataFile = [OPTIONS.DataFile,'.mat'];
end
if ~iscell(OPTIONS.DataFile)
OPTIONS.DataFile = {OPTIONS.DataFile};
end
%GUI.Rank = get(handles.SubspaceSelectionText,'Userdata');
GUI.Rank = []; % Full rank by default
% Select Modality
MEG = get(handles.MEG,'value');
EEG = get(handles.EEG,'value');
OPTIONS.DataTypes = find([MEG,EEG] == 1);
if isempty(OPTIONS.DataTypes ) % No MEG and no EEG selected
errordlg('Please select the proper data type in the Data Viewer window' );
return
end
% Imaging parameters
if isappdata(SrcImagingGUI,'ImagingParam') % Parameters have already been defined
ImagingParam = getappdata(SrcImagingGUI,'ImagingParam');
if ~isfield(ImagingParam,'ComputeKernel')
ImagingParam.ComputeKernel = 0;
end
else
% Define default imaging parameters
ImagingParam.Lamda = str2num(get(handles.cAutomaticParameterSetUp, 'String')); % Spatial regularization
ImagingParam.MNOnAmplitudes = 1; % spatial regularization based on source amplitudes
ImagingParam.LamdaTime = 0;% Do not temporal regularization
ImagingParam.CortNeighborSyst = 0;% Do not Compute a CorticalNeighborhoodSystem on cortical tessellation
ImagingParam.FFNormalization = 1; % Use forward field normalization (ie gain matrix is column-normalized)
ImagingParam.ComputeKernel = 0;
end
OPTIONS.Tikhonov = ImagingParam.Lamda;
OPTIONS.ComputeKernel = ImagingParam.ComputeKernel;
OPTIONS.NNormalization = ImagingParam.NNormalization;
switch ImagingParam.FFNormalization % Check if Forward field normalization is requested
case 1
OPTIONS.FFNormalization = 1;
otherwise
OPTIONS.FFNormalization = 0;
end
% Result file
OPTIONS.ResultFile = 'default';
% Call generic command line function
[Results,OPTIONS] = bst_sourceimaging(OPTIONS);
if ~length(Results)
return
end
clear Results
% end
%ESEN
% Now update result file list in STATIC_TASKBAR and load last result file into Viewer if open
fig = findobj(get(0,'children'),'flat','tag','static_taskbar');
fname = getappdata(fig,'RESULTS');
[tmp,h] = bst_static_taskbar('GET','RESULTS');
% OPTIONS.ResultFileSaved{1} = strrep(OPTIONS.ResultFileSaved{1},fileparts(OPTIONS.ResultFileSaved{1}),'');
% OPTIONS.ResultFileSaved{1} = strrep(OPTIONS.ResultFileSaved{1},'.mat','');
fname{end+1} = OPTIONS.ResultFileSaved{1};
bst_static_taskbar('SET','RESULTS',fname);
set(h,'Value',length(fname))
bst_static_taskbar('popupmenu_RESULTS_Callback',h, [],guihandles(fig), []);
val = get(h,'val'); % where are we in the popup
RESULTS = getappdata(fig,'RESULTS');
global DATAPLOT
if ~ishandle(DATAPLOT) % Viewer is not active
% Load last result file into viewer
% bst_message_window('wrap',...
% {'Loading new results into Viewer',' '});
bst_message_window('wrap',...
{'Processing...',' '});
tmp = bst_static_taskbar('GET','RESULTS');
Users = get_user_directory;
if ~isempty(tmp)
dataplot_cb('loadfile',tmp);
end
end
setappdata(SrcImagingGUI,'ResultFiles',OPTIONS.ResultFileSaved)
set(SrcImagingGUI,'Pointer','Arrow') % That's it !
% ESEN set(handles.HView,'enable','on') % enable View button for result file
sourceimaging('View', SrcImagingGUI)
% sourceimaging('Quit')
%-----------------------------------------------------------------------------------------------------
case 'View'
global DATAPLOT hDataplot % access to Viewer figure handle
figure(DATAPLOT), drawnow
[ResultFile,tmp] = bst_static_taskbar('GET','RESULTS');%getappdata(SrcImagingGUI,'ResultFiles');
if ~length(ResultFile)
bst_message_window('wrap',{'Result file is not available',' '});
return
end
set(DATAPLOT,'pointer','watch')
[DataFile,tmp] = bst_static_taskbar('GET','DATA');%getappdata(SrcImagingGUI,'ResultFiles');
dataplot_cb('loadfile',DataFile);
dataplot_cb('loadfile',ResultFile);
set(hDataplot.SourceViewingSpatialDisplay,'Value',2) % See cortical envelope
dataplot_cb('SourceViewingSpatialDisplay')
dataplot_cb('Plot','measures')
dataplot_cb('Plot','results')
set(DATAPLOT,'pointer','arrow')
%-----------------------------------------------------------------------------------------------------
% case 'selectRegions'
% dataplot_cb('Plot','Measures')
case 'CorticalNeighborhoodSystem' % Design a Cortical Neighborhood for the selected Cortical Surface
% Load information regarding the selected cortical surface
iCortex = get(handles.LHeadModels,'Userdata');
Visu = get(DATAPLOT,'UserData');
if isempty(Visu), errordlg('Please load a dataset using DataManager'), return, end
% Load Tessellation information
try
load(Visu.Tesselation)
catch % Not in proper directory - try something else
cd(Users.SUBJECTS)
load(Visu.Tesselation)
end
VoisinsData.iCortex = iCortex;
VoisinsData.TessellationFile = Visu.Tesselation;
VoisinsData.Vertices = Vertices{iCortex}' ; clear Vertices
VoisinsData.Faces = Faces{iCortex}; clear Faces
VoisinsData.Comment = Comment{iCortex}; clear Comment
MakeVoisins = findobj(0,'tag','MakeVoisins');
if isempty(MakeVoisins)
MakeVoisins = openfig('bst_MakeVoisins.fig');
end
% 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');
VoisinsData.Gain = read_gain(ImageGain{iCortex(get(handles.LHeadModels,'Value'))}{1},GoodChannel,[]);
MV_handles = guihandles(MakeVoisins);
set(MakeVoisins,'Userdata',VoisinsData);
bst_Voisins('create');
case 'changeLambda'
SrcImagingGUI = gcbf;
handles = guihandles(SrcImagingGUI);
if isappdata(SrcImagingGUI,'ImagingParam') % Parameters have already been defined
ImagingParam = getappdata(SrcImagingGUI,'ImagingParam');
end
ImagingParam.Lamda = str2num(get(handles.cAutomaticParameterSetUp, 'String'));
setappdata(SrcImagingGUI,'ImagingParam',ImagingParam);
case 'FFNormalization'
SrcImagingGUI = gcbf;
handles = guihandles(SrcImagingGUI);
if isappdata(SrcImagingGUI,'ImagingParam') % Parameters have already been defined
ImagingParam = getappdata(SrcImagingGUI,'ImagingParam');
end
ImagingParam.FFNormalization = get(handles.FFNormalization, 'Value');
setappdata(SrcImagingGUI,'ImagingParam',ImagingParam);
case 'noiseNormalization'
if length(varargin)>1
SrcImagingGUI = varargin{2};
else
SrcImagingGUI = gcbf;
end
handles = guihandles(SrcImagingGUI);
if isappdata(SrcImagingGUI,'ImagingParam') % Parameters have already been defined
ImagingParam = getappdata(SrcImagingGUI,'ImagingParam');
end
ImagingParam.NNormalization = get(handles.noiseNormalization, 'Value');
if ImagingParam.NNormalization
set (handles.noise_start, 'Enable', 'on');
set (handles.noise_stop, 'Enable', 'on');
set (handles.noise_start,'BackgroundColor', [1 1 1]);
set (handles.noise_stop,'BackgroundColor', [1 1 1]);
else
set (handles.noise_start, 'Enable', 'off');
set (handles.noise_stop, 'Enable', 'off');
set (handles.noise_start,'BackgroundColor', [.753 .753 .753] );
set (handles.noise_stop,'BackgroundColor', [.753 .753 .753]);
end
setappdata(SrcImagingGUI,'ImagingParam',ImagingParam);
case 'noise_start'
SrcImagingGUI = gcbf;
handles = guihandles(SrcImagingGUI);
if ~length(str2num(get(handles.noise_start, 'String')))
bst_message_window('Not a numeric Value')
set(handles.noise_start, 'String', '')
else
OPTIONS = guidata(SrcImagingGUI);
OPTIONS.BaselineSegment(1) = str2num(get(handles.noise_start, 'String'))/1000;%to ms
OPTIONS.Baseline(1) = str2num(get(handles.noise_start, 'String'))/1000;
guidata(SrcImagingGUI, OPTIONS);
end
case 'noise_stop'
SrcImagingGUI = gcbf;
handles = guihandles(SrcImagingGUI);
if ~length(str2num(get(handles.noise_stop, 'String')))
bst_message_window('Not a numeric Value')
set(handles.noise_stop, 'String', '')
else
OPTIONS = guidata(SrcImagingGUI);
OPTIONS.BaselineSegment(2) = str2num(get(handles.noise_stop, 'String'))/1000;%to ms
OPTIONS.Baseline(2) = str2num(get(handles.noise_stop, 'String'))/1000;
guidata(SrcImagingGUI, OPTIONS);
end
case 'data_start'
SrcImagingGUI = gcbf;
handles = guihandles(SrcImagingGUI);
if ~length(str2num(get(handles.data_start, 'String')))
bst_message_window('Not a numeric Value')
set(handles.data_start, 'String', '')
else
OPTIONS = guidata(SrcImagingGUI);
OPTIONS.TimeSegment(1) = str2num(get(handles.data_start, 'String'))/1000;%to ms
OPTIONS.Time(1) = str2num(get(handles.data_start, 'String'))/1000;
guidata(SrcImagingGUI, OPTIONS);
end
case 'data_stop'
SrcImagingGUI = gcbf;
handles = guihandles(SrcImagingGUI);
if ~length(str2num(get(handles.data_stop, 'String')))
bst_message_window('Not a numeric Value')
set(handles.data_stop, 'String', '')
else
OPTIONS = guidata(SrcImagingGUI);
OPTIONS.TimeSegment(2) = str2num(get(handles.data_stop, 'String'))/1000;
OPTIONS.Time(2) = str2num(get(handles.data_stop, 'String'))/1000;
guidata(SrcImagingGUI, OPTIONS);
end
case 'mutincomp_kernel' % Checkbox mutual incompatibility check
SrcImagingGUI = gcbf;
handles = guihandles(SrcImagingGUI);
mutincomp([handles.cBrainStorm, handles.cKernel])
if isappdata(SrcImagingGUI,'ImagingParam') % Parameters have already been defined
ImagingParam = getappdata(SrcImagingGUI,'ImagingParam');
end
if get(handles.cBrainStorm,'Value') % Compute a CorticalNeighborhoodSystem on cortical tessellation
ImagingParam.ComputeKernel = 0;
else
ImagingParam.ComputeKernel = 1;
end
% Store new parameters in sourceimaging GUI appdata
setappdata(SrcImagingGUI,'ImagingParam',ImagingParam);
case 'Tikhonov_regularization'
SrcImagingGUI = gcbf;
handles = guihandles(SrcImagingGUI);
if isappdata(SrcImagingGUI,'ImagingParam') % Parameters have already been defined
ImagingParam = getappdata(SrcImagingGUI,'ImagingParam');
end
if get(handles.TikhonovRegularization, 'Value')
set(handles.cAutomaticParameterSetUp,'Enable', 'on');
set (handles.cAutomaticParameterSetUp,'BackgroundColor', [1 1 1]);
ImagingParam.Lamda = str2num(get(handles.cAutomaticParameterSetUp, 'String'));
set(handles.cAutomaticParameterSetUp, 'String', '10')
else
% set(handles.cAutomaticParameterSetUp,'Value', '');
set(handles.cAutomaticParameterSetUp,'Enable', 'Inactive');
set (handles.cAutomaticParameterSetUp,'BackgroundColor', [.753 .753 .753]);
ImagingParam.Lamda = 0;
end
setappdata(SrcImagingGUI,'ImagingParam',ImagingParam);
% Specifically for LCMV beamformer
case 'output_format_index'
SrcImagingGUI = gcbf;
handles = guihandles(SrcImagingGUI);
set(handles.cBrainStorm, 'Value', 0);
set(handles.normalizedOut, 'Value', 0);
set(handles.noiseNormalization, 'Value', 1);
set(handles.sourcePower, 'Value', 0);
sourceimaging('noiseNormalization', SrcImagingGUI);
% Specifically for LCMV beamformer
case 'output_format_currDensity'
SrcImagingGUI = gcbf;
handles = guihandles(SrcImagingGUI);
% To calculate the current density for LCMV BF compute the spatial
% filter and save it as the Kernel and to get the reconstruction
% map multiply this by the original data Data.F
% Similar to Kernel mode in Min-norm
set(handles.cBrainStorm, 'Value', 1);
set(handles.normalizedOut, 'Value', 0);
set(handles.noiseNormalization, 'Value',0);
set(handles.sourcePower, 'Value', 0);
sourceimaging('noiseNormalization', SrcImagingGUI);
% Specifically for LCMV beamformer
case 'output_format_sourcePower'
SrcImagingGUI = gcbf;
handles = guihandles(SrcImagingGUI);
set(handles.cBrainStorm, 'Value', 0);
set(handles.normalizedOut, 'Value', 0);
set(handles.noiseNormalization, 'Value',0);
set(handles.sourcePower, 'Value', 1);
sourceimaging('noiseNormalization', SrcImagingGUI);
case 'output_format_normalizedOut'
SrcImagingGUI = gcbf;
handles = guihandles(SrcImagingGUI);
set(handles.noiseNormalization, 'Value', 1);
sourceimaging('noiseNormalization', SrcImagingGUI);
set(handles.cBrainStorm, 'Value', 0);
set(handles.normalizedOut, 'Value', 1);
set(handles.noiseNormalization, 'Value',0);
set(handles.sourcePower, 'Value', 0);
case 'Quit'
delete(gcbf);
end
return
set(handles.FT2,'Userdata',source)
set(handles.RPvalue,'String',num2str(lamda,3));
set(handles.Residuals,'String',[num2str(residuals,2),' %']);
%ESEN set(handles.Tstatustxt,'String','3D Visualization on its Way...');
drawnow
newfig = figure;
set(newfig,'Units','Normal','position',[0 .26 .25 .25]);
%pointsource = rezonpatch(source,conffile,ampfile);
gcf = newfig;
[sourceplot] = rezinvolume(source(:,1),conffile,ampfile);
%ESEN [sourceselect,Ierrpot,color] = selekt(sourceplot,G,M(:,current),source(:,1),errpot,handles.Tstatustxt);
set(handles.FTtxt3,'Userdata',sourceplot);
set(handles.sourceorient,'Userdata',sourceselect);
set(saveresult,'Userdata',source)
set(handles.Residuals,'userdata',residuals)
set(handles.status,'Foreg','g');
drawnow
sensorfile = get(handles.sensorfile,'Userdata');
%ESEN set(handles.Tstatustxt,'String','Estimated EEG map');
drawnow
newfig = figure;
set(newfig,'Units','Normal','position',[0 .01 .25 .25]);
gcf = newfig;
[x,y,xs,ys,datarecplot] = carto(newfig,sensorfile,G*source(:,1));
set(handles.FTtxt2,'Userdata',datarecplot);
drawnow
set(handles.RPvalue,'String',num2str(lamda,3));
%ESEN set(handles.Tstatustxt,'String','Ready')
sensitivity = findobj(gcbf,'Tag','Sensitivity');
set(sensitivity,'Visible','on')
set(handles.constraint,'visible','off')
set(handles.fconstraint,'visible','off')
set(handles.mn,'visible','off')
set(handles.mg,'visible','off')
if last ~= first
newfig = figure;
set(newfig,'Units','Normal','position',[.26 .26 .25 .25]);
set(newfig,'Numbertitle','off','Name','handles.Residuals')
plot(time, residuals);
xlabel('time samples')
ylabel('residuals in %')
end
newfig = figure;
set(newfig,'Units','Normal','position',[.52 .26 .25 .25]);
set(newfig,'Numbertitle','off','Name','Sources')
if last > first
xlabel('time samples')
plot(time,source')
else
xlabel('source number')
plot(source')
end
ylabel('source amplitudes')
newfig = figure;
set(newfig,'Units','Normal','position',[.52 .52 .25 .25]);
set(newfig,'Numbertitle','off','Name','handles.data Fit')
if last ~= first
xlabel('time samples')
hest = plot(time,(G*source)');
hold on
horig = plot(time,(M(:,time))','.');
hold off
else
xlabel('sensor label')
bar([M(:,current),(G*source)'],'grouped')
end
ylabel('source amplitudes')
if simu>0
%-- Affichage de la vraie config des sources
datafile = get(handles.datafile,'Userdata');
datafile = [datafile(1:end-4),'.src'];
eval(['cd ',HOME]);
eval(['load ',datafile,' -mat']);
newfig = figure;
set(newfig,'Units','Normal','position',[.26 .01 .2 .2]);
title('Config Originale')
[sourceplot] = rezinvolume(p(:,current),conffile,ampfile);
%ESEN [sourceselect,Ierrpot,color] = selekt(sourceplot,G,M(:,current),p(:,current),errpot,handles.Tstatustxt);
end
clear all
%%----------------------------------------------------------------
%% AFFINNAGE DE LA SENSIBILITE AUX DONNEES
%case 'Sensitivity'
orientations = get(handles.FT3,'Userdata');
if ~isempty(orientations)
delete(orientations)
set(handles.FT3,'Userdata','')
end
M = get(handles.data,'Userdata');
current = str2num(get(handles.CurrentSample,'String'));
M = M(:,current);
G = get(handles.model,'Userdata');
source = get(handles.FT2,'Userdata');
sourceplot = get(handles.FTtxt3,'Userdata');
datarecplot = get(handles.FTtxt2,'Userdata');
% Saisie du seuil de sensibilite
Errpot = findobj(gcbf,'Tag','Errpot');
errpot = str2num(get(Errpot,'String'));
if isempty(errpot)
errpot = 0;
end
% Nettoyage des amplitudes
%ESEN [source, Mrec] = clean(source,M,G,errpot,handles.Tstatustxt);
% Calcul des residus
residuals = 100*norm(M-Mrec)/norm(M);
set(handles.Residuals,'String',[num2str(residuals,2) '%']);
% Affichage
if ~isempty(datarecplot)
eval(['cd ',HOME])
load temp/visucarto.mat
Mes = griddata(x,y,Mrec,xs,ys,'invdist');
set(datarecplot,'Cdata',Mes);
clear Mes
end
eval(['cd ',HOME])
load temp/valpatch
if valpatch == 1
ndip = length(source);
source = reshape(source,3,ndip/3);
magsource = norcol(source);
clear source
maxx = max(magsource(:));
magsource = magsource/maxx;
k = 1;
for dip = 1:length(magsource)
coef = magsource(dip);
if coef > 1/100
color = [coef 0 1-coef];
set(sourceplot(dip),'Markersize',15*coef)
sourceselect(k) = dip;
k = k+1;
set(sourceplot(dip),'markerfacecolor',color);
else
color = 'none';
set(sourceplot(dip),'Markersize',1)
set(sourceplot(dip),'markerfacecolor','none');
end
end
clear magsource
else
%ESEN [sourceselect,Ierrpot,color] = selekt(sourceplot,G,M,source,errpot,handles.Tstatustxt);
end
set(handles.sourceorient,'Userdata',sourceselect);
set(handles.FTtxt3,'Userdata',sourceplot);
set(handles.FTtxt2,'Userdata',datarecplot);
set(saveresult,'Userdata',source);
%% -----------------------------------------------------------------------------
%% Selection des Lead-Fields des Sources les plus actives a priori
%case 'select_lf_gui'
select_lf_gui
%% -----------------------------------------------------------------------------
%% Selection des Sources les plus actives a la resolution precedente
%case 'select_src_cb'
select_src_cb
%% -----------------------------------------------------------------------------
%% Affichage des Orientations des Sources Actives
%case 'handles.sourceorient'
sourceselect = get(handles.sourceorient,'Userdata');
source = get(handles.FT2,'Userdata');
ndip = length(source);
% Amplitude des sources selectionnees
source = reshape(source,3,ndip/3);
source = source(:,sourceselect);
sourceplot = get(handles.FTtxt3,'Userdata');
sourceplot = sourceplot(sourceselect);
X = get(sourceplot,'Xdata');
Y = get(sourceplot,'Ydata');
Z = get(sourceplot,'Zdata');
x = zeros(1,length(sourceplot));
y = x;
z = x;
n = length(sourceplot);
if n > 1
for i = 1:n
x(i) = X{i};
y(i) = Y{i};
z(i) = Z{i};
end
else
x = X;
y = Y;
z = Z;
end
clear X
clear Y
clear Z
figure(get(get(sourceplot(1),'Parent'),'Parent'))
axs = get(sourceplot(1),'Parent');
hold on
source = source*inorcol(source);
orientations = quiver3(x,y,z,source(2,:),source(3,:),source(1,:),1,'.');
set(orientations,'Linewidth',2)
set(handles.FT3,'Userdata',orientations);
hold off
%% --------------------------------------------------------
%case 'saveresult'
source = full(get(saveresult,'Userdata'));
residuals = get(handles.Residuals,'Userdata');
lamda = get(handles.RPvalue,'string');
lamda = str2num(lamda);
lt = str2num(get(handles.valregtemp,'String'));
if isempty(lt)
lt = 0;
end
eval(['cd ',HOME]);
[sourcfile,path] = uiputfile('*.src','Save Sources in...');
if sourcfile ~= 0
eval(['save ',path,sourcfile,' source -ascii'])
eval(['save ',path,sourcfile,'.par residuals lamda lt']);
end
%%------------------------------------------------------
%case 'saveparameters'
%% Sauvetage des parametres en cours
% Fichier de configuration des points source
% matrice de gain
% fichier electrodes
% etc... handles.Tout ce qui est necessaire au lancement
% d'un probleme inverse.
% sauvetage dans [HOME,\temp]
disp('Saving Current Processing Parameters...')
auto = get(handles.RPauto,'value');
if auto == 0
lamda = get(handles.RPvalue,'string');
lamda = str2num(lamda);
else
lamda = '';
end
clear handles.RPauto
% Echantillons temporels
Step = findobj(gcbf,'Tag','Step');
step = str2num(get(Step,'String'));
if step == 0, step = 1; end
if isempty(step), step = 1; end
current = str2num(get(handles.CurrentSample,'String'));
first = str2num(get(handles.From,'String'));
last = str2num(get(handles.To,'String'));
if (isempty(first)|isempty(last)|(first == 0)|(last==0))
first = str2num(get(handles.CurrentSample,'String'));
last = first;
end
% Donnees
datafilename = get(handles.datafile,'String');
datafile = get(handles.datafile,'Userdata');
M = get(handles.data,'Userdata');
sensorfile = get(handles.sensorfile,'Userdata');
sensorfilename = get(handles.sensorfile,'String');
nmes = get(handles.nmes,'Userdata');
% Modele
ndip = get(handles.ndip,'Userdata');
G = get(handles.model,'Userdata');
modelfile = get(handles.modelfile,'Userdata');
modelfilename = get(handles.modelfile,'String');
conffile = get(handles.handles.sourcoepattern,'Userdata');
ampfile = get(handles.ampfile,'Userdata');
% Methode choisie
methode = get(handles.start,'Userdata');
if isempty(methode), methode=''; end
method = get(handles.mn,'value');
Errpot = findobj(gcbf,'Tag','Errpot');
errpot = str2num(get(Errpot,'String'));
if isempty(errpot)
errpot = 0;
end
eval(['cd ',HOME]);
cd temp
save location_parameters
disp('Done')
%%---------------------------------------------------
%case 'loadparameters'
disp('Loading Older Processing Parameters...')
eval(['cd ',HOME]);
cd temp
if exist('load_parameters.mat','file')
errordlg('No File for Parameters Found', 'File not found');
return
end
load location_parameters
% Mise a jour des pointeurs sur GUIs
handles.start = findobj(gcbf,'Tag','handles.start');
handles.ndip = findobj(gcbf,'Tag','ndip');
handles.nmes = findobj(gcbf,'Tag','nmes');
handles.nech = findobj(gcbf,'Tag','nech');
handles.datafile = findobj(gcbf,'Tag','datafile');
handles.modelfile = findobj(gcbf,'Tag','modelfile');
handles.handles.sourcoepattern = findobj(gcbf,'Tag','handles.sourcoepattern');
handles.ampfile = findobj(gcbf,'Tag','ampfile');
handles.sensorfile = findobj(gcbf,'Tag','sensorfile');
handles.model = findobj(gcbf,'Tag','model');
handles.data = findobj(gcbf,'Tag','data');
handles.First = findobj(gcbf,'Tag','handles.First');
handles.Last = findobj(gcbf,'Tag','handles.Last');
handles.SampleSlider = findobj(gcbf,'Tag','handles.SampleSlider');
handles.CurrentSample = findobj(gcbf,'Tag','handles.CurrentSample');
handles.regspace = findobj(gcbf,'Tag','handles.regspace');
handles.regtemp = findobj(gcbf,'Tag','handles.regtemp');
handles.valregtemp = findobj(gcbf,'Tag','handles.valregtemp');
handles.From = findobj(gcbf,'Tag','handles.From');
handles.To = findobj(gcbf,'Tag','handles.To');
handles.Displayhandles.data = findobj(gcbf,'Tag','handles.Displayhandles.data');
% Tableaux de bord des methodes de PI
handles.FrameTikh = findobj(gcbf,'Tag','handles.FrameTikh');
set(handles.FrameTikh,'visible','off')
handles.FT2 = findobj(gcbf,'Tag','handles.FT2');
handles.FT3 = findobj(gcbf,'Tag','handles.FT3');
handles.FT4 = findobj(gcbf,'Tag','handles.FT4');
handles.FT5 = findobj(gcbf,'Tag','handles.FT5');
handles.FTtxt1 = findobj(gcbf,'Tag','handles.FTtxt1');
handles.FTtxt2 = findobj(gcbf,'Tag','handles.FTtxt2');
handles.FTtxt3 = findobj(gcbf,'Tag','handles.FTtxt3');
handles.RPmanual = findobj(gcbf,'Tag','handles.RPmanual');
handles.RPauto = findobj(gcbf,'Tag','handles.RPauto');
handles.RPvalue = findobj(gcbf,'Tag','handles.RPvalue');
handles.Residuals = findobj(gcbf,'Tag','handles.Residuals');
handles.status = findobj(gcbf,'Tag','handles.status');
%ESEN handles.Tstatustxt = findobj(gcbf,'Tag','handles.Tstatustxt');
% if ~isempty(handles.Tstatustxt)
% handles.Tstatustxt = handles.Tstatustxt(1);
% end
%ESEN
handles.constraint= findobj(gcbf,'Tag','constraint');
handles.fconstraint = findobj(gcbf,'Tag','handles.fconstraint');
handles.mn = findobj(gcbf,'Tag','mn');
handles.handles.mg = findobj(gcbf,'Tag','handles.mg');
handles.saveresult= findobj(gcbf,'Tag','saveresult');
% Visu des orientations des dipoles sources
handles.sourceorient = findobj(gcbf,'Tag','handles.sourceorient');
% Visu de la configuration des sources
handles.sourcoepattern = findobj(gcbf,'Tag','sourccoepattern');
%------------------------------------------------------------------------
set(handles.RPauto,'value',auto);
if auto == 0
set(handles.RPvalue,'string',num2str(lamda));
set(handles.RPvalue,'value',1);
else
lamda = '';
set(handles.RPvalue,'value',0);
end
set(handles.datafile,'Userdata',datafile,'String',datafilename,'FontW','demi');
set(handles.data,'Userdata',M);
set(handles.sensorfile,'Userdata',sensorfile);
set(handles.sensorfile,'String',sensorfilename);
[nmes,nech] = size(M);
set(handles.nmes,'Userdata',nmes,'String',num2str(nmes));
set(handles.nech,'Userdata',nech,'String',num2str(nech));
ndip = size(G,2);
set(handles.model,'Userdata',G);
set(handles.handles.sourcoepattern,'Userdata',conffile);
set(handles.ampfile,'Userdata',ampfile);
set(handles.modelfile,'Userdata',modelfile,'string',modelfilename);
set(handles.ndip,'Userdata',ndip,'String',num2str(ndip));
set(handles.handles.sourcoepattern,'Userdata',conffile);
set(handles.ampfile,'Userdata',ampfile);
% Methode choisie
set(handles.start,'Userdata',methode);
Step = findobj(gcbf,'Tag','Step');
if step == 0, step = 1; end
set(Step,'String',num2str(step));
set(handles.First,'String','1');
set(handles.Last,'String',int2str(nech));
set(handles.SampleSlider,'Value',current)
set(handles.SampleSlider,'Min',1,'Max',nech);
set(handles.CurrentSample,'String',int2str(get(handles.SampleSlider,'Value')));
set(handles.From,'String',num2str(first));
set(handles.To,'String',num2str(last));
if (isempty(first)|isempty(last)|(first == 0)|(last==0))
set(handles.CurrentSample,'String',num2str(first));
last = first;
end
set(handles.mn,'value',method);
Errpot = findobj(gcbf,'Tag','Errpot');
set(Errpot,'String',num2str(errpot));
eval(['cd ',HOME]);
eval(['save temp/conffile conffile -mat']);
disp('Done')
% end % switch action
function AAt = bst_ImageGainCovar(ImageGainFile,BLOCK, VERBOSE);
% Block Computation of the ImageGainCovar matrix from gain matrix A
% ImageGainFile : the name of the file containing the ImageGain Matrix
% BLOCK : Number of sources to inlcude in each of the blocks
% VERBOSE : Binary flag to turn on / off the verbose mode during computation
if nargin == 1
BLOCK = 5000;
VERBOSE = 0;
elseif nargin == 2
VERBOSE = 0;
end
fid = fopen(...
IamgeGainFile,...
'r','ieee-be'); %
frewind(fid);
rows = fread(fid,1,'uint32');
AAt = zeros(rows);
if VERBOSE
hwaitbar = waitbar(0,'Creating Gain Matrix Covariance...');
end
for i = 1:BLOCK:nv,
if VERBOSE
waitbar(i/nv,hwaitbar);
end
ndx = [0:(BLOCK-1)]+i;
if(ndx(end) > nv), % last block too long
ndx = [ndx(1):nv];
end
cols = length(ndx); % how many columns to retrieve this time
% 8 bytes per element, find starting point
offset = 4 + (ndx(1)-1)*rows*4;
status = fseek(fid,offset,'bof');
if(status == -1),
error(sprintf('Error reading file at column %.0f',i));
end
temp = fread(fid,[rows,cols],'float32');
AAt = AAt + temp*temp'; % next chunk of correlations
end
if VERBOSE
close(hwaitbar);
end
function [ImageGridAmp, Fsynth] = bst_blk_minnorm(F,AAt,ImageGainFile,nv);
% Begin the minimum norm estimate
% recall the y = Ax, let x = A'c -> y = AA'c, solve for c, the 'coefs' of x.
% c = pinv(AA')*y, and we will regularize the inverse of AA' based on user selection
% Then form x = A'c.
[U,S,V] = svd(AAt);
U = regsubspace(GUI,U,sqrt(S)); % condition based on svd of gain matrix
reg_rank = size(U,2); % reduced rank based on regularization
S = diag(S);
switch deblank(lower(GUI.REG))
case 'tikhonov'
% pinv(A)*b is inv(A'*A)*A'*b. Let A = [A;lamb I]. inv(A'*A + lamb^2 I)*A'*b
Lambda = sqrt(S(1))*GUI.Tikhonov/100; % sing values already squared
Si = 1../(S(1:reg_rank)+Lambda^2); % filtered inverse
otherwise
% do nothing, truncation only
Si = 1../S(1:reg_rank);
end
% the min norm coefs
coefs = V(:,1:reg_rank)*...
((spdiags(Si,0,reg_rank,reg_rank)*U(:,1:reg_rank)') * F);
% let rAAt = regularized(AAt), as found by the above
% So F = AAt *c, c = inv(rAAt), so Fsynth = AAt*c, let's calculate and stick into results
Fsynth = AAt * coefs;
% now calculate the final solution
ImageGridAmp = zeros(nv,size(F,2)); % one column per time slice
hwaitbar = waitbar(0,'Processing final inverse transformation');
disp(sprintf('Processing for %.0f grid points . . .',nv));
fid = fopen(...
IamgeGainFile,...
'r','ieee-be'); %
frewind(fid);
rows = fread(fid,1,'uint32');
for i = 1:BLOCK:nv,
waitbar(i/nv,hwaitbar);
ndx = [0:(BLOCK-1)]+i;
if(ndx(end) > nv), % last block too long
ndx = [ndx(1):nv];
end
cols = length(ndx); % how many columns to retrieve this time
% 8 bytes per element, find starting point
offset = 4 + (ndx(1)-1)*rows*4;
status = fseek(fid,offset,'bof');
if(status == -1),
error(sprintf('Error reading file at column %.0f',i));
end
temp = fread(fid,[rows,cols],'float32');
temp = temp(GoodChannel,:);
ImageGridAmp(ndx,:) = temp'*coefs; % next chunk of solutions
end
close(hwaitbar);
fclose(fid);
return
function VecInd = findclosest(VecGuess, VecRef);
% FINDCLOSEST Find entries of closest elements between two vectors
% function VecInd = findclosest(VecGuess, VecRef);
% Find entries of closest elements between two vectors
% VecGuess is a vector for which one wants to find the closest entries in vector VecRef
% VecInd is the vector of indices pointing atr the entries in vector VecRef that are the closest to VecWin
% VecInd is of the length of VecGuess
%
% In other words, VecRef(VecInd(i)) is the element of VecRef closest to VecGuess(j)
%
% VecRef and VecGuess do not need to be the same length
tmp = repmat(VecRef,1,length(VecGuess));
[minn VecInd] = min(abs(repmat(VecGuess,length(VecRef),1) - tmp));
return
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