%% Load data
cd('file directory')
%% recover data
BWs = dir('BW*.*');
Nfiles = length(BWs);
boo=0;
for i=1:Nfiles
    % loads
    % load pattern
    filename = dir(['BWpattern',num2str(i),'.mat']);
    filename = filename.name;
    load(filename);
    sp = regionprops(BW);
    r = sqrt(sp.Area/pi); % size of the pattern
    r = r+0.1*r; % safety distance
    cpx = sp.Centroid(1); % center pattern x
    cpy = sp.Centroid(2); % center pattern x
    % load cell and golgi centers
    filename = dir(['cell',num2str(i),'.mat']);
    filename = filename.name;
    load(filename);
    % load nucleus centers
    filename = dir(['DIC',num2str(i),'_*.mat']);
    filename = filename.name;
    load(filename);
    Nframe = length(cell);
    if i<12&i~=8 % convert the variable stats
        Nframe2 = length(stats);
        if Nframe2<Nframe
            Nframe=Nframe2;
        end
        for j=1:Nframe2
            center(j).x = stats(j).Centroid(2);
            center(j).y = stats(j).Centroid(1);
        end
        nuc = [];
        nuc = [center.x;center.y];
    else
        nuc = [center.x,center.y];
        nuc = nuc';
        Nframe2 = length(nuc);
        if Nframe2<Nframe
            Nframe=Nframe2;
        end
    end
    
    % calculate the angle and distances
    for j=1:Nframe
        ac(j) = atan2(cell(2,j)-cpy,cell(1,j)-cpx);
        ag(j) = atan2(golgi(2,j)-cpy,golgi(1,j)-cpx);
        an(j) = atan2(nuc(2,j)-cpy,nuc(1,j)-cpx);
        agn(j) = atan2(golgi(2,j)-nuc(2,j),golgi(1,j)-nuc(1,j));
        dc(j) = sqrt((cell(2,j)-cpy)^2+(cell(1,j)-cpx)^2);
        dn(j) = sqrt((nuc(1,j)-cpx)^2+(nuc(2,j)-cpy)^2);
    end
    dcs = smooth(dc,0.5,'loess');
    dns = smooth(dn,0.5,'loess');
    
    % find the escape frame
    FEc = find(dcs>r,1);
    FEn = find(dns>r,1);
        if boo==0%FEn<100
            figure(3)
            clf;
            title(['cell ',num2str(i),' click on the escape RED curve'])
    
            hold on
            plot([1:Nframe],dc,'--c')
            plot([1:Nframe],dn,'--r')
            plot([1:Nframe],dcs,'c')
            plot([1:Nframe],dns,'r')
            plot([1:Nframe],r.*ones(Nframe,1),'k')
            [xi,yi]=ginput(1);
            FEn = round(xi);
        end
        if boo==0%FEc<100
            figure(3)
            clf;
            title(['cell ',num2str(i),' click on the escape CYAN curve'])
    
            hold on
            plot([1:Nframe],dc,'--c')
            plot([1:Nframe],dn,'--r')
            plot([1:Nframe],dcs,'c')
            plot([1:Nframe],dns,'r')
            plot([1:Nframe],r.*ones(Nframe,1),'k')
            [xi,yi]=ginput(1);
            FEc = round(xi);
        end
    %     if isempty(FEc)
    %         figure(3)
    %         clf;
    %         title(['cell ',num2str(i),' click on the escape CYAN curve'])
    %
    %         hold on
    %         plot([1:Nframe],dc,'--c')
    %         plot([1:Nframe],dn,'--r')
    %         plot([1:Nframe],dcs,'c')
    %         plot([1:Nframe],dns,'r')
    %         plot([1:Nframe],r.*ones(Nframe,1),'k')
    %         [xi,yi]=ginput(1);
    %         FEc = round(xi);
    %     end
    %     if isempty(FEn)
    %         figure(3)
    %         clf;
    %         title(['cell ',num2str(i),' click on the escape RED curve'])
    %
    %         hold on
    %         plot([1:Nframe],dc,'--c')
    %         plot([1:Nframe],dn,'--r')
    %         plot([1:Nframe],dcs,'c')
    %         plot([1:Nframe],dns,'r')
    %         plot([1:Nframe],r.*ones(Nframe,1),'k')
    %         [xi,yi]=ginput(1);
    %         FEn = round(xi);
    %     end
%     if isempty(FEn)|FEn<15
%         FEn=15;
%     end
%     if isempty(FEc)|FEc<15
%         FEc=15;
%     end

    % save all required info
    data(i).ac = {radtodeg(ac)};
    data(i).an = {radtodeg(an)};
    data(i).agn = {radtodeg(agn)};
    data(i).dn = {dn};
    data(i).dc = {dc};
    
    data(i).FEc = FEc;
    data(i).FEn = FEn;
    clearvars -except data
    boo=0;
end
save('dataALL.mat','data');


