function [Results,OPTIONS] = bst_sourceimaging(varargin);
%BST_SOURCEIMAGING - Command line call to source imaging routines
% function [Results,OPTIONS] = bst_sourceimaging(varargin);
% function [Results] = bst_sourceimaging(varargin);
% [Results] = bst_sourceimaging(OPTIONS)
% [OPTIONS] = bst_sourceimaging % Returns the default OPTIONS parameters
% [Results,OPTIONS] = bst_sourceimaging(OPTIONS) % Returns also the OPTIONS parameters supplemented by the default settings applied to empty fields in original OPTIONS
%
% INPUTS
%       OPTIONS - a structure of parameters
% OUTPUTS
%       Results - a regular BrainStorm Results structure (see documentation for details)
%
% Details of the structure for OPTIONS
%
%   Data Information:
%          .DataFile   : A string or Cell array of strings of actual data files requested for localization. If isempty, look for data array(s) in .Data
%          .TimeSegment: Vector of min and max time samples (in seconds) giving first and last latencies specifying
%                        the time window for analysis.
%                        (Default is empty : i.e. take all time samples)
%          .BaselineTimeSegment: Vector of min and max time samples (in seconds) giving first and last latencies specifying
%                        the time window for baseline. A specific Results file will be created or .BaselineImageGridAmp and .BaselineImageGridTime fields
%                        will be added to the Results output structure
%                        (Default is empty : i.e. no baseline is used).
%          .BaselineSegment: Again, Vector of min and max time samples (in seconds) giving first and last latencies specifying
%                        the time window for baseline. The baseline will be
%                        used to noise normalize the cortical maps
%          .ChannelFlag: a vector of flags (1 for good, -1 for bad) of channels in the data file to process
%                        (Default is ChannelFlag from data file or all good channels from data array )
%          .DataTypes  : A scalar or a vector of codes to specidy which is the type of data to use:
%                        1 is MEG, 2 is EEG, 3 is Fusion
%                        (Default is all data type available)
%          .FusionType : Code for the type of fusion to apply between MEG and EEG data
%                        1 is linear, 2 is non-linear
%                        (Defaut is 1)
%          .Data       : Specify some data to be passed independently of any data file
%                            .F : Cell array of nChannel x nTime array(s) of data.
%                                 If not empty: uses this array as input for inverse process.
%                            .Channel : Channel structure corresponding to data in .F
%                            .Time : Vector specifying the time instants of data samples in .Data.F
%                        (Default for .Data is empty)
%
%   Source Support Information:
%          .GridLoc    : a string or a 3xN matrix
%                        - if a string : specifies the name of the BrainStorm tessellation file where cortical support is stored
%                        NOT RELEVANT : - if a matrix : contains the coordinates of the cortical locations where to compute the cortical current density
%          .iGrid      : when .GridLoc is a file name, specify which imaging grid to use (Default is 1).
%
%   Forward Field Information :
%          .HeadModelFile: Cell array of strings for filename of the BrainStorm headmodel file(s) to use for imaging. If is empty, look for.GainFile field.
%                          There are only 2 options at this point:
%                                - .HeadModelFile is of length 1, a single headmodel file is used for all data sets passed to the bst_sourceimaging.
%                                - .HeadModelFile is of length the number of data sets passed to bst_sourceimaging (i.e. the length of .DataFile).
%    NOT IMPLEMENTED YET :
%          .GainFile   : Cell array of filenames of image gain matrix (.bin) for the selected cortical support. If is empty, look for .Gain field.
%                        Remarks about the length of .HeadModelFile hold for .GainFile too.
%          .Gain       : Cell array containing the gain matrix of the cortical support (either a single cell or a cell array the length the number of data sets passed to function)
%                        Remarks about the length of .HeadModelFile hold for .GainFile too.
%                        (Default is empty).
%
%   Processing Parameters :
%          .Method     : A string that specifies the imaging method
%                        'Minimum-Norm Imaging': MN estimate of current density using MN priors on source amplitudes (Default)
%                        'Recursive Min. Norm.': Recursive MN estimate based on adaptation of FOCUSS (NOT AVAILABLE)
%                        'SPTF' : Scanning technique constrained on cortical surface (see Denis Schwartz) (UNDER DEVELOPMENT)
%                        'ST-MAP' : Non-linear imaging technique based on sparse-focal image models (NOT AVAILABLE)
%                        'RAP-MUSIC' : link to Mosher's GUI (UNDER DEVELOPMENT)
%                        (Default is 'Minimum-Norm Imaging')
%          .Tikhonov   : Hyperparameter used in Tykhonov regularization
%                        (Default is 10)
%          .FFNormalization : Apply forward field normalization of the gain matrix
%                        (Default is 1)
%          .Rank       : Specify a rank for signal subspace following SVD of the spatio-temporal data matrix
%                        (Default is full rank)
%
%   Output Information:
%          .ResultFile : optional string for filename where to write results
%                        if 'default', filename as used in GUI environment is used
%                        (Default is '' (no Result file created))
%          .ResultFileSaved : Cell array of strings that keeps a log of all result file names that were create during call to this function
%                        (Default is left empty, [])
%          .ComputeKernel : If 1, compute MN kernel to be applied subsequently to data instead of full ImageGridAmp array
%                         CAVEAT: Relevant to Minimum-Norm source mapping only so far
%                         (Default is 0: compute full ImageGridAmp array)
%          .SaveModelData : 1 :Save model part of data (Fsynth) into ResultFile (if any)
%                             0 : do not save (disk space saving option)
%                             Default is 0 (do not save)
%          .Verbose    : Turn ON/OFF verbose mode
%                        (Default is 1 (ON))
%          .GUIHandles : A vector of GUI handles to the SOURCEIMAGING GUI
%                        when operating in graphic mode: passed automatically to bst_sourceimaging by the source imaging switchyard
%                        in command line mode: leave this field empty (default).
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%-------------------------------------------%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%<autobegin> ---------------------- 12-Oct-2004 12:01:01 -----------------------
% --------- Automatically Generated Comments Block Using AUTO_COMMENTS ---------
%
% CATEGORY: Inverse Modeling
%
% Alphabetical list of external functions (non-Matlab):
%   toolbox\bst_message_window.m
%   toolbox\bst_sourceimaging.m  NOTE: Routine calls itself explicitly
%   toolbox\colnorm.m
%   toolbox\data_selector_cb.m
%   toolbox\find_brainstorm_structure.m
%   toolbox\findclosest.m
%   toolbox\get_channel.m
%   toolbox\get_user_directory.m
%   toolbox\good_channel.m
%   toolbox\lcmvbf.m
%   toolbox\load_brainstorm_file.m
%   toolbox\meg4read.m
%   toolbox\minnorm.m
%   toolbox\norcol.m
%   toolbox\read_gain.m
%   toolbox\save_fieldnames.m
%
% At Check-in: $Author: Mosher $  $Revision: 29 $  $Date: 10/12/04 10:23a $
%
% This software is part of BrainStorm Toolbox Version 2.0 (Alpha) 12-Oct-2004
% 
% Principal Investigators and Developers:
% ** Richard M. Leahy, PhD, Signal & Image Processing Institute,
%    University of Southern California, Los Angeles, CA
% ** John C. Mosher, PhD, Biophysics Group,
%    Los Alamos National Laboratory, Los Alamos, NM
% ** Sylvain Baillet, PhD, Cognitive Neuroscience & Brain Imaging Laboratory,
%    CNRS, Hopital de la Salpetriere, Paris, France
% 
% See BrainStorm website at http://neuroimage.usc.edu for further information.
% 
% Copyright (c) 2004 BrainStorm by the University of Southern California
% This software distributed  under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPL
% license can be found at http://www.gnu.org/copyleft/gpl.html .
% 
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%<autoend> ------------------------ 12-Oct-2004 12:01:01 -----------------------



