%Copyright (c) 2010, Luuk J.G.W. van Wilderen
%
%dispersion_IRF_calc.m

function dispersion_IRF_calc
%Calculates cross-correlation between each wavelength and impulse reponse
%function (IRF). The maximum of the derivative of the cross-correlation (of
%each signal) determines the delay between signal and IRF.
%Choose between gaussian and step function (which is a gaussian with a
%width that is smaller than the measured time step, i.e. your time
%resolution). Importantly, the data needs to be resampled to have regularly
%spaced time intervals for the cross-correlation to be properly calculated
%(requirement of 'xcorr2').
%User can select order n of the dispersion function
%mu=mu_0+sum(i->n, mu_i*( (lambda-lambda_c)/100 )^i ), and needs to supply
%starting values in the GUI.
%Remember that the IRF that is fitted here is the optimum value for the
%dispersion determination, and not necessarily the IRF of the data.

%Defined here
global filename
%Defined in plot_raw_data
global ydata_raw xdata_raw
%Defined in load_calibration
global calibfilename 
%Defined in fit_data
global wave
%Defined here
global delay_shift_curve

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Running some checks and loading of data and input parameters
if isempty(filename)==1
    errordlg('Load dataset first!')
    error('Load dataset first!')
end

%Defined in load_dispersion.
global dispersfilename
%Clear dispersion data from memory if 'Fit dispersion' button is pressed.
if exist('dispersfilename')==1 && isempty(dispersfilename)==0
    clear dispersfilename
    set(findobj('Tag','dispersion_filename_text'),'Visible','on','String','Dispersion file:  ');
end
clear xdata_raw ydata_raw wave
raw_data=load (filename);
xdata_raw=raw_data(:,1);

%Check if time points monotonically increase
if isequal(sort(xdata_raw),xdata_raw)==0
    errordlg('Time points not monotonically increasing. Check first column in data.')
    error('Time points not monotonically increasing. Check first column in data.')
end
ydata_raw=raw_data(:,2:end);

Dispers_pars=str2num(get(findobj('Tag','Dispers_box'),'String')); %order (n), time zero (muo), central wavelength (lc), and m(n) for each order n
disp_order=Dispers_pars(1);
if length(Dispers_pars)>1
    if disp_order+3 ~= length(Dispers_pars)
        errordlg(['Number of dispersion parameters does not match dispersion order, should be ',num2str(disp_order+3)])
        error(['Number of dispersion parameters does not match dispersion order, should be ',num2str(disp_order+3)]);
    end
    start_pars_dispers=Dispers_pars(1:disp_order+3);
elseif length(Dispers_pars)==1
    errordlg('Define all dispersion parameters first!')
    error('Define all dispersion parameters first!');
end
Par_fix_disp=get(findobj('Tag','Par_fix_disp_toggle'),'Value');

%Input parameters
%start_pars_dispers: time zero (mu0), order (n), central wavelength (lc), and m(n) for each n orders
mu_pars=[];
mu_pars(1)=start_pars_dispers(2);       %time zero (mu0)
mu_pars(3)=start_pars_dispers(3);       %central wavelength (lc)
mu_pars(2)=start_pars_dispers(4);       %m(1) for first order
if disp_order>1
    for i=5:length(start_pars_dispers)
        mu_pars(i-1)=start_pars_dispers(i); %m(n) for each n orders
    end
end

%Definition of data
xdata=xdata_raw;
ydata=ydata_raw;

if exist('calibfilename')==1 && isempty(calibfilename)==0
    wave=load (calibfilename);
    if length(wave)~=size(ydata,2)
        errordlg('Calibration file does not match datafile!');
        error('Calibration file does not match datafile!');
    end
    if size(wave,2)~=length(wave)
        wave=wave';
    end
elseif exist('calibfilename')==0 || isempty(calibfilename)==1
    wave=1:size(ydata,2);
end

%Wavelength and time range selection
if str2double(get(findobj('Tag','wave_lim_box1'),'String'))==0 && str2double(get(findobj('Tag','wave_lim_box2'),'String'))==0
    wave_limits=[1 length(wave)];
elseif str2double(get(findobj('Tag','wave_lim_box1'),'String'))~=0 && str2double(get(findobj('Tag','wave_lim_box2'),'String'))~=0
    wave_limits(1)=str2double(get(findobj('Tag','wave_lim_box1'),'String'));
    wave_limits(2)=str2double(get(findobj('Tag','wave_lim_box2'),'String'));
