trivialPeakFinderSNR(sigt, ampN, thrSNR)

function [Peak, leftMin, rightMin] = trivialPeakFinderSNR(sigt, ampN, thrSNR)
%
% The function performs peaks detection in TOF spectrum by 
% identifying local maxima above threshold SNR in respect to 
% closest minima and globally. The input spectrum is assumed to be
% already smoothed and baseline corrected.
%
% ALGORITHM BASED ON CROMWEL-TOOLBOX PEAK FINDER:
% Copyright (c) 2003, 2004 UT M.D. Anderson Cancer Center. All rights reserved.
% See the accompanying file "license.txt" for details.
%
% Original version by Kevin Coombes
% $Revision: 2.2 $
% Last modified by $Author: krc $ on $Date: 2004/03/01 21:41:15 $.
% NOTE: SNR threshold revision added by Dariya Malyarenko (CWM) on 11/07/07.
%
%
% INPUT:
%
% sigt - TOF signal vector
% ampN - noise amplitude vector for "sigt" (length of "sigt"): can be
% zeros(length(sigt),1)+constNoiseAmpt, if noise is static
%(estimated,e.g., as mean(abs(w(nopeaks)-mean(w(nopeaks))))~ 3-4 STD of noise)
% thrSNR - SNR detection threshold
%
% OUTPUT:
%
%	Peak - a list of peak locations, in ticks
%	leftMin -	a list of the closest minimum to the left of each peak
%	rightMin - a list of the closest minimum to the right of each peak
%
% USAGE:
% [Peak, leftMin, rightMin] = trivialPeakFinderSNR(sigt, ampN, thrSNR);
% Dependency: none
%

%% Step 1: Find all local maxima and associated local minima.

% variable noise amplitude added pn 06/12/09
if ne(length(sigt), length(ampN)) 
    error('noise amplitude vector should be the same length as signal');
end;

w = sigt;
x = diff(w);		% first difference
xl = x(1:end-1);	% left side
xr = x(2:end);		% right side
neartop = (xl > 0 & xr <= 0) | (xl == 0 & xr < 0);
nearbot = (xl < 0 & xr >= 0) | (xl == 0 & xr > 0);
xx = 1:length(x);
nb = [1 xx(nearbot)+1];	% tentative minimum
nt = xx(neartop)+1;		% tentative maximum

%% Step 2: Coalesce rising plateaus to find the right-most one.
curt = 1;
curb = 1;
Peak = [0];
bots = [1];
numtop = 1;
numbot = 1;
while curt < length(nt) & curb < length(nb),
	while nb(curb)<nt(curt) & curb < length(nb),
   	curb = curb+1;
	end
	bots(numbot) = nb(curb);
	numbot = numbot+1;
	while nt(curt) < nb(curb) & curt < length(nt),
	   curt = curt + 1;
	end
	Peak(numtop) = nt(curt);
   numtop = numtop + 1;
end
clear numtop numbot curt curb nb nt nearbot neartop

% After combining maxima, we then figure out where the closest minima
% are to both left and right. These might not agree for adjacent peaks
% because of the presence of plateaus.
leftMin = [1];	% list of resolvable minima to left of peak.
rightMin = [1];	% similar list to right of peak.
for i = 2:length(Peak),
   temp = Peak(i-1):Peak(i);
   leftMin(i) = max(temp(w(temp)==min(w(temp))));
   rightMin(i) = min(temp(w(temp)==min(w(temp))));
end
rightMin = [rightMin(2:end) length(x)];
% SNR threshold revision added on 11/07/07.

kim=find(abs(w(leftMin)-w(Peak))./ampN(Peak)>=thrSNR | (abs(w(rightMin)-w(Peak)))./ampN(Peak)>=thrSNR) ;

Peak=Peak(kim);
Peak=Peak(find((w(Peak)./ampN(Peak))>=thrSNR));

return;