% Code to assemble the inter-experiment registrations so that all the
% relevant experiments are grouped within a cell

cd('.\datasets\videos\original data')
A = dir('*Green.tif');
count = 1;
group = {};     % Contains the video names that have been cleaned up a bit
sttype = {};    % The stimulation type
vDates = [];    % Dates the data was recorded on
FOV = {};       % The unregistered field of view

g1 = ['13'; '15'; '21']; % IDs for 2 sets of data that were recorded on the same day
g2 = ['11'; '05'; '03'];

% Set up cell arrays for the different days
for fn = 1:length(A)
    ovname = A(fn).name;
    video = bigread2(ovname, 1, 50);
    avgimg = mean(video(:, :, 1:50), 3);
    dind = strfind(ovname, '-');
    vdate = ovname(1:(dind(1)-1));
    tind = strfind(ovname, '30');
    vname = ovname;
    if ~isempty(tind)
     vname((0:1)+tind(end)) = [];
    end

    Dind = strfind(vname, 'D');
    if isempty(Dind)
    else
    zind = strfind(vname(Dind(end):end), '0');
    stim = lower(vname(((1):(zind-2))+Dind));
    if isempty(vDates)
        group{count} = {ovname};
        sttype{count} = {stim};
        FOV{count} = {avgimg};
        vDates = vdate;
    elseif all(vdate == vDates)
        if any(vdate ~= '050820')
            group{count} = [group{count} ovname];
            sttype{count} = [sttype{count} stim];
            FOV{count} = [FOV{count} avgimg];
        else
           for nv = 1:3
               c(nv) = contains(vname(20:end), g2(nv, :));
           end
           if any(c == 1)
              group{g1count}  = [group{g1count} ovname];
              sttype{g1count} = [sttype{g1count} stim];
              FOV{g1count} = [FOV{g1count} avgimg];
           else
               group{count} = [group{count} ovname];
               sttype{count} = [sttype{count} stim];
               FOV{count} = [FOV{count} avgimg];
           end
        end

    elseif any(vdate ~= vDates)
        count = count+1;
        group{count} = {ovname};
        sttype{count} = {stim};
        FOV{count} = {avgimg};
        vDates = vdate; 
        if all(vdate == '050820')
            g1count = count;
            count = count+1;
            group{count} = {};
            sttype{count} = {};
            FOV{count} = {};
        end
        
    end
    end
end

% Do the same thing, but generate cells for the registration parameters and
% ROIS
regdir = 'E:\cluster data\20220707_maxshiftafterestimate_newparameters';
segdir = 'C:\Users\redin\Documents\Segmentation\15+14_segmentation';
mMove = cell(length(group), 1); %Cropping matrix for the registration
rFOV = cell(length(group), 1); %The registered field of view
ROIs = cell(length(group), 1); %the identified ROIs within the field of view
for g = 1:length(group)
    vids = group{g};
    for fn = 1:length(vids)
        cd(regdir)
        V = dir(strcat(vids{fn}(1:end-4), '*_registered.h5'));
        a = h5info(V(1).name);
        dims = a.Datasets(6).Dataspace.Size; 
        video = h5read(V(1).name, '/video3', [1 1 1], [dims(1) dims(2) 100]);
        rIm = mean(video, 3);
        Assoc = dir(strcat(vids{fn}(1:end-4), '*associated.mat'));
        load(Assoc(1).name, 'maxmove');
        rFOV{g} = [rFOV{g} {rIm}];
        mMove{g} = [mMove{g} {maxmove}];
        cd(segdir)
        V = dir(strcat(vids{fn}(1:end-4), '*_masks.mat'));
        if ~isempty(V)
            load(V(1).name)
            ROIs{g} = [ROIs{g} {masks}];
        else
            ROIs{g} = [ROIs{g} {[]}];
        end
    end
end

