function F = meg4read(dsfolder,Trials,TimeIndices,CHANNELS,VERBOSE)
%MEG4READ - READ CTF MEG4 file format
% function F = meg4read(dsfolder,Trials,TimeIndices,CHANNELS,VERBOSE)
% reads directly the meg4 file in a ds folder once the res4 was read using ds2brainstorm
% useful when extracting data by block when using native file format and when dat awas not translated in
% BraInStorm data format (mainly for disk space saving purposes)
% 
% dsfolder : character string containing the name of the original ds folder
% Trials : integer specifying the trial to extract from the data file
% TimeIndices : a vector specifying the indices of the time samples to extract (default is extract entire trial length)
% CHANNELS : indices of channels to extract
% VERBOSE : turn on/off verbose mode (default is off : 0 )
%
% F : a matrix containing the data.

%<autobegin> ---------------------- 26-May-2004 11:30:56 -----------------------
% --------- Automatically Generated Comments Block Using AUTO_COMMENTS ---------
%
% CATEGORY: Utility - General
%
% Alphabetical list of external functions (non-Matlab):
%   toolbox\bst_message_window.m
%   toolbox\ds2brainstorm.m
%
% At Check-in: $Author: Mosher $  $Revision: 14 $  $Date: 5/26/04 9:59a $
%
% This software is part of BrainStorm Toolbox Version 2.0 (Alpha) 24-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> ------------------------ 26-May-2004 11:30:56 -----------------------


DCOffset = 0; % Remove DCoffset based on entire trial length

% Load Res4 information
[PATHSTR,NAME,EXT,VERSN] = fileparts(dsfolder) ;

if ~isempty(PATHSTR)
    cd(PATHSTR)
end

try 
    [path,rootname] = fileparts(dsfolder);
    if isdir('BrainStorm') % Was a BrainStorm directory created in .ds folder ? Yes if ds2brainstorm called from importdata
        res4_mat = fullfile('BrainStorm',[rootname,'_res4.mat']);
    else
        res4_mat = [rootname,'_res4.mat'];
    end
    
    try 
        load(res4_mat)
    catch
        [f,Channel,imegsens,ieegsens,iothersens,irefsens,grad_order_no,no_trials,filter,Time,RunTitle] = ds2brainstorm(dsfolder,VERBOSE,3,[],[],[],1);
        res4_mat = fullfile('BrainStorm',[rootname,'_res4.mat']);
        load(res4_mat)
    end
    
    nref=length(irefsens);
catch
    errordlg(...
        {sprintf('Please make sure that resource file %s is available',[NAME,'_res4.mat']),...
            sprintf('in %s.', dsfolder)},...
        'res4 file missing'...
        )
    return
end

if nargin == 3
    VERBOSE = 0;
    CHANNELS = [1:length(Channel)];
elseif nargin == 4
    CHANNELS = [1:length(Channel)];
end


meg4file = ([NAME,'.meg4']);

[meg,message] = fopen(meg4file,'rb','s'); % Big-endian byte ordering
if meg < 0
    errordlg(message,mfilename)
    return
end

% no_samples = round((TIME(end)-TIME(1))*gSetUp.sample_rate) +1; % Number of time samples in a trial

% Extraction time range (sample indices)
if nargin == 2 % Not defined: take entire trial length
    TimeIndices = 1:length(Time);
else
    if length(TimeIndices) > length(Time)
        errordlg(...
            sprintf('Only %d time samples are available in each trial. Please adjust time extraction window',length(Time)));
        return
    end
end

header = char(fread(meg,8,'char')');

%implicit.sample_rate = 1/(Time(2) - Time(1)); % i.e. the one used in the data file given the time begin_time end period.

no_samples = length(TimeIndices);
no_trials = 1:gSetUp.no_trials;

diff_trials = diff(no_trials)-1; % Number of trials between each selected trials (useful when skipping a few trials)

ByteSizeOfTrial= no_channels*gSetUp.no_samples*4; % Byte size of a single trial of data (Int32 coding) 

samples_skip = gSetUp.no_samples-no_samples; %Number of time samples to skip per channel

LastChannelSkip = (no_channels - max(CHANNELS))*gSetUp.no_samples*4; % Skip data from last channels

channels = [min(CHANNELS):max(CHANNELS)]; % Block of channels to extract.

FirstChannelSkip = (min(CHANNELS)-1)*gSetUp.no_samples*4; % Skip data from first channels
no_channels = length(channels);    

F = zeros(no_channels,length(Trials)*no_samples);

for k = 1:length(Trials)
    itrial = Trials(k); %0
    %itrial = itrial+1;
    if VERBOSE
        bst_message_window({...
                sprintf('ds2brainstorm -> Trial %d / %d',itrial,length(no_trials))...
            })
    end
    
    Ftmp = zeros(no_channels,no_samples);
    
    if itrial == no_trials(1) % Read first trial
        fseek(meg,(itrial-1)*ByteSizeOfTrial + FirstChannelSkip + (TimeIndices(1)-1)*4 ,-1);
    else % just shift from the size of a trial
        fseek(meg,LastChannelSkip + (itrial-1)*ByteSizeOfTrial + FirstChannelSkip + (TimeIndices(1)-1)*4 ,-1);
    end
    
    Ftmp = fread(meg,[no_samples no_channels],[num2str(no_samples),'*int32=>int32'], samples_skip*4)';
    %F = F(CHANNELS-min(CHANNELS)+1,:);
    
    % IMPORTANT NOTICE : Applying Gradient Correction
    % Data are saved in a given nth-order gradient correction
    % Applying gradient correction is only needed for forward model computation
    % or if it is desired to reverse to lower-order gradient correction (see importdata.m for instance).
    
    % Apply channel gains
    Ftmp = diag(1./gain_chan(CHANNELS))*double(Ftmp);
    
    % Removing DC Offset
    if isempty(DCOffset)
        DCOffset = 0;
    end
    switch(DCOffset)
    case {1,2}
        if VERBOSE
            bst_message_window({...
                    'ds2brainstorm -> Removing DC offset'...
                })
        end
        
        if DCOffset == 1 % Based on entire trial length
            [tmp ia ib] = intersect(CHANNELS,imegsens);
            
            if ~isempty(ia)
                Ftmp(ia,:) = Ftmp(ia,:) - repmat(mean(Ftmp(ia,:)')',1,size(Ftmp,2));
            end
            
        else % Based on pretrigger
            TimeNeg = intersect(find(Time(TimeIndices(1):TimeIndices(end)) >= TIME(1)),find(Time(TimeIndices(1):TimeIndices(end))<0)); % Find pretrigger time points
            if ~isempty(TimeNeg)
                [tmp ia ib] = intersect(CHANNELS,imegsens);
                if ~isempty(ia)
                    Ftmp(ia,:) = Ftmp(ia,:) - repmat(mean(Ftmp(ia,TimeNeg)')',1,size(Ftmp,2));
                end
                
            else
                % Do nothing
            end
        end
        
        if VERBOSE
            bst_message_window({...
                    'ds2brainstorm -> Removing DC offset - DONE'...
                })
        end
        
    otherwise
        % Do nothing
    end
    
    F(:,(itrial-1)*no_samples + TimeIndices) = Ftmp;
    
end

if VERBOSE
    bst_message_window({...
            sprintf('ds2brainstorm -> Done')...
        })
end

fclose(meg);