end
if wave_limits(1)>wave_limits(2)
    errordlg('Left boundary of dispersion wavelength range must be smaller than right boundary.');
    error('Left boundary of dispersion wavelength range must be smaller than right boundary.')
end
mu_wave_sel=[];
mu_wave_sel=wave(wave_limits(1):wave_limits(2));

if str2double(get(findobj('Tag','time_lim_box1'),'String'))==0 && str2double(get(findobj('Tag','time_lim_box2'),'String'))==0
    time_limits=[1 length(xdata_raw)];
elseif str2double(get(findobj('Tag','time_lim_box1'),'String'))~=0 && str2double(get(findobj('Tag','time_lim_box2'),'String'))~=0
    time_limits(1)=str2double(get(findobj('Tag','time_lim_box1'),'String'));
    time_limits(2)=str2double(get(findobj('Tag','time_lim_box2'),'String'));
end

if (time_limits(1)==0 && time_limits(2)~=0) || (time_limits(1)~=0 && time_limits(2)==0)
    errordlg('Both time boundaries in dispersion range must be unequal to zero if range should be applied.');
    error('Both time boundaries in dispersion range must be unequal to zero if range should be applied.')
elseif time_limits(1)==0 && time_limits(2)==0
    time_limits(1)=1;
    time_limits(2)=length(xdata);
end
if time_limits(1)>=time_limits(2)
    errordlg('Left boundary of dispersion time range must be smaller than right boundary.');
    error('Left boundary of dispersion time range must be smaller than right boundary.')
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Data manipulation just like in globfit_bio_gui
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Time/wavelength shifter
time_shift=str2double(get(findobj('Tag','time_shift_box'),'String'));
wave_shift=str2double(get(findobj('Tag','wave_shift_box'),'String'));
xdata=xdata+time_shift; %Shift time axis with absolute amount
wave=wave+wave_shift;       %Shift wavelength axis with absolute amount

