function varargout = headmodeler_gui(varargin)
%HEADMODELER_GUI - Application M-file to build head models
% function varargout = headmodeler_gui(varargin)
%    fig = HEADMODELER_GUI launch headmodeler_gui GUI.
%    HEADMODELER_GUI('callback_name', ...) invoke the named callback.

%<autobegin> ---------------------- 08-Jun-2004 15:11:22 -----------------------
% --------- Automatically Generated Comments Block Using AUTO_COMMENTS ---------
%
% CATEGORY: GUI and Related
%
% Alphabetical list of external functions (non-Matlab):
%   toolbox\bst_color_scheme.m
%   toolbox\bst_headmodeler.m
%   toolbox\bst_layout.m
%   toolbox\bst_message_window.m
%   toolbox\bst_win_manager.m
%   toolbox\colnorm.m
%   toolbox\data_manager.m
%   toolbox\find_brainstorm_files.m
%   toolbox\get_channel.m
%   toolbox\get_user_directory.m
%   toolbox\good_channel.m
%   toolbox\headmodeler_gui.m  NOTE: Routine calls itself explicitly
%   toolbox\mutincomp.m
%   toolbox\norlig.m
%   toolbox\view_surface.m
%
% Subfunctions in this file, in order of occurrence in file:
%   varargout = LoadStudy(studyfile) % Run when figure is created - load all Study/Subject information
%   varargout = ComputeSphereParameters_Callback(fig, eventdata, handles, varargin)
%   varargout = MEGForwardMethod_Callback(h, eventdata, handles, varargin)
%   varargout = EEGForwardMethod_Callback(h, eventdata, handles, varargin)
%   varargout = crtxlistboxx_Callback(h, eventdata, handles, varargin)
%   varargout = listboxx_Callback(h, eventdata, handles, varargin)
%   varargout = center_Callback(h, eventdata, handles, varargin)
%   varargout = radii_Callback(h, eventdata, handles, varargin)
%   varargout = conductivity_Callback(h, eventdata, handles, varargin)
%   varargout = SPACING_Callback(h, eventdata, handles, varargin)
%   varargout = ComputeVolumeSourceGrid_Callback(h, eventdata, handles, varargin)
%   varargout = verbose_Callback(h, eventdata, handles, varargin)
%   varargout = HeadModelName_Callback(h, eventdata, handles, varargin)
%   varargout = ref_Callback(h, eventdata, handles, varargin)
%   BEMParamSelectCancel_Callback(h, eventdata, handles, varargin)
%   varargout = HRunHeadModeling_Callback(h, eventdata, handles, varargin)
%   varargout = mutincomp_test(h,eventdata,handles)
%   varargout = mutincomp_basis(h,eventdata,handles)
%   varargout = estimate_memory_load(h,eventdata,handles)
%   varargout = SelectMEG_Callback(h, eventdata, handles, varargin)
%   varargout = SelectEEG_Callback(h, eventdata, handles, varargin)
%   varargout = LoadHeadmodelFile_Callback(h, eventdata, handles, varargin)
%   varargout = HEnvelopeSelect_Callback(h, eventdata, handles, varargin)
%   varargout = HEnvelopeRemove_Callback(h, eventdata, handles, varargin)
%   varargout = pViewScalp_Callback(h, eventdata, handles, varargin)
%   varargout = pViewBrain_Callback(h, eventdata, handles, varargin)
%   varargout = ViewBEMEnvelopes(h, eventdata, handles, varargin)
%   varargout = ViewSensors(h, eventdata, handles, varargin)
%   varargout = QUIT_Callback(h, eventdata, handles, varargin)
%
% Application data and their calls in this file:
%   'BEMenvelopes'
%   'BrainStormDataBase'
%   'OPTIONS'
%   'TileType'
%   
%   setappdata(BEMParamSelect,'TileType','T');
%   setappdata(TASKBAR,'BrainStormDataBase',Users)
%   setappdata(fig,'OPTIONS',OPTIONS)
%   setappdata(fig,'OPTIONS',OPTIONS);
%   setappdata(fig,'TileType','T');
%   setappdata(hFig.Head,'BEMenvelopes',BEMenvelopes)
%   setappdata(hf,'TileType','T');
%   setappdata(hh,'TileType','T');
%   
%   BEMenvelopes = getappdata(hFig.Head,'BEMenvelopes');
%   OPTIONS = getappdata(HMfig,'OPTIONS');
%   OPTIONS = getappdata(fig,'OPTIONS');
%   OPTIONS = getappdata(hFig.Head,'OPTIONS');
%   OPTIONS.BEM.EnvelopeNames = getappdata(fig,'BEMenvelopes');
%   STUDIESfolder = getappdata(getappdata(0,'BrainStormTaskbar'),'BrainStormDataBase');
%   TASKBAR = getappdata(0,'BrainStormTaskbar');
%
% Figure Files opened by this function:
%   'bem_param_select.fig'
%   mfilename
%
%   Format of strings below: Type:Style:Tag, "String", CallBack Type and Call
%   <automatic> callback is <Tag>_Callback by Matlab default
%
% Callbacks by figure bem_param_select.fig
%   uicontrol:checkbox:collocation "Collocation" uses Callback for
%     headmodeler_gui('mutincomp_test',gcbo,[],guidata(gcbo));
%   uicontrol:checkbox:constant "Constant" uses Callback for
%     headmodeler_gui('mutincomp_basis',gcbo,[],guidata(gcbo));
%   uicontrol:checkbox:galerkin "Galerkin" uses Callback for
%     headmodeler_gui('mutincomp_test',gcbo,[],guidata(gcbo));
%   uicontrol:checkbox:linear "Linear" uses Callback for
%     headmodeler_gui('mutincomp_basis',gcbo,[],guidata(gcbo));
%   uicontrol:edit:nvertmax "1000" uses Callback for
%     headmodeler_gui('estimate_memory_load',gcbo,[],guihandles(gcbf));
%   uicontrol:listbox:available "" uses Callback for
%     headmodeler_gui('HEnvelopeSelect_Callback',gcbo,[],guidata(gcbo))
%   uicontrol:listbox:removed "" uses Callback for
%     headmodeler_gui('HEnvelopeRemove_Callback',gcbo,[],guidata(gcbo))
%   uicontrol:pushbutton:pushbutton2 ">>" uses Callback for
%     headmodeler_gui('HEnvelopeSelect_Callback',gcbo,[],guidata(gcbo))
%   uicontrol:pushbutton:pushbutton3 "<<" uses Callback for
%     headmodeler_gui('HEnvelopeRemove_Callback',gcbo,[],guidata(gcbo))
%   uicontrol:pushbutton:pViewBEMEnvelopes "View Envelopes" uses Callback for
%     headmodeler_gui('ViewBEMEnvelopes',gcbo,[],guidata(gcbo))
%   uicontrol:pushbutton:pViewSensors "View Sensors" uses Callback for
%     headmodeler_gui('ViewSensors',gcbo,[],guidata(gcbo))
%   uicontrol:pushbutton:QUIT "Cancel" uses Callback for
%     headmodeler_gui('BEMParamSelectCancel_Callback',gcbo,[],guidata(gcbo))
%   uicontrol:text:Btext10 "Ordered Selection" uses Callback for
%     bem_envelope_select('StaticText1_Callback',gcbo,[],guidata(gcbo))
%   uicontrol:text:Btext12 "Available" uses Callback for
%     bem_envelope_select('StaticText1_Callback',gcbo,[],guidata(gcbo))
%   uicontrol:text:Dtext11 "Select the envelopes for BEM modelling: " uses Callback for
%     bem_envelope_select('StaticText1_Callback',gcbo,[],guidata(gcbo))
%
% Callbacks by figure headmodeler_gui.fig
%   figure::headmodeler "" uses CreateFcn for movegui center;
%   uicontrol:checkbox:ComputeVolumeSourceGrid "Compute" uses Callback for <automatic>
%   uicontrol:checkbox:EEG "EEG" uses Callback for
%     headmodeler_gui('SelectEEG_Callback',gcbo,[],guidata(gcbo))
%   uicontrol:checkbox:MEG "MEG" uses Callback for
%     headmodeler_gui('SelectMEG_Callback',gcbo,[],guidata(gcbo))
%   uicontrol:checkbox:Tverbose "Verbose" uses Callback for
%     headmodeler_gui('verbose_Callback',gcbo,[],guidata(gcbo))
%   uicontrol:edit:center "" uses Callback for <automatic>
%   uicontrol:edit:conductivity "" uses Callback for <automatic>
%   uicontrol:edit:HeadModelName "Default" uses Callback for <automatic>
%   uicontrol:edit:radii "" uses Callback for <automatic>
%   uicontrol:edit:spacing "2" uses Callback for
%     headmodeler_gui('SPACING_Callback',gcbo,[],guidata(gcbo))
%   uicontrol:popupmenu:crtxlistboxx "No Study in Memory" uses Callback for <automatic>
%   uicontrol:popupmenu:EEGForwardMethod "Single Sphere" uses Callback for <automatic>
%   uicontrol:popupmenu:listboxx "No Study in Memory" uses Callback for <automatic>
%   uicontrol:popupmenu:LoadHeadModelFile "<Load HeadModel File>" uses Callback for
%     headmodeler_gui('LoadHeadmodelFile_Callback',gcbo,[],guidata(gcbo))
%   uicontrol:popupmenu:MEGForwardMethod "Single Sphere" uses Callback for <automatic>
%   uicontrol:popupmenu:ref "." uses Callback for <automatic>
%   uicontrol:pushbutton:HComputeSphereParameters "Adjust" uses Callback for
%     headmodeler_gui('ComputeSphereParameters_Callback',gcbf,[],guidata(gcbf))
%   uicontrol:pushbutton:HRunHeadModeling "Run" uses Callback for <automatic>
%   uicontrol:pushbutton:pViewBrain "View" uses Callback for <automatic>
%   uicontrol:pushbutton:pViewScalp "View" uses Callback for <automatic>
%   uicontrol:pushbutton:QUIT "Close" uses Callback for <automatic>
%
% At Check-in: $Author: Mosher $  $Revision: 40 $  $Date: 6/08/04 1:38p $
%
% This software is part of BrainStorm Toolbox Version 2.0 (Alpha) 28-May-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> ------------------------ 08-Jun-2004 15:11:22 -----------------------