% Different tests at registration
[optimizer, metric] = imregconfig('multimodal');
optimizer.InitialRadius = 0.001;
optimizer.Epsilon = 1.5e-4;
optimizer.GrowthFactor = 1.01;
optimizer.MaximumIterations = 500;
for fn = 1:length(group)
    num = length(group{fn});
    NCparamStruct.gridsize = [300 300; 150 150; 100 100];
    NCparamStruct.motuf = [5 5 5];
    NCparamStruct.overlap = [64 32 16];
    NCparamStruct.maxshift = [80 25 20];
    NCparamStruct.interp = 'cubic';
    if num < 2
    elseif num == 2
        im1 = FOV{fn}{1};
        im2 = FOV{fn}{2};
        im1 = bandpass2D(im1, 80, 2);
        im2 = bandpass2D(im2, 80, 2);
        [d1, d2, d3] = size(im1);
        M2 = im2;

        for nn = 1:length(NCparamStruct.motuf)
            options_nr = NoRMCorreSetParms('d1',d1,'d2',d2,'bin_width', 50, ...
                'grid_size',NCparamStruct.gridsize(nn, :),'mot_uf',NCparamStruct.motuf(nn),'correct_bidir',false, ...
                'overlap_pre',NCparamStruct.overlap(nn),'overlap_post', NCparamStruct.overlap(nn),'max_shift',....
                NCparamStruct.maxshift(nn), 'upd_template', false, 'shifts_method', NCparamStruct.interp);
            [M2,shifts2{nn+1},template2] = normcorre(M2,options_nr,im1);
        end
        subplot(1, 2, 1)
        imshowpair(im1, im2)
        title('Unregistered')
        subplot(1, 2, 2)
        imshowpair(im1, M2)
        title('Registered')
        pause(0.01)
    elseif num == 3
        im1 = FOV{fn}{1};
        im2 = FOV{fn}{2};
        im1 = bandpass2D(im1, 80, 2);
        im2 = bandpass2D(im2, 80, 2);

        [d1, d2, d3] = size(im1);
        M2 = im2;
        for nn = 1:length(NCparamStruct.motuf)
            options_nr = NoRMCorreSetParms('d1',d1,'d2',d2,'bin_width', 50, ...
                'grid_size',NCparamStruct.gridsize(nn, :),'mot_uf',NCparamStruct.motuf(nn),'correct_bidir',false, ...
                'overlap_pre',NCparamStruct.overlap(nn),'overlap_post', NCparamStruct.overlap(nn),'max_shift',....
                NCparamStruct.maxshift(nn), 'upd_template', false, 'shifts_method', NCparamStruct.interp);
            [M2,shifts2{nn+1},template2] = normcorre(M2,options_nr,im1);
        end
        subplot(1, 2, 1)
        imshowpair(im1, im2)
        title('Unregistered')
        subplot(1, 2, 2)
        imshowpair(im1, M2)
        title('Registered')
        pause(0.01)
        im1 = FOV{fn}{1};
        im2 = FOV{fn}{3};
        im1 = bandpass2D(im1, 80, 2);
        im2 = bandpass2D(im2, 80, 2);

        [d1, d2, d3] = size(im1);
        M2 = im2;
        for nn = 1:length(NCparamStruct.motuf)
            options_nr = NoRMCorreSetParms('d1',d1,'d2',d2,'bin_width', 50, ...
                'grid_size',NCparamStruct.gridsize(nn, :),'mot_uf',NCparamStruct.motuf(nn),'correct_bidir',false, ...
                'overlap_pre',NCparamStruct.overlap(nn),'overlap_post', NCparamStruct.overlap(nn),'max_shift',....
                NCparamStruct.maxshift(nn), 'upd_template', false, 'shifts_method', NCparamStruct.interp);
            [M2,shifts2{nn+1},template2] = normcorre(M2,options_nr,im1);
        end
        subplot(1, 2, 1)
        imshowpair(im1, im2)
        title('Unregistered')
        subplot(1, 2, 2)
        imshowpair(im1, M2)
        title('Registered')
        pause(0.01)
        im1 = FOV{fn}{2};
        im2 = FOV{fn}{3};
        im1 = bandpass2D(im1, 80, 2);
        im2 = bandpass2D(im2, 80, 2);

        [d1, d2, d3] = size(im1);
        M2 = im2;
        for nn = 1:length(NCparamStruct.motuf)
            options_nr = NoRMCorreSetParms('d1',d1,'d2',d2,'bin_width', 50, ...
                'grid_size',NCparamStruct.gridsize(nn, :),'mot_uf',NCparamStruct.motuf(nn),'correct_bidir',false, ...
                'overlap_pre',NCparamStruct.overlap(nn),'overlap_post', NCparamStruct.overlap(nn),'max_shift',....
                NCparamStruct.maxshift(nn), 'upd_template', false, 'shifts_method', NCparamStruct.interp);
            [M2,shifts2{nn+1},template2] = normcorre(M2,options_nr,im1);
        end
        subplot(1, 2, 1)
        imshowpair(im1, im2)
        title('Unregistered')
        subplot(1, 2, 2)
        imshowpair(im1, M2)
        title('Registered')
        pause(0.01)
    end