% /---Script Author--------------------------------------\
% |                                                      |
% | *** Sylvain Baillet Ph.D.                            |
% | Cognitive Neuroscience & Brain Imaging Laboratory    |
% | CNRS UPR640 - LENA                                   |
% | Hopital de la Salpetriere, Paris, France             |
% | sylvain.baillet@chups.jussieu.fr                     |
% |                                                      |
% \------------------------------------------------------/
%
% Date of creation: October 2002
%
% ----------------------------- Script History ---------------------------------
% SB 23-Dec-2002  Added the .ComputeKernel option
%                 .Data option is now working
% SB 21-Jan-2003  .ComputeKernel is now working and was tested
%                  Added .SaveModelData option
% SB 08-Oct-2003  .ResultFile field is properly returned
% EK 4-Oct-2004   Fusion choise in the GUI and related sections in the code
%                 are removed
% EK 12-Oct-2004  Minor changes
% ----------------------------- Script History ---------------------------------

% Basic parameters
DataType = {'MEG','EEG'}; % Data type labels
Methods = {...
    'Minimum-Norm Imaging',...
    'Recursive Min. Norm.',...
    'ST-MAP',...
    'SPTF',...
    'RAP-MUSIC'...
    }; % Methods Labels

% Default options settings------------------------------------------------------
Def_OPTIONS = struct(...
    'BaselineTimeSegment',[],...
    'ChannelFlag',[],...
    'ComputeKernel',0,...
    'Data',[],...
    'DataFile','',...
    'DataTypes',[],...
    'FFNormalization', 1,...
    'Gain',[],...
    'GainFile','',...
    'GridLoc','',...
    'GUIHandles',[],...
    'iGrid',1,...
    'Method','Minimum-Norm Imaging',...
    'Rank',[],...
    'ResultFile','default',...
    'ResultFileSaved',[],...
    'SaveModelData',1,...
    'Tikhonov',10,...
    'TimeSegment',[],...
    'Verbose',1 ...
    );


% Check arguments passed to and required from function call

if nargin == 0
    if nargout > 1
        varargout{1} = Def_OPTIONS;
        varargout{2} = Def_OPTIONS;
    else
        varargout{1} = Def_OPTIONS;
    end
    return
elseif nargin == 1
    if ~isstruct(varargin{1})
        errordlg('Uncorrect input parameter type, please check the help file'), varargout = cell(nargout,1);
        return,
    end
    OPTIONS = varargin{1};
elseif nargin == 2 % too many input arguments
    errordlg('Wrong number of arguments when calling source imaging tool')
    varargout = cell(nargout,1);
    return
end

% Check field names of passed OPTIONS and fill missing ones with default values
DefFieldNames = fieldnames(Def_OPTIONS);
for k = 1:length(DefFieldNames)
    if ~isfield(OPTIONS,DefFieldNames{k})

        OPTIONS = setfield(OPTIONS,DefFieldNames{k},getfield(Def_OPTIONS,DefFieldNames{k}));

    end
end

clear Def_OPTIONS

