function [HeadModel, Channel, GoodChannel, DataFlag] = extract_channels(HeadModelFilename,StudySubject,GUI);
%EXTRACT_CHANNELS - Extract the channels to be used in this particular run
% function [HeadModel, Channel, GoodChannel, DataFlag] = extract_channels(HeadModelFilename,StudySubject,GUI);
% Given the current HeadModelFilename,as well as the StudySubject and GUI information,
%  extract the Channels, HeadModel, and GoodChannel index to be used
% Will most likely become a local subfunction of rapmusic_gui

%<autobegin> ---------------------- 09-Jul-2004 22:17:03 -----------------------
% --------- Automatically Generated Comments Block Using AUTO_COMMENTS ---------
%
% CATEGORY: Data Processing
%
% Alphabetical list of external functions (non-Matlab):
%   toolbox\bst_message_window.m
%   toolbox\get_channel.m
%   toolbox\get_user_directory.m
%   toolbox\good_channel.m
%
% At Check-in: $Author: Mosher $  $Revision: 22 $  $Date: 7/09/04 8:42p $
%
% This software is part of BrainStorm Toolbox Version 2.0 (Alpha) 09-Jul-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> ------------------------ 09-Jul-2004 22:17:03 -----------------------


% $Date: 7/09/04 8:42p $ $Revision: 22 $

% From the revisions that Sylvain made earlier (2000-2001), 
%  it's not clear to me (JCM) that
%  a file with combined eeg and meg would be properly handled. Therefore reworking
%  the logic, with copious comments.

% ------------------------------------------------------------------------
% JCM 14-May-2002 extracted from rapmusic_gui, trying to understand logic
%                 and the adaptation to multiple head model files
% SB  20-Nov-2002 Updated channel extraction according to new management of 
%                 MEG REF and EEG REF channel types.
% JCM 20-Nov-2002 Problems with reference channels not being present
% JCM 12-Nov-2003 Calculating GoodChannel in all cases
% SB  29-Mar-2004 Changing to Users.STUDIES before calling GET_CHANNEL
% SB  09-Jul-2004 Fixed bug for CTF computation when first MEG channel containing coef
%                 correction information is declared as BAD: store coef
%                 information in first good MEG channel and trim number of
%                 lines of correction matrix to proper number of good MEG
%                 channels only.
% -------------------------------------------------------------------------



% load up the information needed for the gain matrix
Users = get_user_directory; % current database location
cd(Users.STUDIES)
Channel = get_channel;

% Load the HeadModel data
HeadModel = load(fullfile(Users.STUDIES,HeadModelFilename));


% now alter the data according to the bad channels
% Note that Data.ChannelFlag has the original channel flag information on 
%  Bad and not Bad. GUI.ChannelFlag set to one says analyze that channel

% What kind of data is this?
% Form a cell array of just the desired channels
GoodTypes = {Channel(GUI.ChannelFlag == 1).Type};
% What kind are these?
GoodTypes = unique(GoodTypes); % the unique cell array subset,e.g. {'MEG'}


% Are there special reference channel considerations?
%  See Channel.mat structure description in the ParameterDescriptions document.
% Assign the reference information of Ctf channels to the first MEG channel 


for i = 1:length(GoodTypes), % for each desired type
   
   switch GoodTypes{i}
      
   case 'EEG'
      % there are EEG channels being used.
      % We need to make sure that the EEG REF channel is included
      
      DataFlag = GoodTypes{i};
      
      % Reference Channel
      %ref = strmatch('AVERAGE REF',{Channel.Comment});
      EEGChannel = good_channel(Channel,GUI.ChannelFlag,'EEG'); % SB 20-Nov-2002 / Simplified extraction of EEG channels - without reference electrode if any.
      EEGRef = good_channel(Channel,[],'EEG REF'); % Empty if reference is AVERAGE REF
      
      allGoodChannel = [EEGChannel EEGRef];
      GoodChannel = [EEGChannel];
      
      %       if length(ref) > 2 % Average reference
      %          ref = 0; % there is no reference channel, use average instead
      %       else
      %          ref = strmatch('EEG REF',{Channel.Comment}); % which channel is it
      %          GoodChannel = sort([GoodChannel, ref]); % add it into the mix
      %       end
      %       
      
   case 'MEG'
      % The reference channels are embedded in the first MEG channel, by design.
      % But the possibility exists that we have excluded the first MEG channel
      %  from the analysis. If so, map the reference channel information into the
      %  first MEG channel actually used in this analysis
      
      DataFlag = 'MEG';
      
      GoodChannel = good_channel(Channel,GUI.ChannelFlag,'MEG'); % All MEG channels

      MEGRef = good_channel(Channel,[],'MEG REF'); % Indices of MEG reference channels.
      
      % JCM 20-Nov-2002, problem is that we don't have any MEG REF used at the moment, but
      %  we do embed CTF data into the comments of the first usable MEG channel.
      %  The following "if" statement is causing trouble, so I blanked it out for now,
      %  until SB and JCM understand the new conventions.
      % 21 Nov, reactivating with changes in data_manager
      
      if ~isempty(MEGRef),
         % there might possibly be CTF reference channel information
         %good channels indexer, scan all channels that are MEG types
         
         allMEGChannel = good_channel(Channel,[],'MEG'); % SB 20-Nov-2002 / All MEG channels
                  
         if (allMEGChannel(1) ~= GoodChannel(1)), 
            % The case we are worried about. The reference information is not in the first MEG channel
            % Critical case where the reference information should be stored in the first good MEG channel  
            
            if(~isempty(Channel(allMEGChannel(1)).Comment)), % there are reference sensors in Channel(1)
               Channel(GoodChannel(1)).Comment = Channel(allMEGChannel(1)).Comment;
               % reference information mapped
            end
         end
         
         % Trim-out lines of the correction matrix corresponding to bad MEG
         % channels
         Channel(GoodChannel(1)).Comment = Channel(GoodChannel(1)).Comment(GoodChannel- allMEGChannel(1)+1,:);
         
         allGoodChannel = [GoodChannel MEGRef]; % Add indices of MEG REF channels to GoodChannel set.
         GoodChannel = [GoodChannel];
      else
         % no MEG reference channels
         allGoodChannel = [GoodChannel];
         GoodChannel = [GoodChannel];
      end
      
      % result is that GoodChannel(1) now has a loaded RefChannel and grad3 fields
      %  that are stripped down to just the number of good channels
      
   case 'Fusion'
      bst_message_window('Cannot handle type FUSION Data yet in RAP-MUSIC')
      return
      
   end % switch if EEG or MEG or Fusion
     
end

% CHEAT: DataFlag is confused if both EEG and MEG data are in this set

% Now reduce the HeadModel and Channel to just the GoodChannels

HeadModel.Param = HeadModel.Param(GoodChannel);
HeadModel.Function = HeadModel.Function(GoodChannel);

if(isstr(HeadModel.Gain{1})), % it's a filename, load it
   fid = fopen(fullfile(Users.STUDIES,HeadModel.Gain{1}),'rb','ieee-be');
   if(~fid),
      errordlg(sprintf('Unable to open %s',HeadModel.Gain{1}));
      return
   end
   tempRows = fread(fid,1,'uint32');
   HeadModel.Gain{1} = fread(fid,[tempRows, inf],'float32');
   
   fclose(fid);
end

HeadModel.Gain{1} = HeadModel.Gain{1}(GoodChannel,:);

Channel = Channel(allGoodChannel);
 
