
% JOINT CONTAMINATED SURFACE DETECTOR
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% General description:
%%
%% This program performs detection on given observation data file
%%
%% This program consists of 3 steps:
%% 1. Prepare spectrum library for detection
%%  1.1 Determine spectrum range and baseline correction option
%%  1.2 Load library information: (libs, x, clusters, list_corr_pairs)
%%  1.3 Adjust library spectrum range to fit observation data if necessary
%% 2. Load observation data, and extract non-background data
%%  2.1 Read date file
%%  2.2 Remove those abnormal data
%%  2.3 Estimate background spectrum
%%  2.4 Extract non-background data, saved in dataNB, and indices in
%%      globalIdxNB
%% 3. Display results
%%
%% Parameters:
%%    datafilename - datafile name 
%%    PresentChemicalIdx - index of present chemical, used for displaying 
%%                         detection results
%%    flag_clean_datafile - when data file is clean, used for displaying
%%                          detction result
%%    maxpulsenumber - the number of observations used for detection
%%    blocksize - size of data block used for detection
%%    readblocksize - size of observations used for reading data file;
%%    thresholdBkg - correlation threshold to be regarded as the same spectrum
%%    confidence - confidence level in determining non-background threshold
%%
%% Written by: Wei Wang
%% 8/28/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.
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all;
close all;
% % % % 
% specify data file name
datafilename='./asphalt_tmes1.dat'; %37890, detected!
PresentChemicalIdx=3;  % index of present chemical, used for displaying 
                        % detection results
flag_clean_datafile=0;  % =1 when data file is clean, used for displaying
                        % detction result

% flag_clean_datafile=1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 1. Prepare spectrum library information
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 1.1 Determine spectrum range and baseline correction option
%
flag_baseline_corrected=1;
flag_narrow_range=0;
if flag_narrow_range==0
    rangeleft=500;
    rangeright=3200;
else    
    rangeleft=500;
    rangeright=1700;
end
    
%
% 1.2 Load library information: (libs, x, clusters, list_corr_pairs)
%     The library file is obtained by "JCSD_LibraryPreprocessing_main.m"
%
preprocessflag=0;
if flag_baseline_corrected==1 & flag_narrow_range==1
    preprocessflag=1;
    strflag='baseline corrected, narrow spectrum';
    load JCSDLibs_BaselineCorrected_NarrowSpectrum.mat;
elseif flag_baseline_corrected==1 & flag_narrow_range==0
    preprocessflag=2;
    strflag='baseline corrected, broad spectrum';
    load JCSDLibs_BaselineCorrected_BroadSpectrum.mat;    
elseif flag_baseline_corrected==0 & flag_narrow_range==0
    preprocessflag=3;
    strflag='with baseline, broad spectrum';
    load JCSDLibs_BroadSpectrum.mat;
else    
    preprocessflag=4;    
    strflag='with baseline, narrow spectrum';
    load JCSDLibs_NarrowSpectrum.mat;
end

%
% 1.3 Adjust library spectrum range to fit observation data if necessary
%

% get spectrum range of observation data: 'wavenumbers'
[data, wavenumbers, header] = JCSD_Matlab_Read3(datafilename, [100:101]);
 
irange=find(wavenumbers>rangeleft & wavenumbers< rangeright);
wavel=wavenumbers(irange);

m=size(libs,2);
% If spectrum range not matching, expand first, then interpolate
if  length(wavel)~=length(x) | norm(wavel-x)>1e-4 
    
   % expand x to cover wavel
   if x(1) < wavel(1) % 'x' is decreasing 
       x=[wavel(1);x];
       libs=[libs(1,:);libs];
       PeakLocs=PeakLocs+1; % update peak locations
%       libs=[zeros(1,m);libs];
   end       
   if x(end) > wavel(end) % 'x' is decreasing 
       x=[x;wavel(end)];
       libs=[libs;libs(end,:)];
   end       

   % interpolate library to observation data spectrum range
   for i = 1:size(libs,2)
       libsall(:,i)=interp1(x,libs(:,i),wavel);
       
       % update peaks locations
       len=sum(PeakLocs(i,:)>1);
       tPeaks=x(PeakLocs(i,1:len));
       for j=1:len
           tx=tPeaks(j);
           itx=abs(wavel-tx);
           imin=find(itx==min(itx));
           newPeakLocs(i,j)=imin;
       end
   end
   PeakLocs=newPeakLocs;
   
end 
libs=libsall;
JCSDAgentMaxID=30;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 2 Load observation data, and perfrom detection
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 2.1 read observation data from data file, and remove abnormal data
%% 2.2 Extract the backgrounds and non-background data
%% 2.3 Determine threshold based  on background data
%% 2.4 Perform DCCP detection and peak matching on non-background data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

maxpulsenumber=10000;%37890;
blocksize=10; % detection blocksize
readblocksize=500;
nn=floor(maxpulsenumber/readblocksize);
thresholdBkg=0.9; % correlation threshold to be regarded as the same spectrum
confidence=0.95; % confidence level in determining non-background threshold

sst=sprintf('blksize=%d, readblksize=%d, threshold=%4.2f, confidence=%4.2f',...
    blocksize,readblocksize, thresholdBkg, confidence);

