
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% General description:
%%
%% This function performs peak matching between the observations and
%% library spectra.
%%
%% The whole library paritioning consists of the following steps:
%% 1. Call 'JCSD_background_matching.m' to find the matched background in
%%    'backgrounds' array. 
%%
%% 2. Peak recognition for observations(dataNB) and background 
%%
%% 3. Use peak positions of Nitrogen and Oxygen in observations and library
%%    spectrum to estimate a range of peak off.
%%
%% 4. Matching peaks of observations and library spectra for a peakrange,
%%    and calculate the weight of all matched peaks for each library
%%    spectrum and save in "WeightMatchedLib", and the weight of all matched
%%    peaks in observation and save in "WeightMatchedNB".
%%
%%    "WeightMatchedLib" is to measure how much percentage of library
%%    spectrum is matched.
%%    "WeightMatchedNB" is to measure how much contribution of the library
%%    spectrum is in the observations.      
%%
%% 5. Mearge peaks of observations and background library spectra for a 
%%    range of "dt", extract peaks in observations that are not in bacgrond.
%%    Then eliminate library spectrum in the candidate list whose peaks do
%%    not cover the major extracted peaks.
%%
%% 6. Calculate NNLS between "observations" and "background + library spectrum".
%%    Then return the spectrum with largest correlation
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%% Input:
%%    dataNB - observations that not classified as background
%%    backgrounds - backgrounds data
%%    wavel - wavenumber 
%%    libs - spectrum library
%%    PeakLocsLib - Peak locations of libary spectrum 
%%
%% Output:
%%    output - matched library spectrum index. 
%              =0 if no spectrum matchs observations.
%%
%% Written by: Wei Wang
%% 9/21/2007
%%
%% Copyright(c) 2007 by
%% Machine Learning for Signal Processing Laboratory
%% University of Maryland, Baltimore County
%%
%% Permission to use, copy, modify, and distribute this software for any
%% purpose without fee is hereby granted, provided that this entire notice
%% is included in all copies of any software which is or includes a copy
%% or modification of this software and in all copies of the supporting
%% documentation for such software.
%% THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
%% WARRANTY.  IN PARTICULAR, NEITHER THE AUTHOR DOES NOT MAKE ANY
%% REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
%% OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [output]=JCSD_Detector_PeakMatching(dataNB, backgrounds, wavel, libs,PeakLocsLib)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1. Call 'JCSD_background_matching.m' to find the matched background in
%    'backgrounds' array. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if size(backgrounds,2)==1
    BkgIdChosen=1;
