Fluordir = '..\datasets\videos\registered\';
A = dir(strcat(Fluordir,'*Fluorescence.mat'));

%% Plot example trajectories form the 35 registered UNR videos.
% D is from row 3 column 5
count = 1;
for fn = 1:length(A)
    load(strcat(Fluordir, A(fn).name), 'xTraj', 'yTraj')
    if ~isempty(xTraj)
        subplot(7, 5, count) 
        t = (1:size(xTraj, 2))*0.05;
        plot(t, sqrt((xTraj-xTraj(:, 1)).^2+(yTraj-yTraj(:, 1)).^2)')
        xlabel('Times (s)')
        ylabel('Displacement (pixels)')
        count = count+1;
    end
end

%% Plot the components of the simulation for subpanels A-C
figdir = '..\datasets\figure_1-3\';
simdir = '..\datasets\simulations\';
simdata = strcat(figdir, 'simulated data_13-Jul-2022-03-01-37Amplitude131nGanglia27nNeuron100seed3.mat');
regdata = strcat(simdir,'050720-10x-2xbin-Mp50-chatgc3-%5B2.041-0.06%5Defs010_Green_background13-Jul-2022-03-01-37Amplitude131nGanglia27nNeuron100seed3uint_okada1cubic_registered_associated.mat');
load(simdata)
video = h5read(strcat(simdir,'050720-10x-2xbin-Mp50-chatgc3-%5B2.041-0.06%5Defs010_Green_background13-Jul-2022-03-01-37Amplitude131nGanglia27nNeuron100seed3uint8.h5'), '/video');
rgang = 2; 
mask = newgang{rgang};
[row, col] = find(sum(mask, 3)); 
r = [min(row) max(row)];
c = [min(col) max(col)];
smask = sum(mask, 3);
figure
subplot(1, 3, 1)
imagesc(smask(r(1):r(2), c(1):c(2)))
title('Single ganglion')
axis equal
colormap gray
subplot(1, 3, 2)
imagesc(svid(:, :, 1))
title('Randomly generated lattice with internodal strands')
axis equal
colormap gray
subplot(1, 3, 3)
imagesc(video(:, :, 1))
title('Assembled field of view')
axis equal
colormap gray

%% Plot the example simulated trajectory for subpanel E
figure
load(strcat(Fluordir, '031220-10x-2xbin-p40-chatgc3-%5B2.041-0.06%5DA30stim006_Green_okada1cubic_registered_Fluorescence.mat'))
simtraj = sqrt((saveMove(:, :, 1)-saveMove(:, 1, 1)).^2+(saveMove(:, :, 2)-saveMove(:, 1, 2)).^2);
t = (1:length(simtraj))*0.05;
plot(t, simtraj)
xlabel('Time (s)')
ylabel('Displacement (pixels)')
ylim([0 180])
title('Simulated movement')

%% Generate the plot and data used for subpanel G and I
[X, Y, eucdist, MSEeu, timeoff, maxErr] = validateSimulationRegistration(simdata, regdata);
whitered = ones(255, 3);
whitered(1:255, 2) = linspace(1, 0, 255);
whitered(1:255, 3) = linspace(1, 0, 255);
figure
imagesc(eucdist', [0 20])
xlim([0 2400])
xticks(linspace(0, 2400, 7))
xticklabels(num2str(linspace(0, 2400, 7)'*0.05))
xlabel('Time (s)')
ylabel('Cell number')
title('Euclidean distance from ground truth')
colormap(whitered)
colorbar

%Subpanel F
rvideo = h5read(strcat(figdir,'050720-10x-2xbin-Mp50-chatgc3-%5B2.041-0.06%5Defs010_Green_background13-Jul-2022-03-01-37Amplitude131nGanglia27nNeuron100seed3uint_okada1cubic_registered.h5'), '/video3');
load(regdata)
mmove = sum(maxmove(1:3, :), 1);
cvideo = video(mmove(3)+1:end-mmove(4), mmove(1)+1:end-mmove(2), :);
figure
imshowpair(rvideo(:, :, 1123), rvideo(:, :, 1));

%% Plot for K and L, the neuron used in the paper is neuron 8
load(strcat(simdir,'050720-10x-2xbin-Mp50-chatgc3-%5B2.041-0.06%5Defs010_Green_background13-Jul-2022-03-01-37Amplitude131nGanglia27nNeuron100seed3uint_okada1cubic_registered_Fluorescence.mat'))
figure
for fn = 8
    subplot(2, 1, 1)
    t = (1:size(F, 2))*0.05;
    plot(t, Ftr(fn, :), 'k')
    hold on
    plot(t, F(fn, :), 'r')
    legend({'Ground Truth', 'NoRMCorre'})
    yv = ylim;
    ylim([0 ceil(yv(2)/5)*5])
    xlabel('Time (s)')
    ylabel('Fluorescence')
    subplot(2, 1, 2)
    plot(t, Ftr(fn, :)-F(fn, :), 'k')
    hold on
    yline(mad(Ftr(fn, :)', 1)*1.48, 'k--')
    yv = ylim;
    ylim([floor(yv(1)/5)*5 ceil(yv(2)/5)*5])
    xlabel('Time (s)')
    ylabel('GT - NoRMCorre Fluorescence')
    legend({'GT - NoRMCorre', 'Estimated sigma'})
    sgtitle(num2str(fn))
end

%% Panels H and J
% Look into the aggregate results for the validation of NoRMCorre on the
% simulated data
figure
clearvars -except figdir
load(strcat(figdir, '20220809_simulationValidation.mat'))
Amplitude_array=[45 83 131 188];
nGang_array=[5 16 23 27 35];
nNeuron_array=[1 5 8 12 100];

for fn = 1:length(simVid)
    IDv = simVid(fn).name;
    Aind = strfind(IDv, 'Amplitude');
    Gind = strfind(IDv, 'nGanglia');
    Nind = strfind(IDv, 'nNeuron');
    sind = strfind(IDv, 'seed');
    
    amp(fn) = str2num(IDv((Aind+9):(Gind-1)));
    nG(fn) = str2num(IDv((Gind+8):(Nind-1)));
    nN(fn) = str2num(IDv((Nind+7):(sind-1)));
    s(fn) = str2num(IDv(sind+4));
    vidID(fn) = str2num(strcat(num2str([amp(fn) nG(fn) nN(fn) s(fn)])')');
    totalN(fn) = length(MSEeu{fn});
    if fn <= length(MSEeu)
        avgEu(fn) = mean(MSEeu{fn});
        avgtime(fn) = mean(timeoff{fn});
        if ~isempty(eucdist{fn})
            avgbytime(fn, :) = nanmean(eucdist{fn}, 2);
        else
            avgbytime(fn, :) = NaN;
        end
    else
        avgEu(fn) = NaN;
        avgtime(fn) = NaN;
        avgbytime(fn, :) = NaN;
    end
end

[~, ia, ~] = unique(vidID);
nonunique = 1:length(vidID);
nonunique(ia) = [];
tamp = unique(amp);
c = 1;
for fn = 1:length(tamp)
    ind = find(amp == tamp(fn));
    ind(ismember(ind, nonunique)) = [];
    subplot(2, 2, fn)
    t = (1:size(avgbytime, 2))*0.05;
    plot(t, nanmean(avgbytime(ind, :)), 'r')
    hold on
    for trace = 1:length(ind)
    plot(t, avgbytime(ind(trace), :)', 'Color', [0 0 1]*(totalN(ind(trace))/max(totalN)))
    if trace == 1
        legend({'Pooled error', 'Video error'}, 'AutoUpdate', 'off')
    end
    if all(diff(avgbytime(ind(trace), :)) == 0)
        bad(c) = ind(trace);
        c = c+1;
    end
    end
    plot(t, nanmean(avgbytime(ind, :)), 'r')
    xlabel('Time (s)')
    ylabel('Displacement (pixels)')
    title(strcat('Maximum displacement =', 32, num2str(tamp(fn)), 32, '(pixels)'))
    ylim([0 150])
end

A = cat(1, amp, nG, nN, s, avgEu, avgtime);

whitered = ones(255, 3);
whitered(1:255, 2) = linspace(1, 0, 255);
whitered(1:255, 3) = linspace(1, 0, 255);
    
f1 = figure;
f2 = figure;
square = [];
squaretime= [];
ind = find(nN == 100);
for n = 1:length(Amplitude_array)
    for g = 1:length(nGang_array)
        sqind = ind(find(amp(ind) == Amplitude_array(n) & nG(ind) == nGang_array(g)));
        sqind(ismember(sqind, nonunique)) = [];
        square(n, g) = nanmean([MSEeu{sqind}]);
        squaretime(n, g) = nanmean([timeoff{sqind}]);
        neuroninbox(n, g) = sum(totalN(sqind));
    end
end
figure(f1)
subplot(1, 2, 1)
imagesc(square, [0 15])
yticks(1:length(Amplitude_array))
xticks(1:length(nGang_array))
yticklabels(num2str(Amplitude_array'))
xticklabels(num2str(nGang_array'))
ylabel('Maximum movement amplitude (pixels')
xlabel('Number of Ganglia in field of view')
title('100 maximum neurons in ganglia')
colormap(whitered)
colorbar
for n = 1:length(Amplitude_array)
    for g = 1:length(nGang_array)
        text(g-0.2, n-0.1, strcat(num2str(round(square(n, g), 1))))
        text(g-0.2, n+0.1, strcat(num2str(round(neuroninbox(n, g), 1))))
    end
end
figure(f2)
subplot(1, 2, 1)
imagesc(squaretime, [0 40])
yticks(1:length(Amplitude_array))
xticks(1:length(nGang_array))
yticklabels(num2str(Amplitude_array'))
xticklabels(num2str(nGang_array'))
ylabel('Maximum movement amplitude (pixels')
xlabel('Number of Ganglia in field of view')
title('100 maximum neurons in ganglia')
for n = 1:length(Amplitude_array)
    for g = 1:length(nGang_array)
        text(g-0.2, n-0.1, strcat(num2str(round(squaretime(n, g), 1))))
        text(g-0.2, n+0.1, strcat(num2str(round(neuroninbox(n, g), 1))))
    end
end
colormap(whitered)
colorbar

square = [];
squaretime= [];
ind = find(nG == 23);
for n = 1:length(Amplitude_array)
    for g = 1:length(nNeuron_array)
        sqind = ind(find(amp(ind) == Amplitude_array(n) & nN(ind) == nNeuron_array(g)));
        sqind(ismember(sqind, nonunique)) = [];
        square(n, g) = nanmean([MSEeu{sqind}]);
        squaretime(n, g) = nanmean([timeoff{sqind}]);
        neuroninbox(n, g) = sum(totalN(sqind));
    end
end
figure(f1)
subplot(1, 2, 2)
imagesc(square, [0 15])
yticks(1:length(Amplitude_array))
xticks(1:length(nNeuron_array))
yticklabels(num2str(Amplitude_array'))
xticklabels(num2str(nNeuron_array'))
colormap(whitered)
ylabel('Maximum movement amplitude (pixels')
xlabel('Maximum number of neurons in ganglia')
colorbar
title('23 ganglia in field of view')
sgtitle('Average error')
for n = 1:length(Amplitude_array)
    for g = 1:length(nGang_array)
        text(g-0.2, n-0.1, strcat(num2str(round(square(n, g), 1))))
        text(g-0.2, n+0.1, strcat(num2str(round(neuroninbox(n, g), 1))))
    end
end
figure(f2)
subplot(1, 2, 2)
imagesc(squaretime, [0 40])
yticks(1:length(Amplitude_array))
xticks(1:length(nNeuron_array))
yticklabels(num2str(Amplitude_array'))
xticklabels(num2str(nNeuron_array'))
colormap(whitered)
ylabel('Maximum movement amplitude (pixels')
xlabel('Maximum number of neurons in ganglia')
colorbar
title('23 ganglia in field of view')
sgtitle('Average time spent > 8 pixels off')
for n = 1:length(Amplitude_array)
    for g = 1:length(nGang_array)
        text(g-0.2, n-0.1, strcat(num2str(round(squaretime(n, g), 1))))
        text(g-0.2, n+0.1, strcat(num2str(round(neuroninbox(n, g), 1))))
    end
end

%% Fluorescence comparison for panel M
clearvars -except figdir
load(strcat(figdir, '20220809_simulationValidation.mat'))
fluodir = '..\datasets\simulations\';
Seed_array=[0 1 2 3];
Amplitude_array=[45 83 131 188];
nGang_array=[5 16 23 27 35];
nNeuron_array=[1 5 8 12 100];

RMSE = cell(length(simVid), 1);
for fn = 1:length(simVid)
    vidname = simVid(fn).name;
    uind = strfind(vidname, '_');
    fname = dir(strcat(fluodir, '*', vidname(uind+1:end-4),'*Fluorescence.mat'));
    if ~isempty(fname) && length(fname) == 1
        IDv = fname(1).name;
        Aind = strfind(IDv, 'Amplitude');
        Gind = strfind(IDv, 'nGanglia');
        Nind = strfind(IDv, 'nNeuron');
        sind = strfind(IDv, 'seed');
        
        amp(fn) = str2num(IDv((Aind+9):(Gind-1)));
        nG(fn) = str2num(IDv((Gind+8):(Nind-1)));
        nN(fn) = str2num(IDv((Nind+7):(sind-1)));
        s(fn) = str2num(IDv(sind+4));
        vidIDF(fn) = str2num(strcat(num2str([amp(fn) nG(fn) nN(fn) s(fn)])')');
        load(strcat(fluodir, fname(1).name))
        totalNF(fn) = size(F, 1);
        dFF = (F'-Ftr')./repmat(median(Ftr, 2,'omitnan'), 1, size(F, 2))';
        dFFtr = (Ftr'-repmat(median(Ftr, 2,'omitnan'), 1, size(F, 2))')./repmat(median(Ftr, 2,'omitnan'), 1, size(F, 2))';
        RMSE{fn} = sqrt(sum((dFF).^2)/length(F));
        FMAD{fn} = mad(dFFtr, 1);
    else
        amp(fn) = NaN;
        nG(fn) = NaN;
        nN(fn) = NaN;
        s(fn) = NaN;
        vidIDF(fn) = NaN;
        totalNF(fn) = NaN;
        RMSE{fn} = NaN;
        FMAD{fn} = NaN;
        
    end
end

[C, ia, ic] = unique(vidIDF);
nonunique = 1:length(vidIDF);
nonunique(ia) = [];



whitered = ones(255, 3);
whitered(1:255, 2) = linspace(1, 0, 255);
whitered(1:255, 3) = linspace(1, 0, 255);

f3 = figure;
f4 = figure;
square = [];
squaretime= [];
ind = find(nN == 100);
for n = 1:length(Amplitude_array)
    for g = 1:length(nGang_array)
        sqind = ind(find(amp(ind) == Amplitude_array(n) & nG(ind) == nGang_array(g)));
        sqind(ismember(sqind, nonunique)) = [];
        square(n, g) = nanmean([RMSE{sqind}]./([FMAD{sqind}]*1.48));
        squaretime(n, g) = nanmean([RMSE{sqind}]);
        neuroninbox(n, g) = sum(totalNF(sqind));
    end
end
figure(f3)
subplot(1, 2, 1)
imagesc(square, [0 10])
yticks(1:length(Amplitude_array))
xticks(1:length(nGang_array))
yticklabels(num2str(Amplitude_array'))
xticklabels(num2str(nGang_array'))
ylabel('Maximum movement amplitude (pixels')
xlabel('Number of Ganglia in field of view')
title('100 maximum neurons in ganglia')
colormap(whitered)
colorbar
for n = 1:length(Amplitude_array)
    for g = 1:length(nGang_array)
        text(g-0.2, n-0.1, strcat(num2str(round(square(n, g), 2))))
        text(g-0.2, n+0.1, strcat(num2str(round(neuroninbox(n, g), 2))))
    end
end
figure(f4)
subplot(1, 2, 1)
imagesc(squaretime, [0 0.5])
yticks(1:length(Amplitude_array))
xticks(1:length(nGang_array))
yticklabels(num2str(Amplitude_array'))
xticklabels(num2str(nGang_array'))
ylabel('Maximum movement amplitude (pixels')
xlabel('Number of Ganglia in field of view')
title('100 maximum neurons in ganglia')
for n = 1:length(Amplitude_array)
    for g = 1:length(nGang_array)
        text(g-0.2, n-0.1, strcat(num2str(round(squaretime(n, g), 2))))
        text(g-0.2, n+0.1, strcat(num2str(round(neuroninbox(n, g), 2))))
    end
end
colormap(whitered)
colorbar

square = [];
squaretime= [];
ind = find(nG == 23);
for n = 1:length(Amplitude_array)
    for g = 1:length(nNeuron_array)
        sqind = ind(find(amp(ind) == Amplitude_array(n) & nN(ind) == nNeuron_array(g)));
        sqind(ismember(sqind, nonunique)) = [];
        square(n, g) = nanmean([RMSE{sqind}]./([FMAD{sqind}]*1.48));
        squaretime(n, g) = nanmean([RMSE{sqind}]);
        neuroninbox(n, g) = sum(totalNF(sqind));
    end
end
figure(f3)
subplot(1, 2, 2)
imagesc(square, [0 10])
yticks(1:length(Amplitude_array))
xticks(1:length(nNeuron_array))
yticklabels(num2str(Amplitude_array'))
xticklabels(num2str(nNeuron_array'))
colormap(whitered)
ylabel('Maximum movement amplitude (pixels')
xlabel('Maximum number of neurons in ganglia')
colorbar
title('23 ganglia in field of view')
sgtitle('$\frac{RMSE}{Estimated\ \sigma}$', 'Interpreter', 'latex')

for n = 1:length(Amplitude_array)
    for g = 1:length(nGang_array)
        text(g-0.2, n-0.1, strcat(num2str(round(square(n, g), 2))))
        text(g-0.2, n+0.1, strcat(num2str(round(neuroninbox(n, g), 2))))
    end
end
figure(f4)
subplot(1, 2, 2)
imagesc(squaretime, [0 0.5])
yticks(1:length(Amplitude_array))
xticks(1:length(nNeuron_array))
yticklabels(num2str(Amplitude_array'))
xticklabels(num2str(nNeuron_array'))
colormap(whitered)
ylabel('Maximum movement amplitude (pixels')
xlabel('Maximum number of neurons in ganglia')
colorbar
title('23 ganglia in field of view')
sgtitle('RMSE of raw fluorescence')
for n = 1:length(Amplitude_array)
    for g = 1:length(nGang_array)
        text(g-0.2, n-0.1, strcat(num2str(round(squaretime(n, g), 2))))
        text(g-0.2, n+0.1, strcat(num2str(round(neuroninbox(n, g), 2))))
    end
end