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