else
    [BkgIdChosen]=JCSD_background_matching(dataNB, backgrounds, wavel);
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 2. Peak recognition for observations(dataNB) and background 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
yBkg0=backgrounds(:,BkgIdChosen);
yNB0=mean(dataNB,2);
[posNB, ampsNB, yNB, y_smoothed]=extracting_peaks(yNB0',wavel);   
[posBkg, ampsBkg, yBkg, y_smoothed]=extracting_peaks(yBkg0',wavel);
locsBkg=wavel(posBkg)';
locsNB=wavel(posNB)';
weightsBkg=ampsBkg/max(ampsBkg);
weightsNB=ampsNB/max(ampsNB);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 3. Use peak positions of Nitrogen and Oxygen in observations and library
%    spectrum to estimate a range of peak off.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
win_width=20;
locNitrogenLib=2327;
peakNitrogen=locNitrogenLib;
locOxygenLib=1556;
peakOxygen=locOxygenLib;
winNitrogen=[locNitrogenLib-win_width locNitrogenLib+win_width];
winOxygen=[locOxygenLib-win_width locOxygenLib+win_width];
num=length(locsBkg);
for i=1:num
   if locsBkg(i) > winNitrogen(1) & locsBkg(i)< winNitrogen(2)
      peakNitrogen=locsBkg(i);        
   end
   
   if locsBkg(i) > winOxygen(1) & locsBkg(i)< winOxygen(2)
      peakOxygen=locsBkg(i);        
   end
end
diff_PeakNitrogen = peakNitrogen - locNitrogenLib;
diff_PeakOxygen = peakOxygen - locOxygenLib;
middiff=(diff_PeakNitrogen+diff_PeakOxygen)/2;
%WinPeak=-[middiff+15 middiff-15];
WinPeak=[-20 10]; % use fixed range for now, might need to change later

%% Debugging block: 
% [m,k]=size(PeakLocs);
% PeakAll=[];
% for i=1:m
%     len=sum(PeakLocs(i,:)>0);
%     tPeak=wavel(PeakLocs(i,1:len));
%     PeakAll=[PeakAll tPeak];
% end
% % PA=unique(PeakAll);
% pk22=PeakLocsLib(22,1:19);
% figure;plot(wavel(pk22),libs(pk22,22),'o-');grid on;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 4. Matching peaks of observations and library spectra for a peakrange,
%    and calculate the weight of all matched peaks for each library
%    spectrum and save in "WeightMatchedLib", and the weight of all matched
%    peaks in observation and save in "WeightMatchedNB".
%
%    "WeightMatchedLib" is to measure how much percentage of library
%    spectrum is matched.
%    "WeightMatchedNB" is to measure how much contribution of the library
%    spectrum is in the observations.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
RangePeaks=[locsNB+WinPeak(1) locsNB+WinPeak(2)];
lenNBLocs=length(locsNB);
mm=size(libs,2);
for i=1:mm
    len=sum(PeakLocsLib(i,:)>0);    
    tPeaks=PeakLocsLib(i,1:len);
    iflagLib=zeros(1,len);
    for j=1:len
        tpos=wavel(tPeaks(j));
        for k=1:lenNBLocs
            if tpos <= RangePeaks(k,1) | tpos >= RangePeaks(k,2)
                continue;
            end
            iflagLib(j)=k;
            break;            
        end
    end
    if sum(iflagLib)==0
        continue;
    end
    
    pks=PeakLocsLib(i,find(iflagLib>0));
    curlibs=libs(:,i);
    amps=curlibs(tPeaks);
    ampsmatched=curlibs(pks);
    ampsmatched=ampsmatched';
    weightPkLib=ampsmatched.^2/norm(amps)^2;
    
    WeightMatchedLib(i)=sum(weightPkLib);

    tposNB=iflagLib(find(iflagLib>0));
    
    ampMatchedNBpk=ampsNB(tposNB);
    for j=1:length(tposNB)
        curPk=tposNB(j);
        if ampsmatched(j)<1e-4
            dddd=1;
        end
        tmprat=ampMatchedNBpk(j)/ampsmatched(j);     
        ampLibinNB=ampsmatched*tmprat;
        tmpDiff=ampMatchedNBpk-ampLibinNB;
        if sum(tmpDiff < 0)== 0 % there is no negative point
            weightPksNB=sum(ampLibinNB.^2)/norm(ampsNB)^2;
            WeightMatchedNB(i)=weightPksNB;
            break;
        end
    end
    
         
%% Debugging block:
%     rtt=max(yNB)/max(curlibs);
%     curlibs=curlibs*rtt;
%     figure;hold on;plot(wavel,yNB);
%     plot(wavel,curlibs,'k');plot(wavel(pks),ampsmatched,'or');grid on;
%                 
%     
end


% The total weight of matched peaks of library spectrum must be greater
% than: threshold_spectrum
threshold_spectrum = 0.85;
iSpec=find(WeightMatchedLib >threshold_spectrum);

% The total weight of matched peaks of in data > threshold_data
threshold_data = 0.02;
iData=find(WeightMatchedNB >threshold_data);

% The candidate spectrum
listCandidate=intersect(iData, iSpec);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 5. Mearge peaks of observations and background library spectra for a 
%    range of "dt", extract peaks in observations that are not in bacgrond.
%    Then eliminate library spectrum in the candidate list whose peaks do
%    not cover the major extracted peaks.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%5.1 Matching peaks for background and dataNB
%
dt=20;
[locsData,idA,idB, ampsNB_all, ampsBkg_all] ...
    = union_peaks(locsNB, locsBkg, ampsNB, ampsBkg, dt);
% 
% % Inverse mapping indices of idA and idB
% lnall=length(locsData);
% invIdxNB=zeros(lnall,1);
% invIdxNB(idA)=[1:length(idA)];
% %invIdxNB=invIdxNB(find(invIdxNB>0));
% invIdxBkg=zeros(lnall,1);
% invIdxBkg(idB)=[1:length(idB)];
% %invIdxBkg=invIdxBkg(find(invIdxBkg>0));

%
% 5.2 rescale ampsBkg to align with Nitrogen peak at '2327'
%
loc_Nitrogen=2327;
tt=abs(locsData-loc_Nitrogen);
posNi=find(tt<40);
if length(posNi)>0
    tampsNi=ampsNB_all(posNi);
    tpps=find(tampsNi==max(tampsNi));
    pNi=posNi(tpps);
    ampsNB_Nitrogen=ampsNB_all(pNi);
    ampsBkg_Nitrogen=ampsBkg_all(pNi);
    rat=ampsNB_Nitrogen/ampsBkg_Nitrogen;

    ampsBkg_all=ampsBkg_all*rat;
end

%
% 5.3 Extract peaks in observations that are not in bacgrond.
%
diffNB_Bkg=ampsNB_all-ampsBkg_all;
BigPeak_threshold=0.4;
pksRatio=ampsNB/max(ampsNB);
BigPeakLocs = find(pksRatio>BigPeak_threshold);
pksRatioAll=diffNB_Bkg/max(diffNB_Bkg);
BigPeakLocsAll = find(pksRatioAll>BigPeak_threshold);
BigPeakNBLocs=setdiff(BigPeakLocsAll,idB); % Big peaks extracted
                                           % from observations

%
% 5.4 Eliminate library spectrum in the candidate list whose peaks do
%    not cover big extracted peaks.
numBigPeaks=length(BigPeakNBLocs);
numCandidate=length(listCandidate);
flagCandidate=ones(numCandidate,1);
for i=1:numBigPeaks
    curBigPeakLoc=locsData(BigPeakNBLocs(i));
    for k=1:numCandidate
        if flagCandidate(k)==0
            continue;
        end
        ilib=listCandidate(k);       
        len=sum(PeakLocsLib(ilib,:)>0);    
        tPeaks=PeakLocsLib(ilib,1:len);
        tpos=wavel(tPeaks);
        
        ddt=abs(tpos-curBigPeakLoc);
        if ddt > win_width
            flagCandidate(k)=0;
        end
    end
    
end

if sum(flagCandidate)==0
    output=0;
    return;
elseif sum(flagCandidate)==1
    output=listCandidate(find(flagCandidate==1));
    return;
end
        
listCandidate=listCandidate(find(flagCandidate==1));



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 6. Calculate NNLS between "observations" and "background + library spectrum".
%    Then return the spectrum with largest correlation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
len=length(listCandidate);
for i=1:len
    ilib=listCandidate(i);

    len=sum(PeakLocsLib(ilib,:)>0);    
    tPeaks=PeakLocsLib(ilib,1:len);
    locsLib=wavel(tPeaks);
    curlib=libs(:,ilib);
    curlib=curlib*max(yNB)/max(curlib)/2;
    ampsLib=curlib(tPeaks);
% 
%     sit=strcat('libnumber: ',int2str(ilib));
%     figure;
%     hold on;
%     plot(wavel,yNB);
%     plot(wavel,yBkg,'k');
%     plot(wavel,curlib,'r');
%     hold off;
%     title(sit);
    
    
    [locsAll,idxData,idxLib,ppNB,ppLib]...
        = union_peaks(locsData, locsLib, ampsNB_all, ampsLib, dt);
    topNB=max(ppNB);
    
  %  figure;plot(locsAll,ppNB,'o-'); grid on;title('dataNB');
    m=length(locsAll);
    ppBkg=zeros(m,1);
    ppBkg(idxData)=ampsBkg_all;
    ppBkg=ppBkg*(topNB/max(ampsBkg_all));
  %  figure;plot(locsAll,ppBkg,'o-'); grid on;title('background');
    ppLib=ppLib*(topNB/max(ampsLib));
 %   figure;plot(locsAll,ppLib,'o-'); grid on;title(sit);
    
    %% Correlation between 'dataNB-background' and 'libspectrum'
    ddData=ppNB-ppBkg;
    cors(i)=corr2(ddData,ppLib);    
 %   figure;plot(locsAll,ddData,'o-'); grid on;title('ddData');
 
    tlibs=[ppBkg ppLib];
    [p,idx]=NonNegativeLS(tlibs,ppNB);
    if norm(p)==0 
        cc(i)=0;
    else
        cc(i)=corr2(ppNB, tlibs*p);    
    end
    
    [a,b,r]=canoncorr(tlibs,ppNB);
    rr(i)=r(1);
end
maxc=max(cc);

maxcorr_threshold=0.85;
if maxc<maxcorr_threshold
    output=0;
    return;
end
pos=find(cc==maxc);
if length(pos)==0 
    output=0;
    return;
end
output=listCandidate(pos(1));
% 
% figure;plot(WeightMatchedLib); grid on;
% 
% figure;plot(WeightMatchedNB); grid on;