% Last Modified by GUIDE v2.0 11-May-2004 11:50:38

% /---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                                    |
% \------------------------------------------------------/
%  
% Date of creation: May 2002

% Script History -----------------------------------------
% SB  07-Aug-2002 Corrected a few typos in script comments
% SB  05-Feb-2003 Added visualization of adjusted spheres w.r.t. scalp surface or 
%                 sensor locations when clicking on the "Adjust" sphere parameters button.
% SB  03-Jun-2003 Changed GUI layout and display for compliance to new BST tiling of GUIs.
% JCM 08-Sep-2003 Commenting
%---------------------------------------------------------

Users = get_user_directory;

if nargin == 0  % LAUNCH GUI

    % Clean up previous HeadModel GUIs
    tmp = bst_win_manager;
    clear tmp.Head; clear tmp

    % Generate a structure of handles to pass to callbacks, and store it. 
    fig = openfig(mfilename,'reuse');
    bst_color_scheme(fig); 
    
    setappdata(fig,'TileType','T');
    bst_layout('align',fig,1,2,1);
   
    bst_win_manager(fig,'Head');
    
    set(fig,'Visible','on')
    
    handles = guihandles(fig);
    guidata(fig, handles);% Mathworks convention to store handles 
    
    % Code names for MEG and EEG forwards approaches as requested inputs to BST_HEADMODELER.M
    OPTIONS = getappdata(fig,'OPTIONS');
    
    OPTIONS.MEGMethods = {'meg_sphere','meg_os','meg_bem'};
    OPTIONS.MEGMethodsLabel = get(handles.MEGForwardMethod,'String'); % Full name of corresponding forward methods
    
    OPTIONS.EEGMethods = {'eeg_sphere','eeg_3sphere','eeg_3sphereBerg','eeg_os','eeg_bem'};
    OPTIONS.EEGMethodsLabel = get(handles.EEGForwardMethod,'String'); % Full name of corresponding forward methods
    OPTIONS.Verbose = get(handles.Tverbose,'value');
    
    setappdata(fig,'OPTIONS',OPTIONS)
    
    if nargout > 0
        varargout{1} = fig;
    end
    
 elseif ischar(varargin{1}) % INVOKE NAMED SUBFUNCTION OR CALLBACK
    
    try
       if (nargout)
          [varargout{1:nargout}] = feval(varargin{:}); % FEVAL switchyard
       else
          feval(varargin{:}); % FEVAL switchyard
       end
    catch
       disp(lasterr);
    end
    
 end


% --------------------------------------------------------------------
function varargout = LoadStudy(studyfile) % Run when figure is created - load all Study/Subject information

Users = get_user_directory;
Users.CurrentData.StudyFile = studyfile; % Update Current Data information
TASKBAR = getappdata(0,'BrainStormTaskbar');
setappdata(TASKBAR,'BrainStormDataBase',Users)

% Open GUI figure
fig = headmodeler_gui;
%fig = openfig(mfilename,'reuse'); 
handles = guihandles(fig);
%guidata(fig, handles);

OPTIONS = getappdata(fig,'OPTIONS'); % Start filling out the OPTIONS fields which will pass all information to the forward routines
OPTIONS.StudyFile = studyfile;

bst_message_window(...
    'Loading Study / Subject Information...');

% Gather selected study information
DMfig = bst_win_manager;
DMfig = DMfig.DataManager; % Handle to DataManager GUI.
DMhandles = guihandles(DMfig); % Get its handles
%
% FileList = getappdata(DMfig,'FileList');
% 
% StudySelection  = get(DMhandles.listbox_Study,'Value'); % Which study was selected ?
% studyfile = FileList(StudySelection).name;%Users.CurrentData.StudyFile; 
% 
cd(Users.STUDIES)
load(studyfile)
[OPTIONS.studypath,studyfile,ext] = fileparts(studyfile); 
studyfile = [studyfile,ext];

if ~isstr(Name);
   Name = 'Default';
end
StudyName = deblank(Name);
set(handles.Dstudy,'String',StudyName)
OPTIONS.Subject = BrainStormSubject;

rooot = findstr(studyfile,'brainstormstudy.mat');
if isempty(rooot)
    errordlg('Study file name should be of the form ''*brainstormstudy.mat''')
    return
end
rooot = strrep(studyfile,'.mat','');
OPTIONS.rooot = rooot;

cd(OPTIONS.studypath)

[Channel, ChannelFile] = get_channel(fullfile(Users.STUDIES,OPTIONS.studypath,studyfile));
% Store channel information for call to bst_headmodeler
OPTIONS.ChannelFile = fullfile(fullfile(Users.STUDIES,OPTIONS.studypath),ChannelFile);
OPTIONS.Channel = Channel;

if isempty(ChannelFile)
    errordlg(sprintf('Channel file %s is missing. Please have it available it the current study folder.',OPTIONS.ChannelFile), 'Missing Channel File')
    return
end

% channelfile = strrep(studyfile,'brainstormstudy','channel');
% 
% if ~exist(channelfile,'file')
%     h = errordlg(['Channel file ',channelfile,' does not exist; Please selet the proper channel file or import the data again']);
%     uiwait(h)
%     [channelfile,studypath] = uigetfile(['*_channel.mat'], 'Please Select a BrainStorm Channel File');
%     if channelfile == 0, return, end
% end
% 
% load(channelfile)

OPTIONS.EEGndx = good_channel(Channel,[],'EEG');
OPTIONS.MEGndx = good_channel(Channel,[],'MEG');

% find EEG Reference Channel
if ~isempty(OPTIONS.EEGndx)
    RefNdx = good_channel(Channel,[],'EEG Ref'); % Find reference EEG channel in Channel structure
    if isempty(RefNdx) & ~strcmp(Channel(OPTIONS.EEGndx(1)).Comment,'AVERAGE REF') % EEG ref was not defined
        bst_message_window({...
                '',...
                'Attention:',...
                sprintf('EEG reference was not defined - please edit channel file %s',OPTIONS.ChannelFile),...
                'EEG reference was assumed to be ''AVERAGE'''...
            })

        OPTIONS.EEGRef = 'AVERAGE REF';
        
    else % Specific reference channel was specified
        if ~isempty(RefNdx)
            OPTIONS.EEGRef = Channel(RefNdx).Name; 
        end
    end
    
    % Update handles.ref display
    if ~isempty(RefNdx)
        EEG_ChanList  = [{'AVERAGE REF'},{Channel(RefNdx).Name}]; % SB 19-Nov-2002 / now propose only the original EEG reference or AVERAGE reference.
        % Other electrodes are not supposed to be taken as EEG references.
        set(handles.ref,'String',EEG_ChanList,'Value',2)
    else
        EEG_ChanList  = [{'AVERAGE REF'}];
        set(handles.ref,'String',EEG_ChanList,'Value',1);
    end
    
end

    

if get(handles.MEG,'Value') & isempty(OPTIONS.MEGndx) & get(handles.EEG,'Value')==0% User wants to compute MEG but no MEG data is available
    set(handles.MEG,'Value',0)
    msgbox('Sorry - No MEG data is available')
end

if get(handles.EEG,'Value') & isempty(OPTIONS.EEGndx) & get(handles.MEG,'Value')==0 % User wants to compute EEG but no EEG data is available
    set(handles.EEG,'Value',0)
    msgbox('Sorry - No EEG data is available')
end

if isempty(OPTIONS.MEGndx)
    set(handles.MEG,'enable','off','Value',0)
else
    set(handles.MEG,'enable','on','Value',1)
end

if isempty(OPTIONS.EEGndx)   
    set(handles.EEG,'enable','off','Value',0)
else
    set(handles.EEG,'enable','on','Value',1)
end

OPTIONS.SubjectImage = '';
if ~exist(BrainStormSubject,'file') % Make access to files more robust
    tmp = fullfile(Users.SUBJECTS,BrainStormSubject);
else
    tmp = BrainStormSubject;
end

[OPTIONS.subjectpath,file,ext] = fileparts(BrainStormSubject);

if ~exist(tmp,'file')
    errordlg(...
        sprintf('Cannot find the Tessellation file: %s. Please update the Subject Information by running the Subject Editor',tmp));
    return
end
try 
    load(fullfile(Users.SUBJECTS,BrainStormSubject))
catch
    load(BrainStormSubject)
end
set(handles.Dsubject,'String',Name) % Subject's name

% Detect Tessellation files
% SB 31-May-2002: Could use what was done in DataManager but I suggest we keep
% HEADMODELER_GUI rather independent from DM for now.
% Therefore we, again, look for available tess files for this subject.