end


repeats = 2;
Regc = cell(length(group), 1);
regROI = cell(length(group), 1);
regFOV = cell(length(group), 1);
for fn = 1:length(group)
    num = length(group{fn});
    NCparamStruct.gridsize = [600 600; 200 200; 150 150; 50 50];
    NCparamStruct.motuf = [5 5 5 5];
    NCparamStruct.overlap = [64 32 16 8];
    NCparamStruct.maxshift = [80 40 20 15];
    NCparamStruct.interp = 'cubic';
    if num < 2
    elseif num == 2
        im1 = rFOV{fn}{1};
        im2 = rFOV{fn}{2};
        m1 = mMove{fn}{1};
        m2 = mMove{fn}{2};
        [im1, im2, nmaxmove] = convergeMaxMove(im1, im2, m1, m2);

        [d1, d2, d3] = size(im1);
        im1 = mat2gray(im1);
        im2 = mat2gray(im2);
        im1bw = imbinarize(im1,graythresh(im1));
        im2bw = imbinarize(im2,graythresh(im2));
        im1(im1bw == 0) = 0;
        im2(im2bw == 0) = 0;
        M2 = im2;

        for nn = 1:length(NCparamStruct.motuf)
            options_nr = NoRMCorreSetParms('d1',d1,'d2',d2,'bin_width', 50, ...
                'grid_size',NCparamStruct.gridsize(nn, :),'mot_uf',NCparamStruct.motuf(nn),'correct_bidir',false, ...
                'overlap_pre',NCparamStruct.overlap(nn),'overlap_post', NCparamStruct.overlap(nn),'max_shift',....
                NCparamStruct.maxshift(nn), 'upd_template', false, 'shifts_method', NCparamStruct.interp);
            for rep = 1:repeats
                [M2,shifts2{nn, rep},template2] = normcorre(M2,options_nr,im1);
            end
        end
        im1m = ROIs{fn}{1};

        im2m = ROIs{fn}{2};
        [im1m, im2m, nmaxmove] = convergeMaxMove(im1m, im2m, m1, m2);
        im1 = rFOV{fn}{1};
        im2 = rFOV{fn}{2};
        [im1, im2, nmaxmove] = convergeMaxMove(im1, im2, m1, m2);
        for nn = 1:length(NCparamStruct.motuf)
            options_nr = NoRMCorreSetParms('d1',d1,'d2',d2,'bin_width', 50, ...
                'grid_size',NCparamStruct.gridsize(nn, :),'mot_uf',NCparamStruct.motuf(nn),'correct_bidir',false, ...
                'overlap_pre',NCparamStruct.overlap(nn),'overlap_post', NCparamStruct.overlap(nn),'max_shift',....
                NCparamStruct.maxshift(nn), 'upd_template', false, 'shifts_method', NCparamStruct.interp);
            for rep = 1:repeats
                for neuron = 1:size(im2m, 3)
                    im2m(:, :, neuron) = apply_shifts(im2m(:, :,neuron), shifts2{nn, rep}, options_nr);
                end
                im2 = apply_shifts(im2, shifts2{nn, rep}, options_nr);
            end
        end
        Regc{fn} = [Regc{fn} {shifts2}];
        regROI{fn} = [regROI{fn} {{im1m, im2m}}];
        regFOV{fn} = [regFOV{fn} {{im1, im2}}];
        subplot(1, 2, 1)
        imshowpair(sum(im1m, 3), sum(im2m, 3))
        title('Masks')
        subplot(1, 2, 2)
        imshowpair(im1, im2)
        title('FOV ')
        pause(0.01)
    elseif num == 3
        im1 = rFOV{fn}{1};
        im2 = rFOV{fn}{2};
        m1 = mMove{fn}{1};
        m2 = mMove{fn}{2};
        [im1, im2, nmaxmove] = convergeMaxMove(im1, im2, m1, m2);

        im1 = mat2gray(im1);
        im2 = mat2gray(im2);
        im1bw = imbinarize(im1,graythresh(im1));
        im2bw = imbinarize(im2,graythresh(im2));
        im1(im1bw == 0) = 0;
        im2(im2bw == 0) = 0;

        [d1, d2, d3] = size(im1);
        M2 = im2;
        for nn = 1:length(NCparamStruct.motuf)
            options_nr = NoRMCorreSetParms('d1',d1,'d2',d2,'bin_width', 50, ...
                'grid_size',NCparamStruct.gridsize(nn, :),'mot_uf',NCparamStruct.motuf(nn),'correct_bidir',false, ...
                'overlap_pre',NCparamStruct.overlap(nn),'overlap_post', NCparamStruct.overlap(nn),'max_shift',....
                NCparamStruct.maxshift(nn), 'upd_template', false, 'shifts_method', NCparamStruct.interp);
            for rep = 1:repeats
            [M2,shifts2{nn, rep},template2] = normcorre(M2,options_nr,im1);
            end
        end
        im1m = ROIs{fn}{1};

        im2m = ROIs{fn}{2};
        [im1m, im2m, nmaxmove] = convergeMaxMove(im1m, im2m, m1, m2);
        im1 = rFOV{fn}{1};
        im2 = rFOV{fn}{2};
        [im1, im2, nmaxmove] = convergeMaxMove(im1, im2, m1, m2);
        for nn = 1:length(NCparamStruct.motuf)
            options_nr = NoRMCorreSetParms('d1',d1,'d2',d2,'bin_width', 50, ...
                'grid_size',NCparamStruct.gridsize(nn, :),'mot_uf',NCparamStruct.motuf(nn),'correct_bidir',false, ...
                'overlap_pre',NCparamStruct.overlap(nn),'overlap_post', NCparamStruct.overlap(nn),'max_shift',....
                NCparamStruct.maxshift(nn), 'upd_template', false, 'shifts_method', NCparamStruct.interp);
            for rep = 1:repeats
                for neuron = 1:size(im2m, 3)
                    im2m(:, :, neuron) = apply_shifts(im2m(:, :,neuron), shifts2{nn, rep}, options_nr);
                end
                im2 = apply_shifts(im2, shifts2{nn, rep}, options_nr);
            end
        end
        Regc{fn} = [Regc{fn} {shifts2}];
        regROI{fn} = [regROI{fn} {{im1m, im2m}}];
        regFOV{fn} = [regFOV{fn} {{im1, im2}}];
        subplot(1, 2, 1)
        imshowpair(sum(im1m, 3), sum(im2m, 3))
        title('Masks')
        subplot(1, 2, 2)
        imshowpair(im1, im2)
        title('FOV ')
        pause(0.01)
        im1 = rFOV{fn}{1};
        im2 = rFOV{fn}{3};
        m1 = mMove{fn}{1};
        m2 = mMove{fn}{3};
        [im1, im2, nmaxmove] = convergeMaxMove(im1, im2, m1, m2);

        [d1, d2, d3] = size(im1);
        im1 = mat2gray(im1);
        im2 = mat2gray(im2);
        im1bw = imbinarize(im1,graythresh(im1));
        im2bw = imbinarize(im2,graythresh(im2));
        im1(im1bw == 0) = 0;
        im2(im2bw == 0) = 0;
        M2 = im2;
        for nn = 1:length(NCparamStruct.motuf)
            options_nr = NoRMCorreSetParms('d1',d1,'d2',d2,'bin_width', 50, ...
                'grid_size',NCparamStruct.gridsize(nn, :),'mot_uf',NCparamStruct.motuf(nn),'correct_bidir',false, ...
                'overlap_pre',NCparamStruct.overlap(nn),'overlap_post', NCparamStruct.overlap(nn),'max_shift',....
                NCparamStruct.maxshift(nn), 'upd_template', false, 'shifts_method', NCparamStruct.interp);
            for rep = 1:repeats
            [M2,shifts2{nn, rep},template2] = normcorre(M2,options_nr,im1);
            end
        end
         
        im1m = ROIs{fn}{1};

        im2m = ROIs{fn}{3};
        [im1m, im2m, nmaxmove] = convergeMaxMove(im1m, im2m, m1, m2);
        im1 = rFOV{fn}{1};
        im2 = rFOV{fn}{3};
        [im1, im2, nmaxmove] = convergeMaxMove(im1, im2, m1, m2);
        for nn = 1:length(NCparamStruct.motuf)
            options_nr = NoRMCorreSetParms('d1',d1,'d2',d2,'bin_width', 50, ...
                'grid_size',NCparamStruct.gridsize(nn, :),'mot_uf',NCparamStruct.motuf(nn),'correct_bidir',false, ...
                'overlap_pre',NCparamStruct.overlap(nn),'overlap_post', NCparamStruct.overlap(nn),'max_shift',....
                NCparamStruct.maxshift(nn), 'upd_template', false, 'shifts_method', NCparamStruct.interp);
            for rep = 1:repeats
                for neuron = 1:size(im2m, 3)
                    im2m(:, :, neuron) = apply_shifts(im2m(:, :,neuron), shifts2{nn, rep}, options_nr);
                end
                im2 = apply_shifts(im2, shifts2{nn, rep}, options_nr);
            end
        end
        Regc{fn} = [Regc{fn} {shifts2}];
        regROI{fn} = [regROI{fn} {{im1m, im2m}}];
        regFOV{fn} = [regFOV{fn} {{im1, im2}}];
        subplot(1, 2, 1)
        imshowpair(sum(im1m, 3), sum(im2m, 3))
        title('Masks')
        subplot(1, 2, 2)
        imshowpair(im1, im2)
        title('FOV ')
        pause(0.01)
        
        
        im1 = rFOV{fn}{2};
        im2 = rFOV{fn}{3};
        m1 = mMove{fn}{2};
        m2 = mMove{fn}{3};
        [im1, im2, nmaxmove] = convergeMaxMove(im1, im2, m1, m2);

        [d1, d2, d3] = size(im1);
        im1 = mat2gray(im1);
        im2 = mat2gray(im2);
        im1bw = imbinarize(im1,graythresh(im1));
        im2bw = imbinarize(im2,graythresh(im2));
        im1(im1bw == 0) = 0;
        im2(im2bw == 0) = 0;
        M2 = im2;
        for nn = 1:length(NCparamStruct.motuf)
            options_nr = NoRMCorreSetParms('d1',d1,'d2',d2,'bin_width', 50, ...
                'grid_size',NCparamStruct.gridsize(nn, :),'mot_uf',NCparamStruct.motuf(nn),'correct_bidir',false, ...
                'overlap_pre',NCparamStruct.overlap(nn),'overlap_post', NCparamStruct.overlap(nn),'max_shift',....
                NCparamStruct.maxshift(nn), 'upd_template', false, 'shifts_method', NCparamStruct.interp);
            for rep = 1:repeats
            [M2,shifts2{nn, rep},template2] = normcorre(M2,options_nr,im1);
            end
        end
         
        im1m = ROIs{fn}{2};
        im2m = ROIs{fn}{3};
        [im1m, im2m, nmaxmove] = convergeMaxMove(im1m, im2m, m1, m2);
        im1 = rFOV{fn}{2};
        im2 = rFOV{fn}{3};
        [im1, im2, nmaxmove] = convergeMaxMove(im1, im2, m1, m2);
        for nn = 1:length(NCparamStruct.motuf)
            options_nr = NoRMCorreSetParms('d1',d1,'d2',d2,'bin_width', 50, ...
                'grid_size',NCparamStruct.gridsize(nn, :),'mot_uf',NCparamStruct.motuf(nn),'correct_bidir',false, ...
                'overlap_pre',NCparamStruct.overlap(nn),'overlap_post', NCparamStruct.overlap(nn),'max_shift',....
                NCparamStruct.maxshift(nn), 'upd_template', false, 'shifts_method', NCparamStruct.interp);
            for rep = 1:repeats
                for neuron = 1:size(im2m, 3)
                    im2m(:, :, neuron) = apply_shifts(im2m(:, :,neuron), shifts2{nn, rep}, options_nr);
                end
                im2 = apply_shifts(im2, shifts2{nn, rep}, options_nr);
            end
        end
        Regc{fn} = [Regc{fn} {shifts2}];
        regROI{fn} = [regROI{fn} {{im1m, im2m}}];
        regFOV{fn} = [regFOV{fn} {{im1, im2}}];
        subplot(1, 2, 1)
        imshowpair(sum(im1m, 3), sum(im2m, 3))
        title('Masks')
        subplot(1, 2, 2)
        imshowpair(im1, im2)
        title('FOV ')
        pause(0.01)
    end