% Default options settings : DONE---------------------------------------------------------------------------------------

% Load Data Information-------------------------------------------------------------------------------------------------

UseBaseline = ~isempty(OPTIONS.BaselineTimeSegment); % Was a baseline defined for the data sets ? (1 is yes)

if ~iscell(OPTIONS.HeadModelFile) % A single headmodel file was specified in a string - transform to generic case : a cell
    OPTIONS.HeadModelFile = {OPTIONS.HeadModelFile};
end


% Loop operating on every dataset requested for processing

if ~isempty(OPTIONS.DataFile) % Data is stored in regular BsT datafile(s)

    DATA = OPTIONS.DataFile;
    DataArray = 0; % Flag meaning data is not stored in an array but in a file

elseif ~isempty(OPTIONS.Data) % Direct access to data array(s)

    DATA = [];
    DataArray = 1; % Flag meaning data is stored in an array, bot in a file

else % No data was given : return
    if OPTIONS.Verbose
        bst_message_window('No data was passed to BST_SOURCEIMAGING for processing: abort')
    end
    return
end

if ~iscell(DATA) % Transform DATA in cell array
    DATA = {DATA};
end


for kDATA = 1:length(DATA) % Loop on every data set


    if ~isempty(OPTIONS.GUIHandles) & OPTIONS.Verbose & ~DataArray
        bst_message_window({' ',DATA{kDATA}})

    elseif OPTIONS.Verbose & ~DataArray
        disp(char( {' ',DATA{kDATA} } ))

    elseif OPTIONS.Verbose & DataArray
        bst_message_window('Running inverse process on data array. . .')
    end

    if kDATA == 1 % Time segment and ChannelFlag assumed to be the same for all data sets passed to function


        if ~isempty(OPTIONS.TimeSegment) | isempty(OPTIONS.SegmentIndices) % User provided time indices rather than absolute time samples in seconds.
            % Default Time Segment: the whole set of time samples
            if ~DataArray
                load(DATA{kDATA},'Time') % Load available time samples (in seconds)
                TimeInd = findclosest([OPTIONS.TimeSegment(1) OPTIONS.TimeSegment(2)],Time'); %get the indices for start and stop points
                OPTIONS.SampleRate = mean(diff(Time(TimeInd))); % (in sec)
                OPTIONS.TimeSegment = [Time(TimeInd(1)) Time(TimeInd(end))]; % Time segment (min max only, in sec)
                OPTIONS.SegmentIndices = TimeInd(1):TimeInd(end); % Full time segment (in sample indices)

                if OPTIONS.NNormalization
                    BaselineInd = findclosest([OPTIONS.BaselineSegment(1) OPTIONS.BaselineSegment(2)],Time'); %get the indices for start and stop points
                    OPTIONS.SampleNoiseRate = mean(diff(Time(BaselineInd))); % (in sec)
                    OPTIONS.BaselineSegment = [Time(BaselineInd(1)) Time(BaselineInd(end))]; % Time segment (min max only, in sec)
                    OPTIONS.NoiseSegmentIndices = BaselineInd(1):BaselineInd(end); % Full Baseline segment (in sample in
                    if length(OPTIONS.NoiseSegmentIndices)<3
                        bst_message_window('wrap', {'Baseline region doesnot have enough time slices. Select a different noise region'                })
                        Results = []; % If selected noise region is too small give a msg to change it
                        return
                    end
                    
                end

                clear Time
            else
                if ~isempty(OPTIONS.Data.Time)
                    OPTIONS.TimeSegment = OPTIONS.Data.Time;
                else
                    OPTIONS.TimeSegment = [1:size(OPTIONS.Data.F,2)];
                end
                OPTIONS.SegmentIndices = [1:size(OPTIONS.Data.F,2)];
            end

        end

        if ~DataArray & isempty(OPTIONS.ChannelFlag) % Data file is specified but no ChannelFlag entry in OPTIONS
            load(DATA{kDATA},'ChannelFlag')
        elseif DataArray & isempty(OPTIONS.ChannelFlag) % Data matrix is passed to function but no ChannelFlag: use all rows.
            OPTIONS.ChannelFlag = ones(size(DATA{kDATA},1),1);
        end

        % Channel information
        if ~DataArray

            [tmp,channelFileName] = get_channel(DATA{kDATA});
            OPTIONS.ChannelFile = fullfile(fileparts(DATA{kDATA}),channelFileName);

        else
            if isempty(OPTIONS.Data.Channel)
                errordlg('A channel structure needs to be specified in OPTIONS.Data.Channel when data is stored in an array')
            end
            OPTIONS.Channel = OPTIONS.Data.Channel;
        end

        if isempty(OPTIONS.DataTypes) & ~DataArray % Modality not specified - use all available, therefore load Channel information from channel file
            load(OPTIONS.ChannelFile)
            % Now find whether we have EEG and MEG available
            tmpChannelFlag = ones(length(Channel),1);
            MEGChans = good_channel(Channel,ChannelFlag,'MEG');
            EEGChans = good_channel(Channel,ChannelFlag,'MEG');
            if ~isempty(MEGChans)
                OPTIONS.DataTypes = [OPTIONS.DataTypes,1];
            end
            if ~isempty(EEGChans)
                OPTIONS.DataTypes = [OPTIONS.DataTypes,2];
            end
 
        elseif isempty(OPTIONS.DataTypes) & DataArray
            error('Please define the .DataTypes OPTION when computing inverse solution from a data array')
            return
        end
    end

    % Gain / HeadModel information
    % Load Headmodel information
    if ~isempty(OPTIONS.HeadModelFile{1}) % User passed HeadModel file names
        % Now manage situation where user has passed either a single or multiple headmodel files
        if length(OPTIONS.HeadModelFile) == 1 % A single headmodel will be used for all data sets
            iHMFile = 1; % Still use first headmodel
        else
            iHMFile = kDATA; % Use headmodel associated to current data set.
        end

        if exist(OPTIONS.HeadModelFile{iHMFile},'file') % User has passed a valid Headmodelfile name
            HeadModelFileWhos = whos('-file',OPTIONS.HeadModelFile{iHMFile});
            if isempty(find(strcmp(cellstr(char(HeadModelFileWhos(:).name)),'GridLoc')))
                errordlg({...
                    sprintf('%s is an HeadModel file from an older BrainStorm version.',OPTIONS.HeadModelFile{iHMFile}),...
                    sprintf('Please compute it again with new Headmodeling tool.')},...
                    'Headmodel file from older BrainStorm version')
                return
            end

            % Load Image Support Information
            load(OPTIONS.HeadModelFile{iHMFile},'GridLoc');
            % Load Tessellation information
            if length(GridLoc) > 1 %| size(GridLoc{1},1) ~= 3
                errordlg({...
                    sprintf('%s is an HeadModel file from an older BrainStorm version.',OPTIONS.HeadModelFile{iHMFile}),...
                    sprintf('Please compute it again with new Headmodeling tool.')},...
                    'Headmodel file from older BrainStorm version')
                return
            end

            if isempty(OPTIONS.GridLoc)
                OPTIONS.GridLoc = GridLoc{1};
            end

            if ischar(OPTIONS.GridLoc) % Source locations are vertices of a tessellated envelope
                if ~exist(OPTIONS.GridLoc,'file') % Tessellation file does not exist
                    if ~isempty(OPTIONS.GUIHandles) % User plays with GUIs
                        func = 'errordlg'; % Use a error message window
                    else
                        func = 'disp';
                    end
                    DisplayString = {
                        sprintf('Tessellation file %s was not found in Matlab''s path.',OPTIONS.GridLoc),...
                        sprintf('Please update your session''s Matlab path using the GENPATH and ADDPATH commands with the STUDIES and SUBJECTS folders.')...
                        };

                    eval([func,'(char(DisplayString))'])
                    Results = [];
                    return
                end

                % Load Tessellation Information
                load(OPTIONS.GridLoc);
                HeadModelWhos = whos('-file',OPTIONS.HeadModelFile{iHMFile}); % Look for variables in headmodel file
                if ~isempty(find(strcmp(cellstr(char(HeadModelWhos(:).name)),'iGrid'))) % See note above
                    tmp = load(OPTIONS.HeadModelFile{iHMFile},'iGrid');
                    OPTIONS.iGrid = tmp.iGrid; clear tmp
                else
                    OPTIONS.iGrid = 1;
                end
                Faces = Faces{OPTIONS.iGrid}; Vertices = Vertices{OPTIONS.iGrid};

                %             else
                %                 bst_message_window({...
                %                         sprintf('HeadModel File %s was not found',OPTIONS.HeadModelFile{1}),...
                %                         'Inverse computation aborted.'...
                %                     })
                %
            end
        end

    elseif isempty(OPTIONS.Gain)
        % No headmodel file, no gain matrix : error
        error(...
            sprintf('Please specify a HeadModel file or a Gain matrix when calling %s',upper(mfilename))...
            )
    end % User has passed a BrainStorm HeadModel file

    if strcmpi(OPTIONS.ResultFile,'default') % Use default name
        ResultFileNameDefault = 1; % Flag meaning default result file name is assumed for every data set
    else
        ResultFileNameDefault = 0;
    end

    for tmp = OPTIONS.DataTypes % For each selected modality

        OPTIONS.DataType = tmp;

        if OPTIONS.Verbose
            bst_message_window('wrap', {...
                ' ',...
                OPTIONS.Method,...
                sprintf('Processing %s Data. . .',DataType{OPTIONS.DataType})...
                })
        end

        % Name for Result file ------------------------------------------------------------------------------------------------------
        if ResultFileNameDefault == 1;%strcmpi(OPTIONS.ResultFile,'default') % Use default name

            ResultFileNameDefault = 1; % Flag meaning default result file name is assumed for every data set

            switch(OPTIONS.Method)
                case 'Minimum-Norm Imaging'
                    MethodeCode = ['MNE_',DataType{OPTIONS.DataType}];
                case 'LCMV Beamformer'
                    MethodeCode = ['LCMFBF_',DataType{OPTIONS.DataType}];
                case 'Recursive Min. Norm.'
                    MethodeCode = ['RMN_',DataType{OPTIONS.DataType}];;
                case 'RAP-MUSIC'
                    MethodeCode = ['RAP_',DataType{OPTIONS.DataType}];;
                case 'SPTF'
                    MethodeCode = ['SPTF_',DataType{OPTIONS.DataType}];;
                otherwise
                    errordlg('Please choose a method for source localization from the list')
                    return
            end

            % now hash together a results name
            if ~DataArray
                [PATH,NAME,EXT,VER] = fileparts(DATA{kDATA});
                NAME = strrep(NAME,'data','results');
            else
                PATH = pwd;
                NAME = 'DataArray';
            end

            c = clock;
            OPTIONS.ResultFile = fullfile(PATH,[NAME sprintf('_%s_%02.0f%02.0f',MethodeCode,c(4:5)) EXT VER]);
            i = 0;
            while(exist(OPTIONS.ResultFile,'file')),
                i = i+1; % subtract another minute
                c(5) = mod(c(5) - 1,60);
                OPTIONS.ResultFile = fullfile(PATH,[NAME sprintf('_%s_%02.0f%02.0f',MethodeCode,c(4),c(5)) EXT VER]);
            end
            %OPTIONS.ResultFile = OPTIONS.ResultFile;

            if UseBaseline % Associated Baseline Results file
                OPTIONSBaseline = OPTIONS; % Fill out a specific OPTIONS structure for inverse processing of Baseline
                OPTIONSBaseline.ResultFile = strrep(OPTIONS.ResultFile,'.mat','_baseline.mat');
                while(exist(OPTIONSBaseline.ResultFile,'file')),
                    i = i+1; % subtract another minute
                    c(5) = mod(c(5) - 1,60);
                    OPTIONSBaseline.ResultFile = fullfile(PATH,[NAME sprintf('_%s_%02.0f%02.0f',MethodeCode,c(4),c(5)),'_baseline' EXT VER]);
                end
                %OPTIONSBaseline.Results = OPTIONS.ResultFile;
            end

            %         else % ResultFile name is provided by user and is different than 'default'
            %             ResultFileNameDefault = 0; %
        end

        CurrentOPTIONS = OPTIONS; % Duplicate the OPTIONS field before calling imaging routines on current data set
        % Alter fields specific to current data set (HeadModels, GridLoc, etc.)
        CurrentOPTIONS.HeadModelFile = OPTIONS.HeadModelFile{iHMFile};
        if ~DataArray
            CurrentOPTIONS.DataFile = DATA{kDATA};
        else
            CurrentOPTIONS.DataFile = '';
            CurrentOPTIONS.Data = OPTIONS.Data;
            UseBaseline = 0;
        end

        if UseBaseline
            OPTIONSBaseline.HeadModelFile = OPTIONS.HeadModelFile{iHMFile};
            OPTIONSBaseline.DataFile = DATA{kDATA};
            OPTIONSBaseline.TimeSegment = OPTIONS.BaselineTimeSegment;
            OPTIONSBaseline.SegmentIndices = OPTIONS.BaselineSegmentIndices;
        end

        % Switchyard to available imaging method-----------------------------------------------------------------------------------
        switch( OPTIONS.Method )

            case 'Minimum-Norm Imaging'

                % Here we go :
                Results = minnorm(CurrentOPTIONS);
            case 'LCMV Beamformer'
                Results = lcmvbf(CurrentOPTIONS);


        end

        % Compute residuals --------------------------------------------------------------------------------------------------------------
        GoodChannel = good_channel(Results.Channel,OPTIONS.ChannelFlag,DataType{OPTIONS.DataType}); %good channels for selected modality
        if ~DataArray
            load(Results.DataFile,'F'); % Load data
            if ischar(F) % Data is in native file format
                F = meg4read(F,OPTIONS.Trial,OPTIONS.SegmentIndices,OPTIONS.Verbose);
                F = F(GoodChannel,:);
            else
                % Block of analyzed data
                F = F(GoodChannel,OPTIONS.SegmentIndices);
            end
        else
            F = OPTIONS.Data.F;
            % Block of analyzed data
            F = F(GoodChannel,OPTIONS.SegmentIndices);
        end



        clear Channel
        % ESEN Do it if Fsynth is available--not available for LCMV if
        % neural index is calculated
        if length(Results.Fsynth)
            Results.Residuals = 100*sqrt(norcol(F-Results.Fsynth)./norcol(F)); % norcol return the square of the column norms
            rms_res = 100*norm(F-Results.Fsynth,'fro')/norm(F,'fro');
            % Compute basic statistics on residuals
            mean_res = mean(Results.Residuals);
            std_res = std(Results.Residuals);
            min_res = min(Results.Residuals);
            max_res = max(Results.Residuals);
            if OPTIONS.Verbose
                DisplayString = {...
                    'Residuals:',...
                    sprintf('RMS: %3.1f%%',rms_res),...
                    sprintf('Std.: %3.1f%%',std_res),...
                    sprintf('Min.: %3.1f%%',min_res),...
                    sprintf('Max.: %3.1f%%',max_res)...
                    };

                bst_message_window(DisplayString);
            end

        end
        clear F;





        if OPTIONS.SaveModelData == 0
            Results.Fsynth = [];
        end

        if ~isempty(OPTIONS.ResultFile)
            % Save Results into file
            save_fieldnames(Results,OPTIONS.ResultFile);
            bst_message_window('wrap',...
                sprintf('Saved in: %s',strrep(OPTIONS.ResultFile,fileparts(OPTIONS.ResultFile),''))...
                ); % Display only file name without path
        end

        % -------------------------------------------------------------------------------------------------------------------------------

        if UseBaseline % Compute source maps during baseline
            if ~isempty(OPTIONS.GUIHandles)
                set(OPTIONS.GUIHandles.Tstatustxt,'String','Now computing source maps on baseline segment...'), drawnow
            end

            if OPTIONS.Verbose
                bst_message_window({'','Analysing baseline segment. . .'})
            end

            if ~isempty(OPTIONS.ResultFile) % if result file is requested from user
                Results.ZScore.BaselineFile = OPTIONSBaseline.ResultFile; % append baseline information to original Results structure
                save_fieldnames(Results,OPTIONS.ResultFile,'-append');
            end

            % Run estimation of Baseline source map.
            ResultsBaseline = minnorm(OPTIONSBaseline);

            % Compute residuals on Baseline------------------------------------------------------------------------------------------------------
            load(ResultsBaseline.DataFile,'F'); % Load data
            GoodChannel = good_channel(ResultsBaseline.Channel,OPTIONSBaseline.ChannelFlag,DataType{OPTIONSBaseline.DataType}); %good channels for selected modality

            % Block of analyzed data
            F = F(GoodChannel,[OPTIONSBaseline.SegmentIndices]);

            clear Channel
            ResultsBaseline.Residuals = 100*sqrt(norcol(F-ResultsBaseline.Fsynth)./norcol(F));
            Baselinerms_res = 100*norm(F-ResultsBaseline.Fsynth,'fro')/norm(F,'fro');
            clear F;

            % Compute basic statistics on residuals
            Baselinemean_res = mean(ResultsBaseline.Residuals);
            Baselinestd_res = std(ResultsBaseline.Residuals);
            Baselinemin_res = min(ResultsBaseline.Residuals);
            Baselinemax_res = max(ResultsBaseline.Residuals);

            if OPTIONS.SaveModelData == 0
                ResultsBaseline.Fsynth = [];
            end

            if ~isempty(OPTIONSBaseline.ResultFile)
                % Save Results into file
                save_fieldnames(ResultsBaseline,OPTIONSBaseline.ResultFile);
                bst_message_window('wrap',...
                    sprintf('Baseline Results saved in: %s',strrep(OPTIONSBaseline.ResultFile,fileparts(OPTIONSBaseline.ResultFile),''))...
                    ); % Display only file name without path
            end

            if OPTIONS.Verbose & UseBaseline
                DisplayString{end+1} = 'Residuals (baseline)';
                DisplayString{end+1} = sprintf('RMS: %4.2f%% - Std.: %4.2f%%',Baselinerms_res,Baselinestd_res);
                DisplayString{end+1} = sprintf('Min.: %4.2f%% - Max.: %4.2f%%',Baselinemin_res,Baselinemax_res);
                DisplayString{end+1} = '---------------------------------------------------';
            end

            if OPTIONS.Verbose
                bst_message_window('wrap',DisplayString)
            end

        end

        % Display Residuals performances on analysis time window

        % ESEN We don't have Tstatustxt anymore
        %         if ~isempty(OPTIONS.GUIHandles)
        %             set(OPTIONS.GUIHandles.Tstatustxt,'String',DisplayString)
        %         end
        %ESEN

        %---------------------------------------------------------------------

    end % For each selected modality

    if ResultFileNameDefault == 1
        Users = get_user_directory;
        OPTIONS.ResultFileSaved{end+1} = strrep(OPTIONS.ResultFile,[Users.STUDIES,filesep],''); % Keep log of result files
        OPTIONS.ResultFile = 'default'; %for next data set
    end

end % For each dataset

if OPTIONS.Verbose
    bst_message_window({...
        ' ',...
        '---------> Done with Source Imaging'...
        })
end



switch( 0) % Older stuff

    case 'SPTF'
        %Get Study File
        %Filename = (Users.CurrentData.StudyFile);
        StudySubject = find_brainstorm_structure(Users.CurrentData.StudyFile);


        % Now load in the necessary data

        cd(Users.SUBJECTS)
        if(~isempty(StudySubject.Subject)),
            try
                Subject = load(fullfile(Users.SUBJECTS,StudySubject.Subject)); %load all of the information on the subject
            catch
                cd(Users.SUBJECTS)
                Subject = load(fullfile(Users.SUBJECTS,StudySubject.Subject)); %load all of the information on the subject
            end

            cd(Users.SUBJECTS)
            load(Subject.Tesselation,'Vertices');
            nv = size(Vertices{OPTIONS.iGrid}',1); % number of grid points
        else
            Subject = [];
        end

        Study = load_brainstorm_file(StudySubject.Study); % and the study information

        % load up the information needed for the gain matrix
        cd(Users.STUDIES)
        try
            if strcmp(OPTIONS.Column_norm,'n')
                HeadModel = load(StudySubject.HeadModel,'Function','Param',...
                    'ImageGridLoc','ImageGain','ImageGainCovar'); % expensive load
            else
                HeadModel = load(StudySubject.HeadModel,'Function','Param',...
                    'ImageGridLoc','ImageGain','ImageGainCovar','ImageGainCovar_ColNorm'); % expensive load


            end

        catch

            errordlg([StudySubject.HeadModel, '  be incomplete. Please make sure the forward problem (head model) has been solved for this study'])
            return

        end

        try
            if(isempty(HeadModel.ImageGain{OPTIONS.iGrid})),
                msgbox('Your Grid Gain Matrix is empty, please run Head Modeler to build it for the selected cortical surface or try another cortex tessellation available','Error','modal');
                return
            end
            %     if(isempty(HeadModel.ImageGainCovar{OPTIONS.iGrid})),
            %         disp(' ')
            %         disp('Your Grid Gain Correlation matrix is empty, consider running Head Modeler to build');
            %         disp(' ')
            %     end
        catch
            errordlg([StudySubject.HeadModel, ' may be incomplete. Please make sure the forward problem (head model) has been solved for this study'])
            return
        end

        if(~isempty(findstr('_results',deblank(lower(OPTIONS.DataName))))),
            % CHEAT
            % The DataName has the string '_results' in it
            % We are trying to treat independent topographies as min norm data. Eventually this
            %  capability will be moved to a new routine.
            % name has the string '_results' in it
            Results = load(OPTIONS.DataName);
            % now synthesize a Data structure from this
            [ignorm,Anrm] = colnorm(Results.IndepTopo);

            Data = struct('F',Anrm,'ChannelFlag',Results.ChannelFlag,...
                'Time',[1:size(Results.IndepTopo,2)],... % time is now simply an indexer
                'NoiseCov',Results.NoiseCov,'SourceCov',Results.SourceCov,'Projector',Results.Projector,...
                'Comment',sprintf('Synthesized Data Set from Results: %s',Results.Comment));
        else
            % assume it is a standard data file
            Data = load(OPTIONS.DataName,'F','ChannelFlag','Time','Projector','Comment');
        end

        Channel = load(StudySubject.Channel);
        Channel = Channel.Channel;

        % now alter the data according to the bad channels

        GoodChannel = good_channel(Channel,OPTIONS.ChannelFlag,DataType{OPTIONS.DataType}); %good channels for selected modality

        EEGChannel = good_channel(Channel,ones(size(Channel)),'EEG');
        MEGChannel = good_channel(Channel,ones(size(Channel)),'MEG');

        Data.F = Data.F(GoodChannel,:); % good channels only

        %nv = size(Vertices{OPTIONS.iGrid}',1); % number of grid points
        GainFile = fullfile(fileparts(StudySubject.HeadModel),HeadModel.ImageGain{OPTIONS.iGrid}{1});

        switch(OPTIONS.DataType)
            case 1
                HeadModel = load(StudySubject.HeadModel,'Function','Param');
                G = feval(HeadModel.Function{MEGChannel(1)},Vertices{OPTIONS.iGrid},...
                    Channel(MEGChannel),HeadModel.Param(MEGChannel),-1);

                %% cd(Users.STUDIES)
                %% G = read_gain(GainFile,MEGChannel,1:nv); %% Direction constraint
            case 2
                G = read_gain(GainFile,EEGChannel,1:nv);
            case 3
                errordlg('Not implemented yet')
        end

        % Here we go :
        set(handles.Residuals,'String','...')

        set(handles.Tstatustxt,'String',sprintf('Status:\nQuite busy at the moment...'))

        clear ImageGridAmp

        %% ImageGridAmp = sptfd(G, Data.F, [OPTIONS.Segment(1)], [OPTIONS.Segment(end)], 10, 1, 6)' ; %% Direction constraint
        ImageGridAmp = sptfd(G, Data.F, [OPTIONS.Segment(1)], [OPTIONS.Segment(end)], 10, 1, 6)' ; %% without Direction constraint


        Comment = sprintf('SPTF, rank %.0f',OPTIONS.Rank);
        Date = datestr(datenum(now),0);
        Subject = Subject;
        Study = Study;
        Time = [OPTIONS.Segment(1):1:OPTIONS.Segment(end)-10];
        NoiseCov = [];
        SourceCov = [];
        Projector = Visu.Data.Projector;
        SourceLoc = [];
        SourceOrder = [];
        SourceOrientation = cell(1);
        ModelGain = [];
        IndepTopo = [];
        TimeSeries = [];
        PatchNdx = [];
        PatchAmp = [];
        ImageGridTime =  Visu.Data.Time(Time);
        Function = 'sptfa'; % name of this calling routine
        Fsynth = [];%G*ImageGridAmp;

        Results = struct('Comment',Comment,'Function',mfilename,'StudySubject',StudySubject,...
            'GUI',GUI,'Date','','Subject',Subject,'Study',Study,'Time',Time,...
            'ChannelFlag',Visu.Data.ChannelFlag,'NoiseCov',NoiseCov,'SourceCov',SourceCov,...
            'Projector',Projector,'SourceLoc',SourceLoc,'SourceOrder',SourceOrder,...
            'SourceOrientation',[],'ModelGain',ModelGain,...
            'IndepTopo',IndepTopo,'TimeSeries',TimeSeries,'PatchNdx',PatchNdx,...
            'PatchAmp',PatchAmp,'ImageGridAmp',ImageGridAmp,'ImageGridTime',ImageGridTime,'Fsynth',Fsynth);
        Results.SourceOrientation = SourceOrientation; % struct command could kron out

        ChannelFlag = Visu.Data.ChannelFlag;

        if(isfield(GUI,'Results')), %user provided such a string
            if(~isempty(OPTIONS.ResultFile)), % and it's not empty
                cd(Users.STUDIES)
                save(OPTIONS.ResultFile,'Comment','Function','StudySubject','GUI','Date','Subject','Study','Time',...
                    'ChannelFlag','NoiseCov','SourceCov','Projector','SourceLoc','SourceOrder',...
                    'SourceOrientation','ModelGain','IndepTopo','TimeSeries','PatchNdx',...
                    'PatchAmp','ImageGridAmp','ImageGridTime','Fsynth');
                msgbox(sprintf('Wrote results to %s',OPTIONS.ResultFile),'Notice','modal');
            end
        end



        set(handles.Tstatustxt,'String',sprintf('Status:\n\nReady'))

        %---------------------------------------------------------------------

    case 'Recursive Min. Norm.'
        OPTIONS.Lamda = str2num(get(handles.Lamda,'String'));
        errpot = .1;

        % Load Gain Matrix
        StudySubject = find_brainstorm_structure(Users.CurrentData.StudyFile);
        [path,file] = fileparts(StudySubject.Study);
        cd(fullfile(Users.STUDIES,path));

        if(~isempty(StudySubject.Subject)),
            Subject = load(fullfile(Users.SUBJECTS,StudySubject.Subject)); %load all of the information on the subject
        else
            Subject = [];
        end
        Study = load_brainstorm_file(StudySubject.Study); % and the study information

        load(StudySubject.HeadModel,'ImageGain');
        iCortex  = get(handles.LHeadModels,'Userdata');

        G = read_gain(ImageGain{iCortex(get(handles.LHeadModels,'Value'))}{1},GoodChannel,[]);

        [ImageGridAmp,lamda,res] = bst_mnef(G,Visu.Data.F(GoodChannel,[OPTIONS.Segment(1):OPTIONS.Segment(end)]),OPTIONS.Lamda,errpot);

        % now map to results

        Comment = sprintf('FOCAL MIN NORM, rank %.0f',OPTIONS.Rank);
        Date = datestr(datenum(now),0);
        Subject = Subject;
        Study = Study;
        Time = [OPTIONS.Segment(1):OPTIONS.Segment(end)];
        NoiseCov = [];
        SourceCov = [];
        Projector = Visu.Data.Projector;
        SourceLoc = [];
        SourceOrder = [];
        SourceOrientation = cell(1);
        ModelGain = [];
        IndepTopo = [];
        TimeSeries = [];
        PatchNdx = [];
        PatchAmp = [];
        ImageGridTime = Visu.Data.Time(Time);
        Function = 'focuss'; % name of this calling routine
        Fsynth = G*ImageGridAmp;

        Results = struct('Comment',Comment,'Function',mfilename,'StudySubject',StudySubject,...
            'GUI',GUI,'Date','','Subject',Subject,'Study',Study,'Time',Time,...
            'ChannelFlag',Visu.Data.ChannelFlag,'NoiseCov',NoiseCov,'SourceCov',SourceCov,...
            'Projector',Projector,'SourceLoc',SourceLoc,'SourceOrder',SourceOrder,...
            'SourceOrientation',[],'ModelGain',ModelGain,...
            'IndepTopo',IndepTopo,'TimeSeries',TimeSeries,'PatchNdx',PatchNdx,...
            'PatchAmp',PatchAmp,'ImageGridAmp',ImageGridAmp,'ImageGridTime',ImageGridTime,'Fsynth',Fsynth);
        Results.SourceOrientation = SourceOrientation; % struct command could kron out

        ChannelFlag = Visu.Data.ChannelFlag;

        if(isfield(GUI,'Results')), %user provided such a string
            if(~isempty(OPTIONS.ResultFile)), % and it's not empty
                cd(Users.STUDIES)
                save(OPTIONS.ResultFile,'Comment','Function','StudySubject','GUI','Date','Subject','Study','Time',...
                    'ChannelFlag','NoiseCov','SourceCov','Projector','SourceLoc','SourceOrder',...
                    'SourceOrientation','ModelGain','IndepTopo','TimeSeries','PatchNdx',...
                    'PatchAmp','ImageGridAmp','ImageGridTime','Fsynth');
                msgbox(sprintf('Wrote results to %s',OPTIONS.ResultFile),'Notice','modal');
            end
        end

        %-----------------------------------------------------------------------------------------------------


    case 'ST-MAP'
        OPTIONS.LamdaTime = str2num(get(handles.LamdaTime,'String'));
        if isempty(OPTIONS.LamdaTime)
            OPTIONS.LamdaTime = 0;
        end
        OPTIONS.Lamda = str2num(get(handles.Lamda,'String'));
        errpot = .1;

        % Load Gain Matrix
        StudySubject = find_brainstorm_structure(Users.CurrentData.StudyFile);
        [path,file] = fileparts(StudySubject.Study);
        cd(fullfile(Users.STUDIES,path));

        if(~isempty(StudySubject.Subject)),
            Subject = load(fullfile(Users.SUBJECTS,StudySubject.Subject)); %load all of the information on the subject
        else
            Subject = [];
        end
        Study = load_brainstorm_file(StudySubject.Study); % and the study information

        load(StudySubject.HeadModel,'ImageGain');
        iCortex  = get(handles.LHeadModels,'Userdata');

        G = read_gain(ImageGain{iCortex(get(handles.LHeadModels,'Value'))}{1},GoodChannel,[]);

        %Neighborhood System File
        cd(Users.SUBJECTS)
        NSFile = strrep(Visu.Tesselation,'.mat','_NeighSys.mat');
        if exist(NSFile,'file')
            load(NSFile)
        else
            bst_sourceimaging('CorticalNeighborhoodSystem')
        end

        load(Visu.Tesselation,'Vertices');
        Vertices = Vertices{iCortex}';

        method = 'mn' ; % 'mn' or 'mg' minnnom or minnorm gradient
        [ImageGridAmp, lamda, residuals] = ...
            bst_STMAP(G,Visu.Data.F(GoodChannel,[OPTIONS.Segment(1):OPTIONS.Segment(end)]),Vertices,NSFile,OPTIONS.Lamda,OPTIONS.LamdaTime,method,errpot);

    case 'RAP-MUSIC'
        bst_sourceimaging SubspaceSelection;
        data_selector_cb('rap-music')
end
