hmfolder = pwd;
segcorr = "..\datasets\figure_1-3\focus-data-15-videos.mat";
fluoroof = "..\datasets\figure_1-3\050820-10x-2xbin-Mp50-chatgc3-%5B2.041-0.06%5Defs021_Green_okada1cubic_registered_comeintofocus_Fluorescence.mat";
fluorconsensus = "050820-10x-2xbin-Mp50-chatgc3-%5B2.041-0.06%5Defs021_Green_okada1cubic_registered_Fluorescence.mat";

framestemp = 30:35;
framesoof = 850:855;
framesinf = 770:775;
oofind = [96 97 99 100];
framerate = 20;

infocusmaskfolder = '..\datasets\masks\';
cd(infocusmaskfolder)
load('050820-10x-2xbin-Mp50-chatgc3-%5B2.041-0.06%5Defs021_Green_appliedRegistration_projection_consensusmasks.mat');
load('050820-10x-2xbin-Mp50-chatgc3-%5B2.041-0.06%5Defs021_Green_okada1cubic_registered_associated.mat', 'maxmove');

infmaxmove = maxmove;
focusmasks = masks; 


focusfolder = '..\datasets\masks\consensus-labels\';
cd(hmfolder)
cd(focusfolder)
load('050820-10x-2xbin-Mp50-chatgc3-%5B2.041-0.06%5Defs021_Green_appliedRegistration.mat')