% Older approach - Tesselation field of the Subject File
% if exist('Tesselation','var') % Older BrainStorm subject - should make it obsolete in this releasae 
%     OPTIONS.TessellationFile = Tesselation;
%     if ~exist(OPTIONS.TessellationFile,'file') % Make access to files more robust 
%         OPTIONS.TessellationFile = fullfile(Users.SUBJECTS,OPTIONS.TessellationFile);
%     end
% end
% Find tessellation files for this subject
   
[TessPath, file, ext] = fileparts(BrainStormSubject);
TessPath = fullfile(Users.SUBJECTS,TessPath);
[Tess,DataPopup,Leader] = find_brainstorm_files('tess',TessPath);

set(handles.BHeadCompartments,'String', {sprintf('Head Compartments: %d available',length(Tess))})

if isempty(Tess) % No Tessellation file available for this subject
    Tess = 'N/A';
    DispComment = {''};
    OPTIONS.TessellationFile = [];
    set(handles.listboxx,'String',[{'<Best-fitting Sphere>'}]) % Default head model
    OPTIONS.TessFileComment = [];
else
    cd(TessPath)
    for k = 1:length(Tess) % Explore Tess files and read Comments
        load(Tess{k},'Comment');
        DispComment{k} = char(Comment{:});
        % OPTIONS.TessFileComment is a cell array ntess (number of tess files)
        % long. kth cell contains a nsurf (number of surfaces in kth tess file) by two array.
        % First column is the index of the tess file in OPTIONS.TessellationFile and second column 
        % is the index of the surface in the Comment field of the Tess file.        % 
        OPTIONS.TessFileComment{k} = [k*ones(length(Comment),1), [1:length(Comment)]'];
    end
    OPTIONS.TessellationFile = Tess; % Later pass this info to bst_headmodeler
    % Display available envelopes in HEADMODELER GUI
    set(handles.listboxx,'String',[{'<Best-fitting Sphere>'},DispComment])
    set(handles.crtxlistboxx,'String',[{'<None>'},DispComment])
end

% Display available head model files in popupmenu LoadHeadModelFIle
% any file in the current study folder with the keyword "headmodel"
cd(Users.STUDIES)
[HeadModelFiles,DataPopup,Leader] = find_brainstorm_files('headmodel',OPTIONS.studypath);
if isempty(DataPopup) % No headmodel file available
    HeadModelFilesString = 'No HeadModel File Available';
    set(handles.LoadHeadModelFile,'enable','off')
else
    HeadModelFilesString = {'<Available HeadModel File(s)>',HeadModelFiles{:}};
    set(handles.LoadHeadModelFile,'enable','on')
end
set(handles.LoadHeadModelFile,'String',HeadModelFilesString)

bst_message_window(...
    'Loading Study / Subject Information... -> Done')

setappdata(fig,'OPTIONS',OPTIONS)

% --------------------------------------------------------------------
function varargout = ComputeSphereParameters_Callback(fig, eventdata, handles, varargin)

fig = handles.headmodeler;

OPTIONS = getappdata(fig,'OPTIONS');
OPTIONS.Verbose = get(handles.Tverbose,'value');


% What are the selected forward methods ?
MEG = 0;
EEG = 0; % Flags
if get(handles.MEG,'Value') % MEG is requested - check associated forward method
    MEG = 1; 
    MEGMethod = OPTIONS.MEGMethods{get(handles.MEGForwardMethod,'Value')};
end
if get(handles.EEG,'Value') % EEG is requested - check associated forward method
    EEG = 1; 
    EEGMethod = OPTIONS.EEGMethods{get(handles.EEGForwardMethod,'Value')};
end

   
%%%-----------------------Compute Sphere Parameters
if OPTIONS.Verbose
    bst_message_window({...
            '----------------------------------------------------',...
            'Computing.....Best-fitting sphere parameters',...
        });
end

iScalp = get(handles.listboxx,'Value')-1; % First string of envelope listbox is always '<None>'

if iScalp == 0 % <None> (no scalp envelope was selected) - Compute sphere parameters from sensor locations
    
    % if EEG, compute sphere parameters from EEG sensor locations; otherwise use MEG's (which are however further from the scalp, but user can 
    % still alter sphere parameters manually afterwards).
    indx_channels = [];
    if EEG % Provided there are 'enough' EEG electrodes to compute a decent approximation of a spherical head... 
        if length(OPTIONS.EEGndx) < 9 
            ButtonName = questdlg(sprintf(...
                'BrainStorm has found %d sensors only to compute the sphere parameters. Please note these parameters might not be accurate. Would you like to continue ?',length(OPTIONS.EEGndx)),...
                'Few EEG sensors found', 'Quit','Continue','Continue');
            switch(ButtonName)
            case 'Continue'
            case 'Quit'
                return
            end
        end
        indx_channels = OPTIONS.EEGndx;
        dispstr = 'EEG';
        
    elseif MEG
        indx_channels = OPTIONS.MEGndx;
        dispstr = 'MEG';    
    end
    
    if isempty(indx_channels) % No sensor found
        errordlg('No sensor found - PLease make sure you have selected MEG and/or EEG in the HeadModeling window')        
        return
    end
    
        
    if OPTIONS.Verbose
        bst_message_window({...
                'No scalp surface was selected',...
                'Spherical approximation of the head will be derived from',...
                sprintf('%d %s sensor locations instead',length(indx_channels),dispstr),...
            });
    end

    % Fill a virtual Vertices array of locations for sphere-parameter computations
    Vertices = zeros(length(indx_channels),3);
    for k = indx_channels
        Vertices(k-indx_channels(1)+1,:) = OPTIONS.Channel(k).Loc(:,1)';
    end

    clear indx_channels dispstr
    
    [hh] = figure; 
    set(hh,'Name','Check Sphere Parameters','NumberTitle','off','visible','off','color','k')
    setappdata(hh,'TileType','T');
    bst_layout('align',hh,2,2,2);
    hs = scatter3(Vertices(:,1),Vertices(:,2),Vertices(:,3),'filled');
    set(hs,'MarkerFaceColor','b','MarkerSize',10)
    view(120,30)
else % Load envelope information
    
    TessFileComment = cat(1,OPTIONS.TessFileComment{:});
    iTessFile = TessFileComment(iScalp,1); % Index of tessellation file
    iGrid = TessFileComment(iScalp,2); % Index in tessellation cell arrays in this file
   
    try 
        load(OPTIONS.TessellationFile{iTessFile},'Vertices','Comment');
    
    catch
        Users = get_user_directory;
        cd(Users.SUBJECTS)
        load(fullfile(OPTIONS.subjectpath,OPTIONS.TessellationFile{iTessFile}),'Vertices','Comment');
    end
    
    Vertices = Vertices{iGrid}';
    
    % Prepare figure for display of adjusted spheres and scalp surface
    try 
        load(fullfile(OPTIONS.subjectpath,OPTIONS.TessellationFile{iTessFile}),'Faces');
    catch
        cd(Users.SUBJECTS)
        load(fullfile(OPTIONS.subjectpath,OPTIONS.TessellationFile{iTessFile}),'Faces');
    end
                
    Faces = Faces{iGrid};
    [hh,hscalp,hl] = view_surface('Check Sphere Parameters',Faces,Vertices);

    setappdata(hh,'TileType','T');
    bst_layout('align',hh,2,2,2);
    
    set(hh,'visible','off','color','k')
        
end

nscalp = size(Vertices,1); % Numbers of nodes
if nscalp > 500 % 500 points is more than enough to compute scalp's best fitting sphere
   Vertices = Vertices(unique(round(linspace(1,nscalp,500))),:); 
   nscalp = size(Vertices,1); 
end

mass = mean(Vertices); % center of mass of the scalp vertex locations   
R0 = mean(norlig(Vertices - ones(nscalp,1)*mass)); % Average distance between the center of mass and the scalp points

vec0 = [mass,R0];

set(fig,'Pointer','Watch'), drawnow
[minn,brp] = fminsearch('dist_sph',vec0,[],Vertices);
clear Vertices 
set(fig,'Pointer','Arrow'), drawnow

OPTIONS.HeadCenter = minn(1:end-1)'; % 3x1
OPTIONS.R = minn(end);

if OPTIONS.Verbose
    bst_message_window({...
            sprintf('Center of the Sphere(s) : %3.1f %3.1f %3.1f (cm)',100*OPTIONS.HeadCenter),...
            sprintf('Largest Radius : %3.1f (cm)',100*OPTIONS.R),...
            'Computing.....Best-fitting sphere parameters -> Done',...
            '----------------------------------------------------',...
        });
end

