%% Figure 6 A, B, D, E, F, H, I, J, L, M, N, P
window = 1000;
calcdF = @(F, B) (F-B)./medfilt1(F, window, [], 2, 'truncate', 'omitnan');
subdFF = @(F, B) calcdF(F, B)-median(calcdF(F, B), 2); 
folders{1} = '..\datasets\videos\registered\Nestor-Kalinoski\';
folders{2} = '..\datasets\videos\registered\Drumm\';
folders{3} = '..\datasets\videos\registered\Chevalier\';
folders{4} = '..\datasets\videos\registered\Boesmans\';
f1 = figure;
tg1 = uitabgroup;
f2 = figure;
f2.Position(1) = f2.Position(1) - f2.Position(3)/2;
f2.Position(3) = f2.Position(3)*2;
tg2 = uitabgroup;
for fold = 1:length(folders)
    A = dir(strcat(folders{fold},'*registered.h5'));
    it = 3;
    prC = [];
    poC = [];
    for fn = 1:length(A)
        vidname = A(fn).name;
        load(strcat(folders{fold},vidname(1:end-3), '_associated.mat'))
        unmvid1 = h5read(strcat(folders{fold},vidname), '/video3');
        itmove = sum(maxmove((it+1), :),1);
        unmvid1 = unmvid1(itmove(3)+1:end-itmove(4), itmove(1)+1:end-itmove(2), :);
        totalmove = sum(maxmove(1:(it+1), :), 1);
        unmFOV = zeros(d1(1), d2(2));
        unmFOV(totalmove(3)+1:end-totalmove(4), totalmove(1)+1:end-totalmove(2)) = 1;
        oind = strfind(vidname, 'cubic');
        ogname = strcat(vidname(1:(oind-1)), '.h5');
        ogvid = h5read(strcat(folders{fold}, ogname), '/video');
        [row, col] = find(unmFOV == 1);
        r = [min(row) max(row)];
        c = [min(col) max(col)];
        cog = ogvid(r(1):r(2), c(1):c(2), :);
        [rC{fn},mM2f,vM2f] = motion_metrics(unmvid1, 0);
        [oC{fn},mM2f,vM2f] = motion_metrics(cog, 0);
        prC = cat(1, prC, rC{fn});
        poC = cat(1, poC, oC{fn});
        frames = min([120 size(unmvid1, 3)]);

        thistab = uitab(tg2);
        axes('Parent',thistab);

        subplot(1, 2, 1)
        imagesc(max(unmvid1(:, :, 1:frames), [], 3))

        axis equal
        title('Registered max projection')
        subplot(1, 2, 2)
        imagesc(max(cog(:, :, 1:frames), [], 3))
        title('Unregistered max projection')
        axis equal
        sgtitle(ogname)
        name_parts = split(folders{fold},'\');
        thistab.Title = vidname(1:strfind(vidname,"_okada1cubic_registered.h5")-1);
    end
    figure(f1)
    thistab = uitab(tg1);
    thistab.Title = name_parts{end-1};
    axes('Parent',thistab);

    step = 0.02;
    [oofN,oofXedges,oofYedges] = histcounts2(poC, prC, -1:step:1.1, -1:step:1.1, 'Normalization', 'probability');
    ax1=axes;
    hold on
    plot([-1 1], [-1 1], 'k--')
    ax2=axes;
    [M h] = contour(oofXedges(1:end-1)+step/2, oofYedges(1:end-1)+step/2,oofN', 4, 'LineWidth', 1.5);
    ax2.Visible='off';

    xlim([-1 1])
    ylim([-1 1])
    xticks([-1:1])
    yticks([-1:1])

    linkaxes([ax1 ax2]) % this bit is critical, it makes sure that the limits match up
    cvec=linspace(0.4,1,255)';
    z=zeros(255, 1);
    colormap(ax2,[cvec z z])

    xlabel(ax1, 'Unregistered correlation')
    ylabel(ax1, 'Registered correlation')

    figure(f2)
    thistab = uitab(tg2);
    thistab.Title = name_parts{end-1};
    axes('Parent',thistab);

    plot(poC, prC, 'k.')
    hold on
    plot([-1 1], [-1 1], 'k--')
    xlabel('Unregistered correlation')
    ylabel('Registered correlation')
    xlim([-1 1])
    ylim([-1 1])
    xticks([-1:1])
    yticks([-1:1])
end
close(f1)

%% Figure 6O
figure
video = h5read('..\datasets\videos\registered\Chevalier\Chevalier-v3_okada1cubic_registered.h5', '/video3');
load('..\datasets\videos\registered\Chevalier\Chevalier-v3_okada1cubic_registered_associated.mat')
videoc = video(maxmove(4, 3)+1:(end-maxmove(4, 4)), maxmove(4, 1)+1:(end-maxmove(4, 2)), :);
imshowpair(video(:, :, 1), videoc(:, :, 20))

%% Figure 6 C G K
pooledchange = cell(3, 1);
folind = [1 2 4];
for nn = 1:3
    fold = folind(nn);
    mat = dir(strcat(folders{fold},'*Fluorescence.mat'));
    for fn = 1:length(mat)
        figure
        load(strcat(folders{fold},mat(fn).name))
        nframe = size(F, 2);
        startframe = min([50 nframe*0.1]);
        Fr = F;
        Br = B;
        Frmed = median(Fr,2);
        [d1 d2] = size(Fr);
        dFFr = (Fr-Br)./repmat(Frmed, 1, d2);
        dFFr = dFFr-median(dFFr(:, 1:startframe)')';
        dFFrmax = max(dFFr');
        [B dFs] = sort(dFFrmax);
        matu = dir(strcat(folders{fold},mat(fn).name(1:end-4), '_unregistered.mat'));
        load(strcat(folders{fold},matu(1).name))
        Fmed = median(F,2);
        [d1 d2] = size(F);
        dFF = (F-B)./repmat(Fmed, 1, d2);
        dFF = dFF-median(dFF(:, 1:startframe)')';
        pooledchange{nn} = cat(1, pooledchange{nn}, median(Fr-F, 2)./Fmed);
        imshowpair(dFFr(dFs, :), dFF(dFs, :), 'Scaling', 'joint')
        axis normal
        uind = strfind(mat(fn).name, '_');
        title(strcat(mat(fn).name(1:(uind-1)), 32, 'Green registered magenta unregistered'))
        ylabel('Cell number')
        xlabel('Frame number')
        xlim([0 size(F, 2)])
    end
end