%Baseline correction(s): Or whole spectrum at selected time point (or average
%of range of time points), or time trace at selected wavelength (or average
%of range of wavelengths
sub_timepoint=str2num(get(findobj('Tag','time_subtract_box'),'String'));
sub_wavelength=str2num(get(findobj('Tag','wave_subtract_box'),'String'));
nr_waves=length(wave);          %Number of wavelengths
nr_times=length(xdata);         %Number of time points
if length(sub_timepoint)==1 && sub_timepoint>0
    ydata_sub_timepoint=repmat(0,nr_times,nr_waves);    %Parameter initilisation
    for i=1:size(ydata,1)
        ydata_sub_timepoint(i,:)=ydata(i,:)-ydata(sub_timepoint,:);
    end
    clear i
    ydata=ydata_sub_timepoint;
elseif length(sub_timepoint)==2 && sub_timepoint(1)>0 && sub_timepoint(2)>0
    ydata_sub_timepoint=repmat(0,nr_times,nr_waves);    %Parameter initilisation
    ydata_sub_timepoint_ave=sum(ydata(sub_timepoint(1):sub_timepoint(1),:),1)/(abs(sub_timepoint(2)-sub_timepoint(1))+1);
    for i=1:size(ydata,1)
        ydata_sub_timepoint(i,:)=ydata(i,:)-ydata_sub_timepoint_ave;
    end
    clear i
    ydata=ydata_sub_timepoint;
end
clear ydata_sub_timepoint ydata_sub_timepoint_ave
if length(sub_wavelength)==1 && sub_wavelength>0
    ydata_sub_wavelength=repmat(0,nr_times,nr_waves);    %Parameter initilisation
    for j=1:size(ydata,2)
        ydata_sub_wavelength(:,j)=ydata(:,j)-ydata(:,sub_wavelength);
    end
    clear j
    ydata=ydata_sub_wavelength;
elseif length(sub_wavelength)==2 && sub_wavelength(1)>0 && sub_wavelength(2)>0
    ydata_sub_wavelength=repmat(0,nr_times,nr_waves);    %Parameter initilisation
    ydata_sub_wavelength_ave=sum(ydata(:,sub_wavelength(1):sub_wavelength(2)),2)/(abs(sub_wavelength(2)-sub_wavelength(1))+1);
    for j=1:size(ydata,2)
        ydata_sub_wavelength(:,j)=ydata(:,j)-ydata_sub_wavelength_ave;
    end
    clear j
    ydata=ydata_sub_wavelength;
end
clear ydata_sub_wavelength ydata_sub_wavelength_ave
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%End data manipulation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Definition of work dataset
data_sel=timeseries(ydata(time_limits(1):time_limits(2),wave_limits(1):wave_limits(2)),xdata(time_limits(1):time_limits(2)));

%Interpolate testdata such that it contains regularly spaced time points.
%Otherwise the delays (with which the data has to be shifted) that are 
%calculated by the cross-correlation between IRF and testdata have no
%meaning.
tic
time_interp=(data_sel.time(1):(data_sel.time(end)-data_sel.time(1))/length(data_sel.time):data_sel.time(end))';
testdata=interp1q(data_sel.time,data_sel.data,time_interp);
toc
disp('Done resampling')

%Define step function as instrument response function (normalised to
%data, but this is not important for the funcionality of the script).
time_interp_gauss_spacing=time_interp(2)-time_interp(1);
% time_interp_gauss=-10*time_interp_gauss_spacing:time_interp_gauss_spacing:10*time_interp_gauss_spacing;
%Generate normalised IRF function
IRF_norm=gaussmf(time_interp,[0.1*time_interp_gauss_spacing 0])*max(max(abs(testdata)));

tic
%The maximum of the cross correlation between IRF and each signal
%determines their respective delay. Time zero for the IRF is at zero (the
%reference point). Note that IRF and testdata MUST have the same time base.
for i=1:size(testdata,2)
%     [tmp,delay_shift_i]=max(abs(xcorr(testdata(:,i),IRF_norm)));
    [tmp,delay_shift_i]=max(diff(abs(xcorr(testdata(:,i),IRF_norm))));
    delay_shift(i)=delay_shift_i;
end
%By definition the row matrix delay_shift contains more time points than 
%signal actually contains. Correct for this, and make them relative to the 
%reference point (so create pos/neg values).
delay_shift=delay_shift-max([size(testdata,1) length(IRF_norm)]);
clear tmp max_xcor
toc
disp('Autocorrelation calculated.')

% % % Diagnostics; Switch on to check estimated delay at single pixel.
% px=250;
% figure;subplot(3,1,2);
% plot(delay_shift,'b')
% subplot(3,1,1);
% plot(testdata(:,px),'r');
% hold on;
% plot(IRF_norm./max(IRF_norm)*max(testdata(:,px)),'k')
% legend('testdata','IRF')
% title(['pixel ',num2str(px)])
% subplot(3,1,3);
% plot(diff(abs(xcorr(testdata(:,px),IRF_norm))),'k')
% subplot(3,1,2);
% hold on;
% plot(delay_shift*time_interp_gauss_spacing,'g')
% axis tight

%Convert estimated delay shift from 'number of time points' to acutal times.
delay_shift_time=delay_shift'*time_interp_gauss_spacing;

%Fitting routine
if Par_fix_disp==0
    %If the dispersion parameters need to be fitted
    disp('Fitting of dispersion, please wait....')
    
    %Definition of fit function.
    if disp_order==1
        f = fittype('a+b*((x-c)/100)^1');
    elseif disp_order==2
        f = fittype('a+b*((x-c)/100)^1+d*((x-c)/100)^2');
    elseif disp_order==3
        f = fittype('a+b*((x-c)/100)^1+d*((x-c)/100)^2+e*((x-c)/100)^3');
    elseif disp_order==4
        f = fittype('a+b*((x-c)/100)^1+d*((x-c)/100)^2+e*((x-c)/100)^3+f*((x-c)/100)^4');
    elseif disp_order==5
        f = fittype('a+b*((x-c)/100)^1+d*((x-c)/100)^2+e*((x-c)/100)^3+f*((x-c)/100)^4+g*((x-c)/100)^5');
    end
    
    %Calling of fit routine. Has to be a robust method: traces that exhibit
    %no signal, or signals that change sign over the wavelength range,
    %potentially give erratic results when the cross-correlation with an
    %IRF is calculated.
    [fit1] = fit(mu_wave_sel',delay_shift_time,f,'StartPoint',mu_pars,'robust','on','Normalize','off');

    %       Not needed for script.
    %       Only for saving extra data:    
    %       [fit2] = fit(mu_wave_sel',delay_shift',f,'StartPoint',mu_pars,'robust','on','Normalize','off');
    clear mu_wave_sel

elseif Par_fix_disp==1
    %If all dispersion parameters are fixed
    disp('Fixed dispersion parameters:')
    clear fit1;
    a=mu_pars(1);
    b=mu_pars(2);
    c=mu_pars(3);
    fit1.a=a;
    fit1.b=b;
    fit1.c=c;
    if disp_order==1
        fun1=@(x) a+b*((x-c)/100).^1;
    elseif disp_order>1
        d=mu_pars(4);
        fit1.d=d;
        fun1=@(x) a+b*((x-c)/100).^1+d*((x-c)/100).^2;
    elseif disp_order>2
        e=mu_pars(5);
        fit1.e=e;
        fun1=@(x) a+b*((x-c)/100).^1+d*((x-c)/100).^2+e*((x-c)/100).^3;
    elseif disp_order>3
        f=mu_pars(6);
        fit1.f=f;
        fun1=@(x) a+b*((x-c)/100).^1+d*((x-c)/100).^2+e*((x-c)/100).^3+f*((x-c)/100).^4;
    elseif disp_order>4
        g=mu_pars(7);
        fit1.g=g;
        fun1=@(x) a+b*((x-c)/100).^1+d*((x-c)/100).^2+e*((x-c)/100).^3+f*((x-c)/100).^4+g*((x-c)/100).^5;
    end
end

%Displaying of results
disp(['Dispersion order: ', num2str(disp_order)])
disp(['    mu_0: ', num2str(fit1.a)])
disp([' lamda_c: ', num2str(fit1.c)])
disp(['    mu_1: ', num2str(fit1.b)])
if disp_order>1
    disp(['    mu_2: ', num2str(fit1.d)])
elseif disp_order>2
    disp(['    mu_3: ', num2str(fit1.e)])
elseif disp_order>3
    disp(['    mu_4: ', num2str(fit1.f)])
elseif disp_order>4
    disp(['    mu_5: ', num2str(fit1.g)])
end

%Plotting
if isempty(findobj('Name','Dispersion estimation'))==0 %If figure exists already, delete it
    close ('Dispersion estimation');
end
dispersion_data_plot=figure;
set(dispersion_data_plot,'name','Dispersion estimation','Units',...
    'normalized','position',[0.02 0.05 0.97 0.8],'NumberTitle','off',...
    'Toolbar','figure','WindowStyle','docked');

if Par_fix_disp==1
    fit1=fun1;
end

subplot(2,1,1)
surf(wave,xdata(time_limits(1):time_limits(2)),ydata(time_limits(1):time_limits(2),:));     %Raw data
shading interp
view([0 90])
axis tight
hold on
plot(wave,fit1(wave),'k');
line([wave(wave_limits(1)) wave(wave_limits(2))],[data_sel.time(1) data_sel.time(1)],'color','k','linewidth',2)
line([wave(wave_limits(1)) wave(wave_limits(2))],[data_sel.time(end) data_sel.time(end)],'color','k','linewidth',2)
line([wave(wave_limits(1)) wave(wave_limits(1))],[data_sel.time(1) data_sel.time(end)],'color','k','linewidth',2)
line([wave(wave_limits(2)) wave(wave_limits(2))],[data_sel.time(1) data_sel.time(end)],'color','k','linewidth',2)

title([{'Raw data, and selected wavelength/time region (black box) for cross-correlation calculation'},{'Thin black line is the estimated dispersion curve from the plot below'}])
xlabel('Wavelength')
ylabel('Time')
set(gca,'YLim',[data_sel.time(1) data_sel.time(end)])
set(dispersion_data_plot,'renderer','opengl')

subplot(2,1,2)
plot(wave(wave_limits(1):wave_limits(2)),delay_shift_time)
hold on
plot(wave,fit1(wave),'k');
axis tight
line([wave(wave_limits(1)) wave(wave_limits(1))],[get(gca,'yLim')],'color','k','linewidth',2)
line([wave(wave_limits(2)) wave(wave_limits(2))],[get(gca,'yLim')],'color','k','linewidth',2)
title('Delay of time zero calculated by maximum of cross-correlation of IRF and measured signal')
xlabel('Wavelength')
ylabel('Time')

%Generate dispersion curve for all time points and a selection of 
%wavelengths, rounded to nearest time point index
% delay_shift_y_fit_raw_curve=round(fit2(wave(wave_limits(1):wave_limits(2))))';
% disp_save_data=[wave(wave_limits(1):wave_limits(2));delay_shift;fit2(wave(wave_limits(1):wave_limits(2)))';round(fit2(wave(wave_limits(1):wave_limits(2))))';delay_shift_time';fit1(wave(wave_limits(1):wave_limits(2)))'];
% clear disp_save_data

set(findobj('Tag','dispersion_filename_text'),'Visible','on','String','Dispersion file: Dispersion curve available from last dispersion fit. Can be applied, but is not saved (and loaded) yet.');

delay_shift_curve=fit1(wave);