timeDetection=0;
timeRead=0;
timeBkgEst=0;
timeBkgDCCP=0;
count=0;
countAll=0;
for kk=1:nn
    kk
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % 1. read observation data from data file, and remove abnormal data
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    istart=(kk-1)*readblocksize+1;
    iend=kk*readblocksize; 
    tic;
    [data, wavenumbers, header] = JCSD_Matlab_Read3(datafilename, [istart:iend]);
    timeRead=timeRead+toc;
    
    data=data';
    data=data(irange,:);
    
    [data,idx]=JCSD_DataCleaning(data); % Remove those abnormal data

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % 2. Extract the backgrounds and non-background data
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    tic;
    [backgrounds, curIdxNB] = JCSD_BackgroundExtractor(data,wavel);
    timeBkgEst=timeBkgEst+toc;
    
    idxBlock=[istart:iend];
    idxGlobalBlock=idxBlock(idx);
    
    % save non-background and background data
    dataNB=data(:,curIdxNB);
    dataBkg=data;
    dataBkg(:,curIdxNB)=[];
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % 3 Determine threshold based  on background data
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    tic;
    [threshold_DCCP_Bkg(kk),avg_DCCP_Bkg(kk)]=JCSD_threshold_background_DCCP...
        (dataBkg, libs, clusters, blocksize, confidence);
    timeBkgDCCP=timeBkgDCCP+toc;
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % 4 Perform DCCP detection and peak matching on non-background data
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    tic;
    [n,m]=size(dataNB);
    mm=ceil(m/blocksize);
    for i=1:mm
        curStart=(i-1)*blocksize+1;
        curEnd=min(m,i*blocksize);
        x=dataNB(:,curStart:curEnd);    
    
        %
        % 4.1 DCCP detection
        %
        [score_DCCP, targetId_DCCP]=JCSD_detector_DCCP(x,libs,wavel,clusters, list_corr_pairs);

        %
        % 4.2 Peak matching detection
        %         
        if score_DCCP>threshold_DCCP_Bkg(kk)
           targetIdx=JCSD_Detector_PeakMatching(x, backgrounds, wavel, libsall(:,1:JCSDAgentMaxID), PeakLocs);
           
           % Only keep non-zero identification results
           if targetIdx > 0
               count=count+1;
               targetDetected(count)=targetIdx;
               dccpDetected(count)=score_DCCP;
               globalPosDetected(count)=idxGlobalBlock(curStart);
           end
           
           % targetDetectedAll save both zero and non-zero identification
           % results. This is used for displaying resutls from clean
           % data files, in which all zero are correct detection results.
           countAll=countAll+1;
           targetDetectedAll(countAll)=targetIdx;
           dccpDetectedAll(countAll)=score_DCCP;
           globalPosDetectedAll(countAll)=idxGlobalBlock(curStart);
        end
        
    end
    timeDetection=timeDetection+toc;
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 3 Display results
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% prepare for threthold and average from background
for kk=1:nn
    segmentBkg([2*(kk-1)+1 2*kk])=[(kk-1)*readblocksize+1 kk*readblocksize];
    segmentDCCPthreshold([2*(kk-1)+1 2*kk])=threshold_DCCP_Bkg(kk);
    segmentDCCPavg([2*(kk-1)+1 2*kk])=avg_DCCP_Bkg(kk);
end


if flag_clean_datafile==0
    figure;
    subplot(2,1,1);
    hold on;
    plot(segmentBkg,segmentDCCPthreshold,'r-','LineWidth',2);
    plot(segmentBkg,segmentDCCPavg,'k-','LineWidth',2);
    plot(globalPosDetected,dccpDetected,'o','LineWidth',3);
    grid on;
    xlabel('Pulse index');
    ylabel('DCCP Score');
    ylim([0.5 1]);
    box on;
    title(sst);
    hold off;

    subplot(2,1,2);
    hold on;
    tt=find(targetDetected==PresentChemicalIdx);
    plot(globalPosDetected(tt),targetDetected(tt),'o','LineWidth',3);
    ttp=[1:length(globalPosDetected)];
    ttp(tt)=[];
    plot(globalPosDetected(ttp),targetDetected(ttp),'vr','LineWidth',3);
    DetectionRate=length(tt)/length(globalPosDetected);
    s11=sprintf('Correct detection (%d%s',round(DetectionRate*100),'%)');
    legend(s11,'Wrong detection');

    HitRate=length(tt)*blocksize/maxpulsenumber;
    s22=sprintf('Hit: %d%s',round(HitRate*100),'%');
    xpos=maxpulsenumber/10;
    ypos=45;
    text(xpos,ypos,s22,'color','b','fontsize',16);
    xlabel('Pulse index');
    ylabel('Detected spectrum');
    ylim([0 50]);
    xlim([0 maxpulsenumber]);
    hold off;
    grid on;    
    box on;
end
 
% 
% Displaying all non-zero identification resutls. This is for clean
% data files in which all zero are correct detection results.
%
if flag_clean_datafile==1
    figure;
    subplot(2,1,1);
    hold on;
    plot(segmentBkg,segmentDCCPthreshold,'r-','LineWidth',2);
    plot(segmentBkg,segmentDCCPavg,'k-','LineWidth',2);
    plot(globalPosDetectedAll,dccpDetectedAll,'o','LineWidth',3);
    grid on;
    xlabel('Pulse index','fontsize',20);
    ylabel('DCCP Score','fontsize',20);
    ylim([0.5 1]);
    hold off;

    subplot(2,1,2);
    plot(globalPosDetectedAll,targetDetectedAll,'o','LineWidth',3);
    xlabel('Pulse index','fontsize',20);
    ylabel('Detected spectrum','fontsize',20);
    ylim([0 50]);
    xlim([0 maxpulsenumber]);
    hold off;
    grid on;
end



