
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% General description:
%
% This function performs baseline corrrection by flatting out non-peack
% spectrum
%
% input:
%       Peak - peak locations
%       y - spectrum
%
% output: 
%       yy - spectrum after baseline correction
%
% Written by: Wei Wang
% Machine Learning for Signal Processing Lab
% University of Maryland, Baltimore County
% 8/21/2006
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function yy=baseline_correction(Peak, y)

scalexy=length(y)/(max(y)-min(y)); %  y value is multiplied by scalexy when 
                                   %  calculating angle related quantities
scalexy=scalexy/1.3;  % 1.3 is an approximted ratio figure scope
y=y*scalexy;

%% Step 1: Find all local maxima and associated minimums
x = diff(y);		% 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);
nt = xx(neartop)+1;		% tentative maximum

% After combining maxima, we then figure out where the closest minima
% are to both left and right. 
leftMinAll = [1];	% list of resolvable minima to left of peak.
rightMinAll = [1];	% similar list to right of peak.
for i = 2:length(nt),
   temp = nt(i-1):nt(i);
   leftMinAll(i) = max(temp(y(temp)==min(y(temp))));
   rightMinAll(i) = min(temp(y(temp)==min(y(temp))));
end
tt=[nt(length(nt)):length(y)];
rmin=min(tt(y(tt)==min(y(tt))));
tt=[1:nt(1)];
lmin=max(tt(y(tt)==min(y(tt))));
leftMinAll(1)=lmin;
rightMinAll = [rightMinAll(2:end) rmin];

%% Step 2: Find neighbour mininums for chosen peaks
lenp=length(Peak);
for i = 1:lenp
   leftMin(i) = leftMinAll(nt==Peak(i));
   rightMin(i) = rightMinAll(nt==Peak(i));   
end

plotflag=0;
%% display in figure
if plotflag==1    
figure;
plot(Peak,y(Peak),'o',leftMin,y(leftMin),'d',rightMin,y(rightMin),'*');
hold on;plot(y);
end

%% Step 2: Combine two peaks when their baseline angle is less than a threshold
tcosangle=cos(120*pi/180); % angle has to larger than 120 to be separate peaks
peakflag=ones(1,lenp);
tolpks=length(y)/20;
for i = 1:lenp-1
    if (Peak(i+1)-Peak(i))> tolpks
        continue;
    end
    
    x1=[leftMin(i) y(leftMin(i))];
    x2=[rightMin(i) y(rightMin(i))];
    x3=[leftMin(i+1) y(leftMin(i+1))];
    x4=[rightMin(i+1) y(rightMin(i+1))];
    
    u=x2-x1;
    v=x4-x3;
    crossproduct=u(1)*v(2)-u(2)*v(1);
    if crossproduct>0
        continue;
    end
    
    l1=norm(x3-x4);
    l4=norm(x1-x2);
    l3=norm(x1-x2+x3-x4);
    cosangle=(l1^2+l4^2-l3^2)/(2*l1*l4);
%    plot([leftMin(i) rightMin(i)],[y(leftMin(i)) y(rightMin(i))],'y');
   % http://mathworld.wolfram.com/Line-LineAngle.html
%    cosangle=(x2-x1)*(x4-x3)'/(norm(x2-x1)*norm(x4-x3));    
    
    if cosangle>tcosangle 
        peakflag(i+1)=0; % right peak merge to left
    end
end
count=0;
% reset neibough minimums information
for i = 1:lenp
   if peakflag(i)==1
       count=count+1;       
       leftMinNew(count)=leftMin(i);
       rightMinNew(count)=rightMin(i);
       peakNew(count)=Peak(i);
   else
       rightMinNew(count)=rightMin(i);
       peakNew(count)=Peak(i);
   end
end
leftMin=leftMinNew;
rightMin=rightMinNew;

plotflag=0;
%% display in figure
if plotflag==1    
figure;
plot(peakNew,y(peakNew),'o',leftMin,y(leftMin),'d',rightMin,y(rightMin),'*');
hold on;plot(y);
for i=1:lenp
    plot([leftMin(i) rightMin(i)],[y(leftMin(i)) y(rightMin(i))],'r','Linewidth',2);
end
end

%% Step 3: Find the left and right minimums directly connected below peak
lenp=length(leftMin);
for i=1:lenp
   lmin=leftMin(i);
   rmin=rightMin(i);
   flag=1;
   while flag
       ly=y(lmin);
       ry=y(rmin);
       
       if ly<=ry
           f=(lmin+1-lmin)*ry/(rmin-lmin)+(lmin+1-rmin)*ly/(lmin-rmin);
           if y(lmin+1)-f<0
              lmin=lmin+1;              
           else
               for k=1:rmin-lmin-1
                   ff(k)=(lmin+k-lmin)*ry/(rmin-lmin)+(lmin+k-rmin)*ly/(lmin-rmin);
               end
               if min(y(lmin+1:rmin-1)-ff)<0
                   lmin=lmin+1;
               else
                   flag=0;
               end
               ff=0;
           end
       else
           f=(rmin-1-lmin)*ry/(rmin-lmin)+(rmin-1-rmin)*ly/(lmin-rmin);
           if y(rmin-1)-f<0
              rmin=rmin-1;
           else
               for k=1:rmin-lmin-1
                   ff(k)=(lmin+k-lmin)*ry/(rmin-lmin)+(lmin+k-rmin)*ly/(lmin-rmin);
               end
               if min(y(lmin+1:rmin-1)-ff)<0
                   rmin=rmin-1;
               else
                   flag=0;
               end        
               ff=0;
           end

       end
   end
   leftMin(i)=lmin;
   rightMin(i)=rmin;
end

plotflag=0;
%% display in figure
if plotflag==1    
figure;
plot(Peak,y(Peak),'o',leftMin,y(leftMin),'d',rightMin,y(rightMin),'*');
hold on;plot(y);
for i=1:lenp
    plot([leftMin(i) rightMin(i)],[y(leftMin(i)) y(rightMin(i))],'r','Linewidth',2);
end
end

%% Step 4: Revmove baseline
yy=zeros(1,length(y));
for i=1:lenp
   lmin=leftMin(i);
   rmin=rightMin(i);
   
   ly=y(lmin);
   ry=y(rmin);
   
   for k=lmin+1:rmin-1
       f=(k-lmin)*ry/(rmin-lmin)+(k-rmin)*ly/(lmin-rmin);
       yy(k)=y(k)-f;
   end
end
%figure;plot(y-yy);
yy=yy/scalexy;



