clear all;

% Function written by Alvaro H Crevenna
% LMU Munich, 2013

cycle_time = 1; % Acquisition rate in seconds e.g. 1 frame every 2s is 2 cycle time
PIX = 98; % Pixel size in nm

filename = 'Kymo_test '; % FIlename to save data
Iraw = imread('test.tif'); % Input the name of the kymo to analyse
Ird = double(Iraw);
% Options
check = 1; %displays an overlay of raw image and found ends
Barbed_end_first =0; % 0 or 1 % Whether the b-end appears first from left to right
Mean = 1; % 0 or 1; If the mean elongation rate should be printed
Depol = 0; % Is the change negative? i.e. the filament end is depolymerizing
Double_end = 0; % 0 or 1; are both ends needed to analyze or only one?
Emin = 1; % Threshold to calculate pausing probability
NP = 11; % Number of points to average to estimate running mean

if Barbed_end_first == 0
    Ird = flipdim(Ird,2);
end


% Barbed end

% This next section is to manually draw a line close to the filament end
% to be analysed

figure; imshow(Ird,[]);
p = impoly(gca,[]);
position = wait(p);
close


Time_axis_b = [min(position(:,2)):max(position(:,2))]';
Length_axis_b = interp1(position(:,2),position(:,1),Time_axis_b);

s = floor(Time_axis_b(end)-Time_axis_b(1)+1);

% Here the end will be fitted with an error function
w = 5;
for i=1:s
    lb = floor(Length_axis_b(i)-w); up = floor(Length_axis_b(i)+w);
    m = [lb:up];
    Line = Ird(floor(Time_axis_b(i)),lb:up);
    
    pars2=iCDFfit(m',Line', [mean(Line(1:3)) Length_axis_b(i) -1 120],[20 Length_axis_b(i)-10 -5 1],[2000 Length_axis_b(i)+10 10 500]);
    XX(i) = pars2(2);
    rowc(i) = XX(i);
    StDevGa(i) = abs(pars2(3));    YYFit = surGaussFunct(pars2,m);
    %     plot(m,Line,'o',m,YYFit,'r')
    
end

% To check visualize raw + detection
if size(Time_axis_b,1) > max(size(rowc))
    Time_axis_b = Time_axis_b(1:end-1);
end

if check==1
    figure; imshow(Ird,[]); hold on; plot(rowc, Time_axis_b(1:end)','ro'); hold off
end

Fil_Lum = abs(rowc-rowc(1))*(0.001*PIX);
if Depol ==1
    Fil_Lum = -Fil_Lum;
end
Fil_Lsub_b = Fil_Lum*370;
Time_b = Time_axis_b(1:end).*(cycle_time);

% Moving average X points
NP = 11; % Number of points to average
LengthRM_b = runmean(Fil_Lsub_b, NP);

figure;
subplot(2,1,1)
plot(Time_b,Fil_Lsub_b,'.',Time_b,LengthRM_b,'r','LineWidth',2)
xlabel('time (s)')
ylabel('\DeltaLength (sub)')
title(filename)
dF_b = diff(LengthRM_b)./diff(Time_b)';
[n_b,x_b] = hist(dF_b,[floor(min(dF_b)):1:ceil(max(dF_b))]);
P_b = 1- sum((dF_b>Emin)/length(dF_b));

subplot(2,1,2)
plot(Time_b(2:end),dF_b,'')
xlabel('time (s)')
ylabel('Elongation Rate (sub/s)')
text(Time_b(1),0.75*max(dF_b),['Pause probability =', num2str(P_b)] )

if Mean ==1
    mdF_b = mean(dF_b(dF_b<0))
end


% Pointed-end
if Double_end == 1;
    Ird_p = flipdim(Ird,2);
    Depol = 1;
    figure; imshow(Ird_p,[]);
    p_p = impoly(gca,[]);
    position_p = wait(p_p);
    Time_axis_p = [min(position_p(:,2)):max(position_p(:,2))]';
    Length_axis_p = interp1(position_p(:,2),position_p(:,1),Time_axis_p);
    
    s_p = floor(Time_axis_p(end)-Time_axis_p(1)+1);
    
    for i=1:s_p
        lb_p = floor(Length_axis_p(i)-w); up = floor(Length_axis_p(i)+w);
        m_p = [lb_p:up];
        Line_p = Ird_p(floor(Time_axis_p(i)),lb_p:up);
        
        pars2_p=iCDFfit(m_p',Line_p', [mean(Line_p(1:3)) Length_axis_p(i) -1 120],[20 Length_axis_p(i)-5 -5 1],[2000 Length_axis_p(i)+5 10 500]);
        XX_p(i) = pars2_p(2);
        rowc_p(i) = XX_p(i);
        StDevGa_p(i) = abs(pars2_p(3));    YYFit = surGaussFunct(pars2_p,m_p);
        
    end
    
    
    
    % To check visualize raw + detection
    if size(Time_axis_p,1) > max(size(rowc_p))
        Time_axis_p = Time_axis_p(1:end-1);
    end
    
    if check ==1
        figure; imshow(Ird_p,[]); hold on; plot(rowc_p, Time_axis_p(1:end)','ro'); hold off
    end
    
    Fil_Lum_p = abs(rowc_p-rowc_p(1))*(0.001*PIX);
    if Depol ==1
        Fil_Lum_p = Fil_Lum_p;
    end
    Fil_Lsub_p = Fil_Lum_p*370;
    Time_p = Time_axis_p(1:end).*(cycle_time);
    
    % Moving average X points
    NP = 21; % Number of points to average
    LengthRM_p = runmean(Fil_Lsub_p, NP);
               
    figure;
    subplot(2,1,1)
    plot(Time_p,Fil_Lsub_p,'.',Time_p,LengthRM_p,'r','LineWidth',2)
    xlabel('time (s)')
    ylabel('\DeltaLength (sub)')
    title(filename)
    
    dF_p = diff(LengthRM_p)./diff(Time_p)';
    [n_p,x_p] = hist(dF_p,floor(min(dF_p)):1:ceil(max(dF_p)));
    P_p = 1- sum((dF_p>Emin)/length(dF_p));
    
    subplot(2,1,2)
    plot(Time_p(2:end),dF_p,'')
    xlabel('time (s)')
    ylabel('Elongation Rate (sub/s)')
    text(1,0.75*max(dF_p),['Pause probability =',num2str(P_p)] )
    
    
    
    figure;
    plot(Time_b,Fil_Lsub_b,'^r')
    hold on;
    plot(Time_b,LengthRM_b,'Color',[0.3 0.3 0.3],'LineWidth',2)
    plot(Time_p,-Fil_Lsub_p,'o','Color',[0.2 0.5 1])
    plot(Time_p,-LengthRM_p,'Color',[0.3 0.3 0.3],'LineWidth',2)
    hold off;
    xlabel('time (s)')
    ylabel('\DeltaLength (sub)')
    title(filename)
end


save(filename)

