clear variables
tic
% delete(gcp('nocreate'))
addpath(genpath('../../MATLAB functions_KG'));
set_default_figure_parameters(10);

RS_range = [-40:20:3400]';
% IMPORTANT: the spectral range should to be the same as the training data
% Otherwise the material recognition will likely fail

x_bin = 10;

d_pixel = 6.5e-6; % [m] pixel size, PCO
d_pixel_real = d_pixel/25/1.333*1e6; % [um] actual pixel size
D_l = 0.2; % [m] focal length of lens
r_dD = d_pixel/D_l;


shift_factor1 = 609.6e-6*1.333*25/d_pixel; 
shift_factor2 = 681.2e-6*1.333*25/d_pixel*sind(38); 
shift_factor_x = 681.2e-6*1.333*25/d_pixel*cosd(38); 
% [pixel] multiplier to Scanner value for calculating spectrum shift


wl_laser = 660e-3; % [um] Raman laser wavelength, for caculating Raman shift


p_fit = [3.333333333333333   0.0999   0.000003848220834];% 7
%% read beam

[x_im,y_im] = meshgrid(1:2560,1:2160);
y_im = y_im-size(y_im,1)/2;

x_line = 1e3*[0.4508    2.2479]; % 1-4
y_line = [485.8696  456.4533]; % 5-7
beam_centre = (y_line(2)-y_line(1))/(x_line(2)-x_line(1))*(-x_line(1):size(y_im,2)-1-x_line(1))'+y_line(1);
beam_centre = beam_centre-size(y_im,1)/2;
centre_fit = [ones(size(y_im,2),1),(1:size(y_im,2))']\beam_centre;

beam_centre = (y_line(2)-y_line(1))/(x_line(2)-x_line(1))*(-x_line(1):size(y_im,2)-1-x_line(1))'+y_line(1);
beam_centre = beam_centre-size(y_im,1)/2;
centre_fit = [ones(size(y_im,2),1),(1:size(y_im,2))']\beam_centre;

%% read background spectrum
folder = ['7'];
bg_folder = '7/';
bg_tag = 'bg_PCO Edge_';
im_bg = 0;
bg_name = dir([bg_folder bg_tag '*']);
length(bg_name)
for i_bg = 1:length(bg_name)
im_bg = im_bg+double(imread([bg_folder bg_name(i_bg).name]));
end
im_bg = im_bg/length(bg_name);


%% sample spectrum
spectrum_all = [];

scantag = 'y';
spectrum_files = dir([folder '/' scantag '*.tif']);
N_spectrum = length(spectrum_files);
x_range = 1:256;
Sy = NaN([N_spectrum,1]);
S2 = Sy;

x_im_um = x_im*d_pixel_real;
for i_spectrum = 1:N_spectrum
    temp = textscan(spectrum_files(i_spectrum).name,'%s%s%f%f%f%f%f','delimiter','_');
    Sy(i_spectrum) = temp{3};
    S2(i_spectrum) = temp{6};
end

[Sy, i_sort] = sort(Sy);
y_plot = Sy*1e3;

spectrum_files = spectrum_files(i_sort);
x_plot = x_range*d_pixel_real*x_bin;
[x_plot,y_plot] = meshgrid(x_plot,y_plot);


for i_spectrum = 1:5%N_spectrum

    wl_im = xy2wavelength(p_fit,x_im,y_im,centre_fit,r_dD);
    RS_im = (1/wl_laser-1./wl_im)*1e4;

    im_spectrum = double(imread([folder '/' spectrum_files(i_spectrum).name]))-im_bg;
    im_spectrum = imrotate(im_spectrum,-1,'bilinear','crop');
    
    RS_im = RS_im(:,round(x_bin/2):x_bin:end);
    RS_im = RS_im(:,x_range);
    im_spectrum = im_binning(im_spectrum,[1,x_bin]);
    
    figure
    pcolor(repmat(x_range,[size(RS_im,1),1]),RS_im,im_spectrum);   
    shading flat
%     caxis([1e4,4e4])
    
    im_spectrum = im_spectrum(:,x_range);
    
    RS_range_mat = repmat([RS_range;RS_range(end)*2-RS_range(end-1)],1,size(RS_im,1));
    spectrum_out = NaN(length(RS_range),length(x_range));

    parfor i_x = 1:length(x_range)

        RS_x_mat = repmat(RS_im(:,i_x),1,length(RS_range))';
        T_mat = (RS_x_mat>=RS_range_mat(1:end-1,:)).*(RS_x_mat<RS_range_mat(2:end,:));
        spectrum_out(:,i_x) = T_mat*(im_spectrum(:,i_x))./sum(T_mat,2);
        hold on

    end
    i_out = find(~isnan(spectrum_out(round(end/2),:)));   

    spectrum_all = [spectrum_all,spectrum_out(:,i_out)];
    
end

spectrum_all(isnan(spectrum_all)) = 1;
save(['spectrum list ' folder ' ' scantag ' 20bin.mat'],'x_plot','y_plot','spectrum_all','RS_range');
toc

size(spectrum_all)
figure
set(gcf,'units','normalized','outerposition',[0 0 1 1]);
imagesc([0,1],RS_range,spectrum_all)