end

c = 1;
for fn = 1:length(group)
    FOVcell = regFOV{fn};
    for nn = 1:length(FOVcell)
        FOVcorr(fn, nn) = corr(FOVcell{nn}{1}(:), FOVcell{nn}{2}(:));
        subplot(2, 1, 1)
        imshowpair(FOVcell{nn}{1}, FOVcell{nn}{2})
        title(num2str(FOVcorr(fn, nn)))
        subplot(2, 1, 2)
        plot(FOVcorr)
        hold on
        plot(fn, FOVcorr(fn, nn), '.')
        xlabel('Group number')
        legend({'Pair 1', 'Pair 2', 'Pair 3'})
        ylabel('Correlation')
        yline(0.2)
        hold off
        pause
        c = c+1;
    end
end

matches = cell(length(group), 1);
for fn = 1:length(group)
    tROI = regROI{fn};
    for nn = 1:length(tROI)
        if FOVcorr(fn, nn) > 0.2
            m2 = tROI{nn}{2};
            m2(m2 < 0.2) = 0;
            [M] = matchNeurons(tROI{nn}{1}, m2, 0.5);
            matches{fn} = [matches{fn} {M}];
            subplot(2, 1, 1)
            
            imshowpair(sum(tROI{nn}{1}, 3), sum(tROI{nn}{2}, 3))
            title('All ROIs')
            subplot(2, 1, 2)
            
            imshowpair(sum(tROI{nn}{1}(:, :, M(:, 1)), 3), sum(tROI{nn}{2}(:, :, M(:, 2)), 3))
            title('Matched ROIs')
            pause
        end
    end
end

 save(strcat(datestr(now,"yyyy.mm.dd_HH.MM.SS"),'_registeredFOVs.mat'), 'Regc', 'regROI', 'regFOV', 'ROIs', 'rFOV', 'mMove', 'group', 'sttype', 'vDates', 'FOV', 'matches', 'FOVcorr', '-v7.3');