vidfolder = '..\datasets\videos\registered\';
cd(hmfolder)
cd("..\datasets\figure_1-3\")
load('050820-10x-2xbin-Mp50-chatgc3-%5B2.041-0.06%5Defs021_Green_okada1cubic_registered_comeintofocus_masks.mat')
cd(hmfolder)
cd(vidfolder)
load('050820-10x-2xbin-Mp50-chatgc3-%5B2.041-0.06%5Defs021_Green_okada1cubic_registered_associated.mat', 'maxmove');
video = h5read('050820-10x-2xbin-Mp50-chatgc3-%5B2.041-0.06%5Defs021_Green_okada1cubic_registered.h5', '/video3');

oofmasks = masks;
oofmasks = fliplr(permute(flipud(oofmasks), [2 1 3]));

lims = [0 100];
[~, video, ~] = convergeMaxMove(focusmasks, video, infmaxmove, maxmove);
[focusmasks, oofmasks, ~] = convergeMaxMove(focusmasks, oofmasks, infmaxmove, maxmove);
subplot(1, 3, 1)
imagesc(mean(video(:, :, framestemp), 3), lims)
title('Template image')

colormap gray
axis equal
subplot(1, 3, 2)
mask2image_contour(focusmasks(:, :, oofind), mean(video(:, :, framesoof), 3), lims)
title('Out of focus neurons')
axis equal
subplot(1, 3, 3)
mask2image_contour(oofmasks, mean(video(:, :, framesinf), 3), lims)
title('New in focus neurons')
axis equal

figure
subplot(2, 1, 1)
load(fluorconsensus)
dFF = (F-repmat(median(F, 2),1 , size(F, 2)))./repmat(median(F, 2),1 , size(F, 2));
plotOpticalTraces(dFF', oofind, framerate)
yv = ylim;
hold on
for fn = 1:size(focusTable, 1)
    if focusTable(fn, 6) == 1
        plot([focusTable(fn, 1) focusTable(fn, 2)]/framerate, [yv(2) yv(2)], 'k', 'Linewidth', 2)
    end
end
xlabel('Time (s)')
ylabel('dFF')
title('Neurons that go out of focus')

subplot(2, 1, 2)

cd(hmfolder)
load(fluoroof)
dFF = (F-repmat(median(F, 2),1 , size(F, 2)))./repmat(median(F, 2),1 , size(F, 2));
plotOpticalTraces(dFF', 1:size(F, 1), framerate)
yv = ylim;
hold on
for fn = 1:size(focusTable, 1)
    if focusTable(fn, 3) == 1
        plot([focusTable(fn, 1) focusTable(fn, 2)]/framerate, [yv(2) yv(2)], 'k', 'Linewidth', 2)
    end
end
xlabel('Time (s)')
ylabel('dFF')
title('Neurons that go into focus')

%% Second example
vidname = '030220-10x-2xbin-p60-mnnosgc3-%5B2.041-0.06%5DA30stim008_Green_okada1cubic_registered.h5';
fluorconsensus = '030220-10x-2xbin-p60-mnnosgc3-%5B2.041-0.06%5DA30stim008_Green_okada1cubic_registered_Fluorescence.mat';
framestemp = 30:35;
framesoof = 875:880;
framesinf = 1210:1215;
oofind = [16 17 18 19];
framerate = 20;

cd(hmfolder)
cd(infocusmaskfolder)
load('030220-10x-2xbin-p60-mnnosgc3-%5B2.041-0.06%5DA30stim008_Green_appliedRegistration_projection_consensusmasks.mat')
load('030220-10x-2xbin-p60-mnnosgc3-%5B2.041-0.06%5DA30stim008_Green_okada1cubic_registered_associated.mat', 'maxmove');
infmaxmove = maxmove;
focusmasks = masks; 

cd(hmfolder)
cd(focusfolder)
load('030220-10x-2xbin-p60-mnnosgc3-%5B2.041-0.06%5DA30stim008_Green_appliedRegistration.mat')
cd(hmfolder)
cd("..\datasets\figure_1-3\")
load('030220-10x-2xbin-p60-mnnosgc3-%5B2.041-0.06%5DA30stim008_Green_okada1cubic_registered_comeintofocus_masks.mat')
cd(hmfolder)
cd(vidfolder)
load('030220-10x-2xbin-p60-mnnosgc3-%5B2.041-0.06%5DA30stim008_Green_okada1cubic_registered_associated.mat', 'maxmove');
oofmasks = masks;
oofmasks = fliplr(permute(flipud(oofmasks), [2 1 3]));
video = h5read(vidname, '/video3');
lims = [0 100];
[focusmasks_rs, video, ~] = convergeMaxMove(focusmasks, video, infmaxmove, maxmove);
[focusmasks, oofmasks, ~] = convergeMaxMove(focusmasks, oofmasks, infmaxmove, maxmove);
subplot(1, 3, 1)
imagesc(mean(video(:, :, framestemp), 3), lims)
title('Template image')

colormap gray
axis equal
subplot(1, 3, 2)
mask2image_contour(focusmasks(:, :, oofind), mean(video(:, :, framesoof), 3), lims)
title('Out of focus neurons')
axis equal
subplot(1, 3, 3)
mask2image_contour(oofmasks, mean(video(:, :, framesinf), 3), lims)
title('New in focus neurons')
axis equal

figure
subplot(2, 1, 1)
load(fluorconsensus)
dFF = (F-repmat(median(F, 2),1 , size(F, 2)))./repmat(median(F, 2),1 , size(F, 2));
plotOpticalTraces(dFF', oofind, framerate)
yv = ylim;
hold on
for fn = 1:size(focusTable, 1)
    if focusTable(fn, 6) == 1
        plot([focusTable(fn, 1) focusTable(fn, 2)]/framerate, [yv(2) yv(2)], 'k', 'Linewidth', 2)
    end
end
xlabel('Time (s)')
ylabel('dFF')
title('Neurons that go out of focus')

subplot(2, 1, 2)

cd(hmfolder)
load(fluoroof)
dFF = (F-repmat(median(F, 2),1 , size(F, 2)))./repmat(median(F, 2),1 , size(F, 2));
plotOpticalTraces(dFF', 1:size(F, 1), framerate)
yv = ylim;
hold on
for fn = 1:size(focusTable, 1)
    if focusTable(fn, 4) == 1
        plot([focusTable(fn, 1) focusTable(fn, 2)]/framerate, [yv(2) yv(2)], 'k', 'Linewidth', 2)
    end
end
xlabel('Time (s)')
ylabel('dFF')
title('Neurons that go into focus')




load(segcorr)
oofcorr = [];
infcorr = [];
oofssim = [];
infssim = [];
medlength = 80;
for vidind = 1:15
correlation = focus_layers(vidind).correlation;
ssim = focus_layers(vidind).ssim;
focusTable = focus_layers(vidind).focusTable;
quad = focus_layers(vidind).quadrant; 

order = [2 1 3 4];
oofind = cell(4, 1);
for period = 1:size(focusTable, 1)
    for q = 1:4
        if focusTable(period, q+2) == 1
            oofind{q} = cat(2, oofind{q}, focusTable(period, 1):focusTable(period, 2));
        end
    end
end
for ganglia = 1:size(correlation, 1)
    if quad(ganglia) ~= 0 && ~isnan(quad(ganglia))
    corrt = movmedian(normalize(correlation(ganglia, :),'center','median'),medlength,'omitnan');
    oofcorrelation = corrt(oofind{quad(ganglia)});
    corrt(oofind{quad(ganglia)}) = [];
    ssimt =  movmedian(normalize(ssim(ganglia, :),'center','median'),medlength,'omitnan');
    oofssimt = ssimt(oofind{quad(ganglia)});
    ssimt(oofind{quad(ganglia)}) = [];
    oofcorr = cat(2, oofcorr, oofcorrelation);
    infcorr = cat(2, infcorr, corrt);
    oofssim = cat(2, oofssim, oofssimt);
    infssim = cat(2, infssim, ssimt);
    end
end


end
    
[oofN,oofXedges,oofYedges] = histcounts2(oofcorr, oofssim, -2:0.01:2, -2:0.01:2, 'Normalization', 'probability');
[infN,infXedges,infYedges] = histcounts2(infcorr, infssim, -2:0.01:2, -2:0.01:2, 'Normalization', 'probability');
B0 = -1.08; %From the logistic regression fit
B1 = -19.848;
B2 = -8.6055;
T = 0.39;

figure
ax1=axes;
contour(oofXedges(1:end-1), oofYedges(1:end-1),oofN', 5, 'LineWidth', 1.5);
hold on
X1 = -1:0.1:1;
plot(X1, (1/B2)*(log(T/(1-T))-B0-B1*X1), 'k')
ax2=axes;
hold on
contour(infXedges(1:end-1), infYedges(1:end-1),infN', 5, 'LineWidth', 1.5);
ax2.Visible='off';



linkaxes([ax1 ax2]) % this bit is critical, it makes sure that the limits match up
xlim([-0.5 0.5])
ylim([-0.5 0.5])
xticks(-0.5:0.5:0.5)
yticks(-0.5:0.5:0.5)
cvec=linspace(0.4,1,255)';
cvec2 = linspace(0.4,1,255)';
z=zeros(255, 1);
z2=zeros(255,1);
colormap(ax1,[cvec z z])
colormap(ax2,[z2 z2 cvec])

xlabel(ax1, 'Normalized Correlation')
ylabel(ax1, 'Normalized ssim')
legend(ax1, {'OOF', 'LR'})
legend(ax2, {'INF'}, 'Location', 'South')