OPTIONS.Radii = [OPTIONS.R/1.14 OPTIONS.R/1.08 OPTIONS.R];
set(handles.center,'String',sprintf('%3.1f\n%3.1f\n%3.1f',OPTIONS.HeadCenter'*100)) % in cm
set(handles.radii,'String',sprintf('%3.1f\n%3.1f\n%3.1f',OPTIONS.Radii*100));
set(handles.conductivity,'String',sprintf('%3.2f\n%5.4f\n%3.2f',[.33 .0042 .33]));

setappdata(fig,'OPTIONS',OPTIONS)

figure(hh), axis off, hold on
[X,Y,Z] = sphere(20);
[TH,PHI,R] = cart2sph(X,Y,Z);
for k = 1:3
    [X,Y,Z] = sph2cart(TH,PHI,R*OPTIONS.Radii(k));
    hsph = surf(X + OPTIONS.HeadCenter(1),...
        Y + OPTIONS.HeadCenter(2),...
        Z + OPTIONS.HeadCenter(3));
    set(hsph,'Facecolor','none','edgecolor',[.8 .8 .8]/(4-k),'linewidth',2,'Tag','Sphere');
end
axis normal, axis auto
set(hh,'visible','on','Tag','SphereVisuFigure')
hold off
rotate3d on
axis vis3d, axis equal

% --------------------------------------------------------------------
function varargout = MEGForwardMethod_Callback(h, eventdata, handles, varargin)

fig = handles.headmodeler;
% Update OPTIONS.Method with the user-defined parameters 
OPTIONS = getappdata(fig,'OPTIONS');

% Forward Methods
MEG = 0; % Flag
if get(handles.MEG,'Value')
    MEG = 1;
    meg_method = OPTIONS.MEGMethods{get(handles.MEGForwardMethod,'Value')};
    OPTIONS.Method{1} = meg_method; % Enforce MEG Method to be first entry of the .Method cell
    if  get(handles.MEGForwardMethod,'Value') == 3 % Turn BEM subwindow ON
        BEMParamSelect = openfig('bem_param_select.fig','reuse'); % Handle to BEM parameter GUI
        set(BEMParamSelect,'Visible','off')
        bst_color_scheme(BEMParamSelect);
        setappdata(BEMParamSelect,'TileType','T');
        bst_layout('align',BEMParamSelect,1,2,2);
        bst_win_manager(BEMParamSelect,'Other');
        set(BEMParamSelect,'Visible','on')
        
        
        hBEM = guihandles(BEMParamSelect);
        
        % Get Envelope Names
        if iscell(OPTIONS.TessellationFile) % Several tessellation files available in subject's folder
            % Load all envelopes from every file
            for k = 1:length(OPTIONS.TessellationFile)
                Users = get_user_directory; cd(Users.SUBJECTS)
                load(fullfile(OPTIONS.subjectpath,OPTIONS.TessellationFile{k}),'Comment')
                if exist('Comment','var') % Valid tessellation file
                    for kk = 1:length(Comment) % Different surfaces from different TessFiles mught have the same name
                        % Add name of file to display
                        Comment{kk} = [strrep(Comment{kk},' ',''),' - (',strrep(OPTIONS.TessellationFile{k},'.mat',''),')'];
                    end
                    Komment{k} = char(Comment);
                    clear Comment
                end
            end
            Comment = cellstr(char(Komment));
        else
            load(OPTIONS.TessellationFile,'Comment')
        end
        
        set(hBEM.available,'String',Comment)
        set(hBEM.available,'Value',1,'Max',length(Comment))
        set(hBEM.removed,'String','','Value',1)

        % Apply Default BEM parameters values
        OPTtmp = bst_headmodeler;
        test = strcmpi(OPTtmp.BEM.Basis,'linear');
        set(hBEM.linear,'value',test);
        set(hBEM.constant,'value',~test);

        test = strcmpi(OPTtmp.BEM.Test,'galerkin');
        set(hBEM.galerkin,'value',test);
        set(hBEM.collocation,'value',~test);
        
        test = OPTtmp.BEM.ISA == 1;
        set(hBEM.isa,'Value',test);

        % Look for BEM xfer matrix files
        User = get_user_directory; 
        [DataDir,DataPopup,Leader] = find_brainstorm_files('xfer',fullfile(User.STUDIES,OPTIONS.studypath));
        test = OPTtmp.BEM.ForceXferComputation == 1 & isempty(DataDir);
        set(hBEM.ForceXferComputation,'Value',test);
        if isempty(DataDir)
            set(hBEM.ForceXferComputation,'enable','off')
        else
            set(hBEM.ForceXferComputation,'enable','on')
        end
   
        test = (OPTtmp.BEM.ForceXferComputation == 1) & isempty(DataDir) ;
        set(hBEM.ForceXferComputation,'Value',test);
        
        set(hBEM.nvertmax,'String',OPTtmp.BEM.NVertMax);
        
        
        if exist('BEMParamSelect','var')
            set(BEMParamSelect,'visible','on')
        end
    
        
    else
        % Delete BEMParameter GUI
        figs = bst_win_manager; 
        if ishandle(figs.Other)
            delete(figs.Other)
            bst_win_manager([],'Other');
        end
        
    end
    
end
setappdata(fig,'OPTIONS',OPTIONS);

if OPTIONS.Verbose
    bst_message_window({...
        'Updated',...
        sprintf('MEG Forward Method : %s',OPTIONS.MEGMethodsLabel{get(handles.MEGForwardMethod,'Value')}),...
        '------------------'});
end


% --------------------------------------------------------------------
function varargout = EEGForwardMethod_Callback(h, eventdata, handles, varargin)

fig = handles.headmodeler;
% Update OPTIONS.Method with the user-defined parameters 
OPTIONS = getappdata(fig,'OPTIONS');

if get(handles.EEG,'Value')
    EEG = 1;
    eeg_method = OPTIONS.EEGMethods{get(handles.EEGForwardMethod,'Value')};
    if isfield(OPTIONS,'Method')
        OPTIONS.Method{2} = eeg_method; % Enforce EEG Method to be second entry of the .Method cell
    else
        OPTIONS.Method{1} = [];
        OPTIONS.Method{2} = eeg_method;
    end
   
    if  get(handles.EEGForwardMethod,'Value') == 5 % Turn BEM subwindow ON
        
        BEMParamSelect = openfig('bem_param_select.fig','reuse'); % Handle to BEM parameter GUI
        set(BEMParamSelect,'Visible','off')
        bst_color_scheme(BEMParamSelect);
        setappdata(BEMParamSelect,'TileType','T');
        bst_layout('align',BEMParamSelect,1,2,2);
        bst_win_manager(BEMParamSelect,'Other');
        set(BEMParamSelect,'Visible','on')
        
        hBEM = guihandles(BEMParamSelect);
        %movegui(BEMParamSelect,'northeast')
        
        % Get Envelope Names
        if iscell(OPTIONS.TessellationFile) % Several tessellation files available in subject's folder
            % Load all envelopes from every file
            for k = 1:length(OPTIONS.TessellationFile)
                Users = get_user_directory; cd(Users.SUBJECTS)
                load(fullfile(OPTIONS.subjectpath,OPTIONS.TessellationFile{k}),'Comment')
                if exist('Comment','var') % Valid tessellation file
                    for kk = 1:length(Comment) % Different surfaces from different TessFiles mught have the same name
                        % Add name of file to display
                        Comment{kk} = [strrep(Comment{kk},' ',''),' - (',strrep(OPTIONS.TessellationFile{k},'.mat',''),')'];
                    end
                    Komment{k} = char(Comment);
                    clear Comment
                end
            end
            Comment = cellstr(char(Komment));
        else
            load(OPTIONS.TessellationFile,'Comment')
        end
        
        set(hBEM.available,'String',Comment)
        set(hBEM.available,'Value',1,'Max',length(Comment))
        set(hBEM.removed,'String','','Value',1)
        
        % Apply Default BEM parameters values
        OPTtmp = bst_headmodeler;
        test = strcmpi(OPTtmp.BEM.Basis,'linear');
        set(hBEM.linear,'value',test);
        set(hBEM.constant,'value',~test);
        
        test = strcmpi(OPTtmp.BEM.Test,'galerkin');
        set(hBEM.galerkin,'value',test);
        set(hBEM.collocation,'value',~test);
        
        test = OPTtmp.BEM.ISA == 1;
        set(hBEM.isa,'Value',test);
        
        User = get_user_directory; 
        [DataDir,DataPopup,Leader] = find_brainstorm_files('xfer',fullfile(User.STUDIES,OPTIONS.studypath));
        test = OPTtmp.BEM.ForceXferComputation == 1 & isempty(DataDir);
        set(hBEM.ForceXferComputation,'Value',test);
        if isempty(DataDir)
            set(hBEM.ForceXferComputation,'enable','off')
        else
            set(hBEM.ForceXferComputation,'enable','on')
        end
        
        set(hBEM.nvertmax,'String',OPTtmp.BEM.NVertMax);
        
        if exist('BEMParamSelect','var')
            set(BEMParamSelect,'visible','on')
        end
        
    else
        % Delete BEMParameter GUI
        figs = bst_win_manager; 
        if ishandle(figs.Other)
            delete(figs.Other)
            bst_win_manager([],'Other');
        end
    end
   
end

setappdata(fig,'OPTIONS',OPTIONS);



if OPTIONS.Verbose
    bst_message_window({...
        'Updated',...
        sprintf('EEG Forward Method : %s',OPTIONS.EEGMethodsLabel{get(handles.EEGForwardMethod,'Value')}),...
        '------------------'});
end


% --------------------------------------------------------------------
function varargout = crtxlistboxx_Callback(h, eventdata, handles, varargin)

fig = handles.headmodeler;

% Update OPTIONS.GridName with the user-defined parameters 
OPTIONS = getappdata(fig,'OPTIONS');
% ScalpSurface
iCrtx = get(handles.crtxlistboxx,'Value');
if iCrtx > 1 %  A cortical surface was selected
    EnvelopeNames = get(handles.crtxlistboxx,'String'); 
    GridName = deblank(EnvelopeNames{iCrtx}); 
    % Determine Tessellation file for this envelope
    TessFileComment = cat(1,OPTIONS.TessFileComment{:});
    iTessFile = TessFileComment(iCrtx-1,1); % Index of tessellation file
    iGrid = TessFileComment(iCrtx-1,2); % Surface index in tessellation file

    %Update these fields for later call to forward model routine
    OPTIONS.Cortex.FileName = fullfile(OPTIONS.subjectpath,OPTIONS.TessellationFile{iTessFile}); 
    OPTIONS.Cortex.iGrid = iGrid;
  
    dispstr = sprintf('      %s (from file ''%s'')',...
        GridName,strrep(OPTIONS.TessellationFile{iTessFile},'.mat',''));
  
    % Display number of vertices and faces
    Users = get_user_directory;
    load(fullfile(Users.SUBJECTS,OPTIONS.subjectpath,OPTIONS.TessellationFile{iTessFile}))
    nfaces = size(Faces{iGrid},1); clear Faces
    nvertx = size(Vertices{iGrid},2); clear Vertices
    set(handles.VerticesFaces,'String',sprintf('%d Vertices - %d Faces',nvertx, nfaces))
    
else
    
    OPTIONS = rmfield(OPTIONS,'Cortex');
    dispstr = 'None';
    set(handles.VerticesFaces,'String','')
    
end
setappdata(fig,'OPTIONS',OPTIONS);

if OPTIONS.Verbose
    bst_message_window({...
        'Updated',...
        sprintf(' ---> Cortical Image Support :'),...
        dispstr,...
        '------------------'});
end

% --------------------------------------------------------------------
function varargout = listboxx_Callback(h, eventdata, handles, varargin)

fig = handles.headmodeler;

% Update OPTIONS.GridName with the user-defined parameters 
OPTIONS = getappdata(fig,'OPTIONS');

% Clear existing sphere parameters
OPTIONS.HeadCenter = [];
OPTIONS.R = [];
OPTIONS.Radii = [];
OPTIONS.Conductivity = [];
set(handles.center,'String','') 
set(handles.radii,'String','');
set(handles.conductivity,'String','');

% ScalpSurface
iScalp = get(handles.listboxx,'Value');
if iScalp > 1 %  A scalp surface was selected
    EnvelopeNames = get(handles.listboxx,'String'); 
    ScalpSurfaceName = deblank(EnvelopeNames{iScalp});

    TessFileComment = cat(1,OPTIONS.TessFileComment{:});
    iTessFile = TessFileComment(iScalp-1,1); % Index of tessellation file
    iGrid = TessFileComment(iScalp-1,2); % Index of tessellation in the Vertices/Faces/Comments cell array
     
    dispstr = sprintf('      %s (from file ''%s'')',...
        ScalpSurfaceName,strrep(OPTIONS.TessellationFile{iTessFile},'.mat',''));
    
    %Update these fields for later call to forward model routine
    OPTIONS.Scalp.FileName = fullfile(OPTIONS.subjectpath,OPTIONS.TessellationFile{iTessFile}); 
    OPTIONS.Scalp.iGrid = iGrid;
    
    Users = get_user_directory;
    load(fullfile(Users.SUBJECTS,OPTIONS.Scalp.FileName))
    nfaces = size(Faces{iGrid},1); clear Faces
    nvertx = size(Vertices{iGrid},2); clear Vertices
    set(handles.ScalpVerticesFaces,'String',sprintf('%d Vertices - %d Faces',nvertx, nfaces))

    
else
    OPTIONS = rmfield(OPTIONS,'Scalp');
    dispstr = 'None';
end
setappdata(fig,'OPTIONS',OPTIONS);
if OPTIONS.Verbose
    bst_message_window({...
        'Updated',...
        sprintf(' ---> Scalp Envelope :'),...
        dispstr,...
        '------------------'});
end

% --------------------------------------------------------------------
function varargout = center_Callback(h, eventdata, handles, varargin)

fig = handles.headmodeler;

% Update OPTIONS.HeadCenter with the user-defined parameters 
OPTIONS = getappdata(fig,'OPTIONS');
OPTIONS.HeadCenter = str2num(get(h,'String'))/100;
setappdata(fig,'OPTIONS',OPTIONS);

if OPTIONS.Verbose
    bst_message_window({...
        'Updated',...
        sprintf('Center of the Sphere : %3.1f %3.1f %3.1f (cm)',100*OPTIONS.HeadCenter),...
        '------------------'});
end

% Visualization of spheres w.r.t sensor locations or scalp surface
hh = findobj(0,'Tag','SphereVisuFigure'); % Handle to sphere visualization  figure
if isempty(hh)
    return
end

figure(hh), hold on
% Find existing sphere surfaces
set([hh,fig],'pointer','watch')
hsphs = findobj(hh,'Tag','Sphere');
delete(hsphs)

% Draw new spheres
[X,Y,Z] = sphere(20);
[TH,PHI,R] = cart2sph(X,Y,Z);
for k = 1:3
    [X,Y,Z] = sph2cart(TH,PHI,R*OPTIONS.Radii(k));
    hsph = surf(X+OPTIONS.HeadCenter(1),Y+OPTIONS.HeadCenter(2),Z+OPTIONS.HeadCenter(3));
    set(hsph,'Facecolor','none','edgecolor',[.8 .8 .8]/(4-k),'linewidth',2,'Tag','Sphere');
end
drawnow
set(hh,'visible','on','Tag','SphereVisuFigure')
hold off
rotate3d on
set([hh,fig],'pointer','arrow')


% --------------------------------------------------------------------
function varargout = radii_Callback(h, eventdata, handles, varargin)

% Update OPTIONS.Radii with the user-defined parameters 
fig = handles.headmodeler;
OPTIONS = getappdata(fig,'OPTIONS');
OPTIONS.Radii = str2num(get(h,'String'))/100;
setappdata(fig,'OPTIONS',OPTIONS);

if OPTIONS.Verbose
    bst_message_window({...
        'Updated',...
        sprintf('Sphere Radii : %3.1f %3.1f %3.1f (cm)',100*OPTIONS.Radii),...
        '------------------'});
end

% Visualization of spheres w.r.t sensor locations or scalp surface
hh = findobj(0,'Tag','SphereVisuFigure'); % Handle to sphere visualization  figure
if isempty(hh)
    return
end

figure(hh), hold on
% Find existing sphere surfaces
set([hh,fig],'pointer','watch')
hsphs = findobj(hh,'Tag','Sphere');
delete(hsphs)

% Draw new spheres
[X,Y,Z] = sphere(20);
[TH,PHI,R] = cart2sph(X,Y,Z);
for k = 1:3
    [X,Y,Z] = sph2cart(TH,PHI,R*OPTIONS.Radii(k));
    hsph = surf(X+OPTIONS.HeadCenter(1),Y+OPTIONS.HeadCenter(2),Z+OPTIONS.HeadCenter(3));
    set(hsph,'Facecolor','none','edgecolor',[.8 .8 .8]/(4-k),'linewidth',2,'Tag','Sphere');
end
drawnow
set(hh,'visible','on','Tag','SphereVisuFigure')
hold off
rotate3d on
set([hh,fig],'pointer','arrow')


% --------------------------------------------------------------------
function varargout = conductivity_Callback(h, eventdata, handles, varargin)

% Update OPTIONS.HeadCenter with the user-defined parameters 
fig = handles.headmodeler;
OPTIONS = getappdata(fig,'OPTIONS');
OPTIONS.Conductivity = str2num(get(h,'String'));
setappdata(fig,'OPTIONS',OPTIONS);

if OPTIONS.Verbose
    bst_message_window({...
        'Updated',...
        sprintf('Tissue conductivity : %3.2f %5.4f %3.2f',OPTIONS.Conductivity),...
        '------------------'});
end


% --------------------------------------------------------------------
function varargout = SPACING_Callback(h, eventdata, handles, varargin)
% Update OPTIONS.VolumeSourceGridSpacing with the user-defined parameters 
fig = handles.headmodeler;
OPTIONS = getappdata(fig,'OPTIONS');
OPTIONS.VolumeSourceGridSpacing = str2num(get(handles.spacing,'String'))/100;
setappdata(fig,'OPTIONS',OPTIONS);

if OPTIONS.Verbose
    bst_message_window({...
        'Updated',...
        sprintf('Volume Source Grid Spacing : %3.2f cm',100*OPTIONS.VolumeSourceGridSpacing),...
    });
end

% --------------------------------------------------------------------
function varargout = ComputeVolumeSourceGrid_Callback(h, eventdata, handles, varargin)
% Update OPTIONS.VolumeSourceGrid with the user-defined parameters 
fig = handles.headmodeler;
OPTIONS = getappdata(fig,'OPTIONS');
OPTIONS.VolumeSourceGrid = get(handles.ComputeVolumeSourceGrid,'Value');
setappdata(fig,'OPTIONS',OPTIONS);

dispstr = {'Off','On'};

if OPTIONS.Verbose
    bst_message_window({...
        'Updated',...
        sprintf('Compute Volume Source Grid : %s',dispstr{OPTIONS.VolumeSourceGrid+1}),...
        '------------------'});
end

% --------------------------------------------------------------------
function varargout = verbose_Callback(h, eventdata, handles, varargin)

% Update OPTIONS.VolumeSourceGrid with the user-defined parameters 
fig = handles.headmodeler;
OPTIONS = getappdata(fig,'OPTIONS');
OPTIONS.Verbose = get(handles.Tverbose,'Value');
setappdata(fig,'OPTIONS',OPTIONS);

dispstr = {'Off','On'};

bst_message_window({...
        'Updated',...
        sprintf('Verbose mode : %s',dispstr{OPTIONS.Verbose+1}),...
        '------------------'});

% --------------------------------------------------------------------
function varargout = HeadModelName_Callback(h, eventdata, handles, varargin)

% Update OPTIONS.HeadModelName with the user-defined parameters 
fig = handles.headmodeler;
OPTIONS = getappdata(fig,'OPTIONS');
OPTIONS.HeadModelName = get(handles.HeadModelName,'String');
setappdata(fig,'OPTIONS',OPTIONS);

bst_message_window({...
        'Updated',...
        sprintf('HeadModel Name : %s',OPTIONS.HeadModelName),...
        '------------------'});


% --------------------------------------------------------------------
function varargout = ref_Callback(h, eventdata, handles, varargin)

% Update Channel File and OPTIONS.EEGRef with new value for EEG reference channel
fig = handles.headmodeler;
OPTIONS = getappdata(fig,'OPTIONS');
STUDIESfolder = getappdata(getappdata(0,'BrainStormTaskbar'),'BrainStormDataBase');
cd(fullfile(STUDIESfolder.STUDIES,OPTIONS.studypath))

load(OPTIONS.ChannelFile)

EEGRef = get(handles.ref,'Value');

if EEGRef == 1 % Average Reference
    [Channel(OPTIONS.EEGndx).Comment] = deal('AVERAGE REF');
    OPTIONS.EEGRef = 'AVERAGE REF';
    
else % Original reference electrode
    [Channel(OPTIONS.EEGndx).Comment] = deal('');
    
    % UPdate Channel structure
    RefString = get(handles.ref,'String');
    OPTIONS.EEGRef = RefString{EEGRef};
    RefInd = strmatch(RefString{EEGRef},(char(Channel(:).Name)));
    Channel(RefInd).Type = 'EEG REF'; % SB 18-Nov-2002
    
end
save(OPTIONS.ChannelFile,'Channel','-append')

bst_message_window({...
        'Updated',...
        sprintf('EEG Reference : %s', OPTIONS.EEGRef),...
        sprintf('Channel File  : %s',OPTIONS.ChannelFile),...
        '------------------'});

setappdata(fig,'OPTIONS',OPTIONS)
% --------------------------------------------------------------------
function BEMParamSelectCancel_Callback(h, eventdata, handles, varargin)
% Cancel selection of BEM parameters

% Close BEMParam window
close(gcbf)

% Set values of forward model pulldown menus to 1 
% i.e. deselect BEM
handles = bst_win_manager;
handles = guihandles(handles.Head);
fig = handles.headmodeler;
set(handles.MEGForwardMethod,'Value',1)
set(handles.EEGForwardMethod,'Value',1)
MEGForwardMethod_Callback(h, eventdata, handles, varargin)
EEGForwardMethod_Callback(h, eventdata, handles, varargin)


% --------------------------------------------------------------------
function varargout = HRunHeadModeling_Callback(h, eventdata, handles, varargin)
% Run HeadModel Calculation
fig = handles.headmodeler;
set(fig,'pointer','watch')

OPTIONS = getappdata(fig,'OPTIONS');

% Prepare call to BST_HEADMODELER.M
% Remove spurious fields
OPTIONS = rmfield(OPTIONS,{'Channel','ChannelFile','rooot'});

% Get the Methods field right
if get(handles.MEG,'Value') == 0 & get(handles.EEG,'Value') == 0 % No forward problem selected
    errordlg('Please check either of the MEG or EEG boxes.','No forward problem selected')
    return
end

if get(handles.MEG,'Value') % MEG is requested - check associated forward method
    MEGMethod = OPTIONS.MEGMethods{get(handles.MEGForwardMethod,'Value')};
    OPTIONS.Method = {MEGMethod};
end
if get(handles.EEG,'Value') % EEG is requested - check associated forward method
    EEGMethod = OPTIONS.EEGMethods{get(handles.EEGForwardMethod,'Value')};
    if get(handles.MEG,'Value')
        OPTIONS.Method = {MEGMethod,EEGMethod};
    else
        OPTIONS.Method = {EEGMethod};
    end
end

% Compute volume Source Grid for scanning techniques ?
OPTIONS.VolumeSourceGrid = get(handles.ComputeVolumeSourceGrid,'value');
OPTIONS.VolumeSourceGridSpacing = str2num(get(handles.spacing,'string'));

% Source model order 
OPTIONS.SourceModel = [];
if get(handles.CD,'Value') % CD source model requested
    OPTIONS.SourceModel(1) = -1;
end
if get(handles.CME,'Value') % CME source model requested
    OPTIONS.SourceModel(end+1) = 1;
end
if isempty(OPTIONS.SourceModel)
    errordlg('No source model was selected - please pick one.')
    return
end

%Verbose head modeling ?
OPTIONS.Verbose = get(handles.Tverbose,'Value');
OPTIONS.HeadModelFile = 'Default';
OPTIONS.ImageGridFile = 'Default';

if ~isempty(findstr([OPTIONS.Method{:}],'bem')) % Request for BEM - assign specific paramters
    % Look for BEMParameters GUI
    BEMParamSelect = findobj(0,'Type','fig','Tag','BEMParamSelect'); % Handle to BEM parameter GUI
    hBEM = guihandles(BEMParamSelect);
        
    if get(hBEM.linear,'value')
        OPTIONS.BEM.Basis = 'Linear';
    else
        OPTIONS.BEM.Basis = 'Constant';
    end
    if get(hBEM.galerkin,'value')
        OPTIONS.BEM.Test = 'Galerkin';
    else
        OPTIONS.BEM.Test = 'Collocation';
    end

    set(hBEM.galerkin,'enable','off')% prevent last-minute changes
    set(hBEM.linear,'enable','off')
    set(hBEM.ForceXferComputation,'enable','off') 
    
    OPTIONS.BEM.NVertMax = str2num(get(hBEM.nvertmax,'String')); set(hBEM.nvertmax,'enable','off')
    OPTIONS.BEM.ISA = get(hBEM.isa,'Value'); set(hBEM.isa,'enable','off') 
    OPTIONS.BEM.ForceXferComputation = get(hBEM.ForceXferComputation,'Value');
    
    %     % Now make ordered selection of envelopes
    %     bem_env_fig = bem_envelope_select('create',[],[],[],OPTIONS);
    %     uiwait(bem_env_fig);
    
    OPTIONS.BEM.EnvelopeNames = getappdata(fig,'BEMenvelopes');
    if isempty(OPTIONS.BEM.EnvelopeNames) % No BEM envelopes were selected - abort
        bst_message_window('wrap','No head tissue envelopes were selected - BEM modeling aborted')
        return
    end
    
end

%Compute HeadModel
if (OPTIONS.VolumeSourceGrid ==0) & ~isfield(OPTIONS,'Cortex')
    errordlg(...
        {'Please either check the ''Compute Volume Source Grid'' option',...
            'or select a Cortical Envelope as an imaging support.'},'No source location was specified');
    return
end
    
[G,OPTIONS2] = bst_headmodeler(OPTIONS.StudyFile,OPTIONS);

% Update parameters display
set(handles.center,'String',sprintf('%3.1f\n%3.1f\n%3.1f',OPTIONS2.HeadCenter'*100)) % in cm
set(handles.radii,'String',sprintf('%3.1f\n%3.1f\n%3.1f',OPTIONS2.Radii*100));
set(handles.conductivity,'String',sprintf('%3.2f\n%5.4f\n%3.2f',OPTIONS2.Conductivity));
set(handles.HeadModelName,'String',OPTIONS2.HeadModelName);

% Refresh DataManager display with new files available
GUIs = bst_win_manager;
data_manager('listbox_Study_Callback',GUIs.DataManager);

% Refresh headmodeler_gui display with new files available
Users = get_user_directory; 
headmodeler_gui('LoadStudy',Users.CurrentData.StudyFile)

set(fig,'pointer','arrow')

% --------------------------------------------------------------------
function varargout = mutincomp_test(h,eventdata,handles)
% Mutual incompatibility between checkboxes
handles = guihandles(gcbf);
mutincomp([handles.collocation,handles.galerkin]);

% --------------------------------------------------------------------
function varargout = mutincomp_basis(h,eventdata,handles)
handles = guihandles(gcbf);
mutincomp([handles.constant,handles.linear]);

headmodeler_gui('estimate_memory_load',gcbo,[],guihandles(gcbf));


% --------------------------------------------------------------------
function varargout = estimate_memory_load(h,eventdata,handles)
% estimate memory load of current BEM model

handles = guihandles(gcbf);

HMfig = findobj(0,'Tag','headmodeler'); if isempty(HMfig), return, end
OPTIONS = getappdata(HMfig,'OPTIONS'); % Modelling options

% Get total number of channels for computation
nchan = length([OPTIONS.MEGndx,OPTIONS.EEGndx]);

% Number of envelopes in BEM model
nenv = length(get(handles.removed,'String'));
% Now give estimation of memory load
switch(get(handles.constant,'Value'))
case 1 % Constant 
    Nfaces = 2*str2num(get(handles.nvertmax,'String'))-4;
    PhysMemLoad = (nenv*Nfaces)^2*8/1e6; % in MB / Number of Faces is the parameter 
               
    Nftot = nenv*Nfaces;
    VirtMemLoad = 3*PhysMemLoad + ... % H,L,U Matrices        
        + (3*nchan*Nftot ....         % A Matrix
        + 11*Nftot*3 ....           % Misc Workspace Matrices
        + 12*Nftot*1)*8*1.0e-6;     % Misc Workspace Matrices
    
    
    
case 0 % Linear
    
    Nverts = str2num(get(handles.nvertmax,'String'));
    Nvtot = nenv*Nverts;
    
    PhysMemLoad = (nenv*Nverts)^2*8/1e6; % in MB 
    

    VirtMemLoad = 3*PhysMemLoad ....        % H,L,U Matrices 
        + (3*nchan*Nvtot ....         % A Matrix
        + 11*Nvtot*3 ....           % Misc Workspace Matrices
        + 12*Nvtot*1)*8*1.0e-6;     % Misc Workspace Matrices
    
end

set(handles.MemoryLoad,'String',{...
        sprintf('Physical: %d MB',ceil(PhysMemLoad)),...
        sprintf('Virtual : %d MB',ceil(VirtMemLoad))...
    })

% Large EBM: warning message
if PhysMemLoad > 500 
    warndlg('This BEM geometry uses too many parameters and may trigger memory overload. Please consider reducing the number of vertices per surface.',...
        'Large geometry: possible memory limitations.')
end


% --------------------------------------------------------------------
function varargout = SelectMEG_Callback(h, eventdata, handles, varargin)
% Check whether MEG is selected for forward computation
fig = handles.headmodeler;

if get(h,'value') % MEG is selected
    set(handles.MEGForwardMethod,'enable','on') % Authorize selection of a MEG forward method
    set(handles.CME,'Value',1,'enable','on')
else
    set(handles.MEGForwardMethod,'enable','off')
    set(handles.CME,'Value',0,'enable','off')
end

% --------------------------------------------------------------------
function varargout = SelectEEG_Callback(h, eventdata, handles, varargin)
% Check whether EEG is selected for forward computation
fig = handles.headmodeler;

if get(h,'value') % EEG is selected
    set(handles.EEGForwardMethod,'enable','on') % Authorize selection of a EEG forward method
    set(handles.ref,'enable','on') % same for selection of a reference electrode
else
    set(handles.EEGForwardMethod,'enable','off')
    set(handles.ref,'enable','off') 
end

% --------------------------------------------------------------------
function varargout = LoadHeadmodelFile_Callback(h, eventdata, handles, varargin)
% Load a headmodel file available on disk for this study
fig = handles.headmodeler;
if get(handles.Tverbose,'Value')
    bst_message_window({...
            '----------------------------------------------------',...
            ' ',...
            'Load HeadModel file information',...
            ' '...
        })
end

    

OPTIONS = getappdata(fig,'OPTIONS');

% Load selected headmodel
Users = get_user_directory;
cd(fullfile(Users.STUDIES,OPTIONS.studypath))
filenum = get(h,'Value'); 
if filenum == 1, return, end % because first line of the popupmenu is just a message to user
HeadModelFile = get(h,'String');
HeadModelFile = HeadModelFile{filenum}; clear filenum

% Test for HeadModelFile version 
HeadModel_whos = whos('-file',HeadModelFile);  

DMfig = bst_win_manager;
DMfig = DMfig.DataManager; % Handle to DataManager GUI.
DMhandles = guihandles(DMfig); % Get its handles

if (isempty(strmatch('GridName',{HeadModel_whos.name})) | ...
        isempty(strmatch('HeadModelName',{HeadModel_whos.name}))) % Not valid or older headmodel file 
    
    if isempty(strmatch('Function',{HeadModel_whos.name})),
        bst_message_window({' ',...
                'Not a valid HeadModel file',...
                ' '})
    else
        set(handles.HeadModelName,'String','File is from older version of BrainStorm')
        bst_message_window({' ',...
                'File is from older version of BrainStorm',...
                ' '})
    end
else
    % GridName and HeadModelName exist, load HeadModel information from file
    HeadModel = load(HeadModelFile);
    set(handles.HeadModelName,'String',HeadModel.HeadModelName)
    
    switch HeadModel.SourceOrder
    case -1 % Current Dipole
        set(handles.CD,'Value',1)
        set(handles.CME,'Value',0)
    case 1 % CME
        set(handles.CME,'Value',1)
        set(handles.CD,'Value',0)
    end
    
    % Do the study and headmodelfile contain MEG ?
    if ~isempty(OPTIONS.MEGndx)
        if isfield(HeadModel,'MEGMethod') % 
            megFunction = strmatch(HeadModel.MEGMethod,OPTIONS.MEGMethods);
            set(handles.MEGForwardMethod,'Value',megFunction,'enable','on')
            set(handles.MEG,'Value',1)
        else
            set(handles.MEG,'Value',0)
        end
        set(handles.center,'string',num2str(100*HeadModel.Param(OPTIONS.MEGndx(1)).Center,'%3.2f')) % CHEAT - taking only Param from 1st sensor (will not be valid for overlapping sphere models)
        set(handles.radii,'string',num2str(100*HeadModel.Param(OPTIONS.MEGndx(1)).Radii',3))
        set(handles.conductivity,'string',num2str(HeadModel.Param(OPTIONS.MEGndx(1)).Conductivity',3))    
    else
        set(handles.MEG,'Value',0)
        set(handles.MEGForwardMethod,'Value',1,'enable','off')
    end
    % is EEG available in this headmodel ?
    if ~isempty(OPTIONS.EEGndx) % Yes - Identify forward method used for EEG
        if isfield(HeadModel,'EEGMethod') % 
            eegFunction = strmatch(HeadModel.EEGMethod,OPTIONS.EEGMethods);
            set(handles.EEGForwardMethod,'Value',eegFunction,'enable','on')
            set(handles.EEG,'Value',1)
        else
            set(handles.EEG,'Value',0)
        end
        set(handles.center,'string',num2str(100*HeadModel.Param(OPTIONS.EEGndx(1)).Center,'%3.2f')) % CHEAT - taking only Param from 1st sensor (will not be valid for overlapping sphere models)
        set(handles.radii,'string',num2str(100*HeadModel.Param(OPTIONS.EEGndx(1)).Radii',3))
        set(handles.conductivity,'string',num2str(HeadModel.Param(OPTIONS.EEGndx(1)).Conductivity',3))
    else
        set(handles.EEG,'Value',0)
        set(handles.EEGForwardMethod,'Value',1,'enable','off')
        set(handles.ref,'Value',1,'enable','off')
    end

    % Headmodel type: ImageGrid or SearchGrid ?
    switch(HeadModel.HeadModelType)
    case 'SearchGrid'
        % Compute spacing between source points
        diffLoc = diff(HeadModel.GridLoc{1},1,2); % Difference between consecutive grid points
        diffLoc = 100*mean(colnorm(diffLoc)); % in cm
        set(handles.spacing,'String',num2str(diffLoc,3)) 
        set(handles.ComputeVolumeSourceGrid,'Value',1)
        set(handles.crtxlistboxx,'Value',1);
        %set(handles.spacing,'enable','on')
        % Note that this spacing may not correspond with the one originally entered
        % when computing the headmodel because of the parsing step used during gridloc design based on subspace correlations
        % of corresponding gain matrices.
    case 'ImageGrid'
       % Find GridLoc support (i.e. the name of the cortical tessellation that was used as image support)
       [tmp,ImageSupport] = strtok(HeadModel.GridName{1},':');
       ImageSupport = (ImageSupport(3:end));
       % Now look in the tessellation listbox for this surface name
       crtxList = get(handles.crtxlistboxx,'String');
       icrtx = min(strmatch(ImageSupport,crtxList));
       if ~isempty(icrtx)
           set(handles.crtxlistboxx,'Value',icrtx);
           crtxlistboxx_Callback(h, eventdata, handles);
       else
           set(handles.crtxlistboxx,'Value',1,'String',{'<Cortical surface was not found>',crtxList(2:end)});
       end
       set(handles.spacing,'String','N/A') 
       set(handles.ComputeVolumeSourceGrid,'Value',0)
       %set(handles.spacing,'enable','off')
       
    end
    
end % if GridName and HeadModelName exist

if get(handles.Tverbose,'Value')
    bst_message_window({...
            '----------------------------------------------------'...
        })
end

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

function varargout = HEnvelopeSelect_Callback(h, eventdata, handles, varargin)

handles = guihandles(h);

removeID = get(handles.available,'Value');
IDs = get(handles.available,'String');
if isempty(IDs)
    return
end
chan = [1:length(IDs)];

set(handles.available,'String', IDs(setdiff(chan,removeID)));
if isempty(get(handles.removed,'String'))
    set(handles.removed,'String', [IDs(removeID)] );
else
    strtmp = ([cellstr(get(handles.removed,'String'));IDs(removeID)]);
    set(handles.removed,'String', cellstr(strtmp));
end
set(handles.available,'Value',1,'Max',length(get(handles.available,'String')))
set(handles.removed,'Value',1,'Max',max([1,length(get(handles.removed,'String'))]))

% Names of selected envelopes for BEM modeling
BEMenvelopes = get(handles.removed,'String'); % Send them to the main heamodeling GUI.

% Find the headmodeler window
hFig = bst_win_manager;
if isempty(hFig.Head) 
    errordlg('Please keep HeadModeling window open')
    return
end

for k = 1:length(BEMenvelopes) % Remove string indicating the name of the tessellation file (if any - i.e. when there are several surfaces per tessellation file)
    [BEMenvelopes{k}.TessName, BEMenvelopes{k}.TessFile] = strread(BEMenvelopes{k},'%s - (%s');
    BEMenvelopes{k}.TessFile = char(strrep(BEMenvelopes{k}.TessFile,')',''));
    BEMenvelopes{k}.TessName = char(BEMenvelopes{k}.TessName);
end

setappdata(hFig.Head,'BEMenvelopes',BEMenvelopes)

headmodeler_gui('estimate_memory_load',gcbo,[],guihandles(gcbf));

% --------------------------------------------------------------------
function varargout = HEnvelopeRemove_Callback(h, eventdata, handles, varargin)

handles = guihandles(h);

availableID = get(handles.removed,'Value');
IDs = get(handles.removed,'String');

if isempty(IDs)
    return
end

chan = [1:length(IDs)];

set(handles.removed,'Max',max([1,length(setdiff(chan,availableID))]))
set(handles.removed,'Value',1,'String',IDs(setdiff(chan,availableID)));

if isempty(get(handles.available,'String'))
    set(handles.available,'String', [IDs(availableID)]);
else
    strtmp = ([cellstr(get(handles.available,'String'));IDs(availableID)]);
    set(handles.available,'String', cellstr(strtmp));
end
set(handles.available,'Value',1,'Max',length(get(handles.available,'String')))

% Names of selected envelopes for BEM modeling
BEMenvelopes = get(handles.removed,'String'); % Send them to the main heamodeling GUI.

% Find the headmodeler window
hFig = bst_win_manager;
if isempty(hFig.Head) 
    errordlg('Please keep HeadModeling window open')
    return
end

for k = 1:length(BEMenvelopes) % Remove string indicating the name of the tessellation file (if any - i.e. when there are several surfaces per tessellation file)
    [BEMenvelopes{k}.TessName, BEMenvelopes{k}.TessFile] = strread(BEMenvelopes{k},'%s - (%s');
    BEMenvelopes{k}.TessFile = char(strrep(BEMenvelopes{k}.TessFile,')',''));
    BEMenvelopes{k}.TessName = char(BEMenvelopes{k}.TessName);
end

setappdata(hFig.Head,'BEMenvelopes',BEMenvelopes)

headmodeler_gui('estimate_memory_load',gcbo,[],guihandles(gcbf));

% --------------------------------------------------------------------
function varargout = pViewScalp_Callback(h, eventdata, handles, varargin)

% 3D rendering of scalp image support

fig = handles.headmodeler;
OPTIONS = getappdata(fig,'OPTIONS');
if ~isfield(OPTIONS,'Scalp')
    return
end

Users = get_user_directory;

load(fullfile(Users.SUBJECTS,OPTIONS.Scalp.FileName))
Faces = Faces{OPTIONS.Scalp.iGrid};
Vertices = Vertices{OPTIONS.Scalp.iGrid};
clear Comment

% Visualization of surface
figname = 'Cortical Image Support';
hcrtx = findobj(0,'Type','figure','Name',figname);
hcrtx = findobj(hcrtx,'Type','patch');
if ~isempty(hcrtx)
    set(hcrtx,'Facecolor','g','facealpha',1,'edgecolor','none');
else
    figname = 'Scalp Envelope';
end

% Visualization of surface
[hf,hs,hl] = view_surface(figname,Faces,Vertices);
set(hf,'visible','off')
setappdata(hf,'TileType','T');
bst_layout('align',hf,2,2,2);

set(hf,'color','k','visible','on')
set(hs,'Backfacelighting','unlit')

if ~isempty(hcrtx)
    set(hs,'Facecolor','r','Facealpha',.3,'edgecolor','none')
end

camlight
material dull
rotate3d on

% --------------------------------------------------------------------
function varargout = pViewBrain_Callback(h, eventdata, handles, varargin)

% 3D rendering of cortical image support

fig = handles.headmodeler;
OPTIONS = getappdata(fig,'OPTIONS');
if ~isfield(OPTIONS,'Cortex')
    return
end

Users = get_user_directory;
load(fullfile(Users.SUBJECTS,OPTIONS.Cortex.FileName))
Faces = Faces{OPTIONS.Cortex.iGrid};
Vertices = Vertices{OPTIONS.Cortex.iGrid};
clear Comment

% Visualization of surface
figname = 'Scalp Envelope';
hscalp = findobj(0,'Type','figure','Name',figname);
hscalp = findobj(hscalp,'Type','patch');
if ~isempty(hscalp)
    set(hscalp,'Facecolor','r','facealpha',.4,'edgecolor','none');
else
    figname = 'Cortical Image Support';
end

[hf,hs,hl] = view_surface(figname,Faces,Vertices);
set(hf,'visible','off')
setappdata(hf,'TileType','T');
bst_layout('align',hf,2,2,2);

set(hf,'color','k','visible','on')
set(hs,'Backfacelighting','unlit')
if ~isempty(hscalp)
    set(hs,'facecolor','g','edgecolor','none')
end

figure(hf)
camlight
material dull
rotate3d on

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

function varargout = ViewBEMEnvelopes(h, eventdata, handles, varargin)

% Find the headmodeler window
hFig = bst_win_manager;
if isempty(hFig.Head) 
    errordlg('Please keep HeadModeling window open')
    return
end
handles = guihandles(h);
BEMenvelopes = getappdata(hFig.Head,'BEMenvelopes');

if isempty(BEMenvelopes), return, end
User = get_user_directory;
OPTIONS = getappdata(hFig.Head,'OPTIONS');

SurfaceColor = {'b','g','r','y','c','m','k','w'};

for k = 1:length(BEMenvelopes)
    
    try 
        load(fullfile(User.SUBJECTS,OPTIONS.subjectpath,BEMenvelopes{k}.TessFile),'Comment');
    catch % Maybe user is using command line call to function with absolute-referenced files OPTIONS.*.TessFile 
        load(BEMenvelopes{k}.TessFile,'Comment');
    end
    
    Comment = strrep(Comment,' ','');
    % find surface in current tessellation file
    BEMenvelopes{k}.SurfId = find(strcmpi(BEMenvelopes{k}.TessName,Comment));
    IDs(k) = BEMenvelopes{k}.SurfId;
    if isempty(BEMenvelopes{k}.SurfId)
        errordlg(...
            sprintf('Surface %s was not found in file %s',...
            BEMenvelopes{k}.TessName, BEMenvelopes{k}.TessFile))
        return
    end
    
    % Load vertex locations
    try 
        tmp = load(fullfile(User.SUBJECTS,OPTIONS.subjectpath,BEMenvelopes{k}.TessFile),'Vertices');
    catch % Maybe user is using command line call to function with absolute-referenced files OPTIONS.*.TessFile 
        tmp = load(BEMenvelopes{k}.TessFile,'Vertices');
    end
    Vertices{k} = tmp.Vertices{BEMenvelopes{k}.SurfId}';
    
    % Load faces
    try 
        tmp = load(fullfile(User.SUBJECTS,OPTIONS.subjectpath,BEMenvelopes{k}.TessFile),'Faces');
    catch% Maybe user is using command line call to function with absolute-referenced files OPTIONS.*.TessFile 
        tmp = load(BEMenvelopes{k}.TessFile,'Faces');
    end
    
    Faces(k) = tmp.Faces(BEMenvelopes{k}.SurfId);
    
    % Visualization of surfaces
    for k = 1:length(Faces)
        [hf,hs(k),hl] = view_surface('Selected BEM envelopes',Faces{k},Vertices{k});
        delete(hl)
        setappdata(hf,'TileType','T');
        bst_layout('align',hf,2,2,1);
    end
    
    set(hs(k),'FaceAlpha',1/(k+1)^1.5,'edgealpha',1,'edgecolor',SurfaceColor{k},'facecolor',SurfaceColor{k})
    
end

material dull 
rotate3d on


set(handles.pViewSensors,'enable','on')


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

function varargout = ViewSensors(h, eventdata, handles, varargin)

% Find the headmodeler window
hFig = bst_win_manager;
if isempty(hFig.Head) 
    errordlg('Please keep HeadModeling window open')
    return
end
OPTIONS = getappdata(hFig.Head,'OPTIONS');
handles = guihandles(hFig.Head);

figname = 'Selected BEM envelopes';
hfig = findobj(0,'Type','figure','Name',figname);

figure(hfig)
hold on

% What are the selected forward methods ?
MEG = 0;
EEG = 0; % Flags
if get(handles.MEG,'Value') % MEG is requested - check associated forward method
    MEG = 1; 
    MEGMethod = OPTIONS.MEGMethods{get(handles.MEGForwardMethod,'Value')};
end
if get(handles.EEG,'Value') % EEG is requested - check associated forward method
    EEG = 1; 
    EEGMethod = OPTIONS.EEGMethods{get(handles.EEGForwardMethod,'Value')};
end

% if EEG, compute sphere parameters from EEG sensor locations; otherwise use MEG's (which are however further from the scalp, but user can 
% still alter sphere parameters manually afterwards).
indx_channels = [];
if EEG % Provided there are 'enough' EEG electrodes to compute a decent approximation of a spherical head... 
    indx_channels = OPTIONS.EEGndx;
    dispstr = 'EEG';
elseif MEG
    indx_channels = OPTIONS.MEGndx;
    dispstr = 'MEG';    
end

if isempty(indx_channels) % No sensor found
    errordlg('No sensor found - Please make sure you have selected MEG and/or EEG in the HeadModeling window')        
    return
end

% Fill a virtual Vertices array of locations for sphere-parameter computations
Vertices = zeros(length(indx_channels),3);
for k = indx_channels
    Vertices(k-indx_channels(1)+1,:) = OPTIONS.Channel(k).Loc(:,1)';
end

clear indx_channels dispstr

hs = scatter3(Vertices(:,1),Vertices(:,2),Vertices(:,3),'filled');
set(hs,'MarkerFaceColor','b','MarkerSize',10)



% --------------------------------------------------------------------
function varargout = QUIT_Callback(h, eventdata, handles, varargin)

delete(gcbf)
delete(findobj(0,'tag','BEMParamSelect','type','figure')) % close BEM parameters window


