function varargout=dopFindRawMarkers(eventData,sampleRates,varargin)
%% examine external channels
% check for markers - return number
% eMarkers=[3 4 5 6];

% place a row of zeros at the top to offset the markers
% take the difference between points
% find theses differences greater than some threshold
%
% [evMarkers found trimScreen]=dopFindRawMarkers(eventData,[epoch],[trim])
%
% notes:
% - the trim option uses epoch information to cut off the ends of the run
% data.  It bases this on the earliest and latest markers independent of
% channel.  That is, if there are two marker channels, data will be trimmed
% based on the earliest and latest from all the markers collectively.
%

% -------------------------------------------------------------------------
%     'dopOSCCI' summarises functional transcranial Doppler ultrasonography
%     (fTCD) data.  The primary function of the software is to assess the
%     hemispheric lateralization of cognitive function.
%
%     Copyright (C) 2011 the Chancellor, Masters and Scholars of the
%     University of Oxford.
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>
%     or write to the Free Software Foundation, Inc., 51 Franklin Street,
%     Fifth Floor, Boston, MA 02110-1301, USA.
%
%     Authors: Nic Badcock, Georgina Holt, Anneka Holden, & Dorothy Bishop.
%       contact: nicholas.badcock@gmail.com
%
%     Type 'dopOSCCIlicense' for license details or see dopOSCCIlicense.txt
%
try
    mkHt=1000; % threshold, greater than 1000 = marker
    % this may depended upon individual doppler settings - probably should be
    % an option
    foundMarkers=0;
    trim=0;
    epoch=[];
    epochTrim=ones(ceil(size(eventData,1)/floor((sampleRates(1)/sampleRates(2)))),1); % keep all data
    
    if ~isempty(varargin)
        for i=1:length(varargin)
            tv=varargin{i};
            if isnumeric(tv)
                % is it the epoch variable
                if size(tv,2)>1
                    epoch=tv;
                    
                    % is it the trim variable?
                elseif length(tv)==1 && tv<=1 && i==2
                    trim=tv;
                else
                    mkHt=tv;
                end
                clear tv
            end
        end
        %     switch length(varargin)
        %         case 2
        %             epoch=varargin{1};
        %             trim=varargin{2};
        %     end
    end
    
%     markers=find(diff(vertcat(zeros(1,size(eventData,2)),eventData(1:size(eventData,1)-1,:)))>abs(mkHt))-1;
%     % included absolute shift here for Marietta's data - need to get her to
%     % adjust it but it will actually work.
%     markers=fix(markers/(sampleRates(1)/sampleRates(2)));
    
    % set to downsample marker values
    if exist('downsample','file') == 2
        downMkData = downsample(eventData,floor(sampleRates(1)/sampleRates(2)));
    else
        downMkData = dopDownsample(eventData,floor(sampleRates(1)/sampleRates(2)));
    end
    markers = [];
    for i = 1 : size(eventData,2)
        cell_markers{i} = find(diff(vertcat(zeros(1,size(downMkData(:,i),2)),downMkData(1:size(downMkData(:,i),1)-1,i)))>abs(mkHt))-1;
        markers = vertcat(markers,cell_markers{i});
    end
    
    
    %% > create marker channel
    evCh=zeros(ceil(size(eventData,1)/floor((sampleRates(1)/sampleRates(2)))),size(eventData,2));
    if ~isempty(markers) && sum(markers)
        foundMarkers=1;
        numCh=size(evCh,2);
        for i=1:numCh
%             tmpFilter=and(markers>size(evCh,1)*(i-1),markers<size(evCh,1)*i);
%             evCh(markers(and(markers~=0,tmpFilter)))=120+(i-1)*20;
            tmpMarkers = cell_markers{i}; tmpMarkers = tmpMarkers(tmpMarkers>0);
            evCh(tmpMarkers,i)=120+(i-1)*20;
        end
    else
        fprintf('\n\t%s\n','=> No markers found...');
    end
    
    %% trim by epoch
    if ~isempty(markers) && numel(markers) > 1 && trim
        fprintf('%s\n',['Trimming data by epoch: ', num2str(epoch)]);
        dataPoints = [min(markers(markers~=0)) max(markers)]+epoch;
%         %     epochTrim=[];
% %         dataPoints=[find(evCh,1,'first') find(evCh,1,'last')]+epoch;
%         dataPoints=[find(evCh,1,'first') find(evCh,1,'last')]+epoch;
%         % but it would be possible for the final point from an earlier channel
%         % to actually be the last point in time
%         % > correct for this but searching for the 'final' point within each
%         % channel of markers:
%         % >> is anything in channel 1 longer than final selected? etc.
%         if numCh >1
%             tmpMax = find(evCh,1,'last');
%             for j = 1 : numCh-1
%                 tmpLimits = [1+length(evCh)*(j-1) 1+length(evCh)*(j)];
%                 tmpMarkers = markers(and(markers>tmpLimits(1),markers<tmpLimits(2)))-(length(evCh)*(j-1));
%                 if tmpMax(end)-(length(evCh)*(j)) < max(tmpMarkers)
%                     tmpMax(end+1) = max(tmpMarkers);
%                 end
%             end
%             dataPoints=[find(evCh,1,'first') tmpMax(end)]+epoch;
%         end
        fprintf('\tdata points =  %u %u\n', dataPoints);
        % this find method returns references independent of rows or columns so
        % need to correct for this.
        for i=1:2 % for each data point
%             j=0;
%             while dataPoints(i)>size(evCh,1) && i==1 || ... % unlikely...
%                     dataPoints(i)>size(evCh,1) && i==2 && j<numCh-1
%                 dataPoints(i)=dataPoints(i)-size(evCh,1);
%                 j=j+1; % can't correct more than the number of channels
%             end
            % correct if outside data range
            if i==1 && dataPoints(i)<1 % first can't be less than 1
                fprintf('\t\t adjusting lower trim: from %u to %u\n',dataPoints(1),1)
                dataPoints(i)=1;
            end
            if i==2 && dataPoints(i)>size(evCh,1) || i==2 && dataPoints(i)<size(evCh,1)
                % second can't be greater than available rows
                fprintf('\t\t adjusting upper trim: from %u to %u\n',dataPoints(2),size(evCh,1))
                dataPoints(i)=size(evCh,1);
            end
        end
        epochTrim(1:end)=0;
        epochTrim(dataPoints(1):dataPoints(2))=1;epochTrim=logical(epochTrim);
        fprintf('\t trim times = %3.2f %3.2f minutes (%u %u seconds)\n',dataPoints*(1/sampleRates(2))/60,fix(dataPoints*(1/sampleRates(2))));
        %     evCh=evCh(epochTrim,:);
    end
    %% set outputs:
    varargout{1}=evCh;
    varargout{2}=foundMarkers;
    varargout{3}=logical(epochTrim);
    %% clear data
    clear mkHt markers eventData foundMarkers trim epoch evCh epochTrim i j numCh
catch err
    %% catch dopOSCCI error
    save(dopOSCCIdebug(mfilename));rethrow(err);
end