function focusCorrelation(regfolder,segfolder,maskfolder,focusfolder)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Generates the focus mat file to use for training the logistic regression
% Inputs:
%     - regfolder: directory where the registered data is stored
%     - segfolder: directory with the segementation .mat files
%     - maskfolder: directory with the mask .mat files (same as segfolder?)
%     - focusfolder: directory with annotated consensus focus data
% Output: (none)
%     - 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

cd(regfolder)
vName = dir('*registered.h5');
count = 1;
for fn = 1:length(vName)
    vidname = vName(fn).name;
    greenind = strfind(vidname, 'Green')+4;
    cd(maskfolder)
    B = dir(strcat(vidname(1:greenind), '*masks.mat'));
    cd(focusfolder)
    fT = dir(strcat(vidname(1:greenind), '*tion.mat'));
    cd(regfolder)
    if ~isempty(fT)
        cd(regfolder)

    vidname = vName(fn).name;
    greenind = strfind(vidname, 'Green')+4;
    A = dir(strcat(vidname(1:greenind), '*.h5'));
    regvid = h5read(A(1).name, '/video3');
    assocdat = dir(strcat(vidname(1:greenind), '*associated.mat'));
    load(assocdat.name, 'maxmove');
    regmax = maxmove; 

    matname = strcat(maskfolder, B(1).name);
    [gangmask, idx] = makeGanglia_kmeans(matname);
    pmask = [];
    for mask = 1:length(gangmask)
        ganglia = gangmask{mask};
        pmask(:, :, mask) = sum(ganglia, 3);
    end

    cd(segfolder)

    assocdat = dir(strcat(vidname(1:greenind), '*associated.mat'));
    load(assocdat.name, 'maxmove');
    segmax = maxmove; 
    [d1 d2 d3] = size(pmask);
    pmaskc = pmask(maxmove(4, 3)+1:end-maxmove(4, 4), maxmove(4, 1)+1:end-maxmove(4, 2), :);
    regvid = regvid(maxmove(4, 3)+1:end-maxmove(4, 4), maxmove(4, 1)+1:end-maxmove(4, 2), :);

    remind = [];
    for mask = 1:size(pmaskc, 3)
        relmask = pmaskc(:, :, mask);
        if sum(relmask(:)) == 0
            remind(mask) = 1;
        else
            remind(mask) = 0;
        end
    end
    pmaskc(:, :, remind==1) = [];
    template = mean(regvid(:, :, 1:300), 3);

    ntemplate = (template-mean(template(:)))/std(template(:));
    matchdat = dir(strcat(vidname(1:greenind), '*matchFocusSegmentation.mat'));
    load(matchdat(1).name)
    fullquad = zeros(d1, d2);

    fullquad(segmentationcrop(3):segmentationcrop(4), segmentationcrop(1):segmentationcrop(2)) = quadrants;
    fullquad = fullquad(maxmove(4, 3)+1:end-maxmove(4, 4), maxmove(4, 1)+1:end-maxmove(4, 2));

    clear regvid1
    vcorr = [];
    vcorrquad = [];
    quadloc = [];
    vssim = [];
    vsnr = [];
    vdiff = [];
    vgrad = [];
    vIOU = [];
    vMAD = [];
    vkurt = [];
    for mask = 1:4
        maskind = find(fullquad == mask);
        quadloc(mask) = mask;
        [r c] = ind2sub(size(regvid), maskind);
        smalltemp = ntemplate(min(r):max(r), min(c):max(c));
        smallvid = regvid(min(r):max(r), min(c):max(c), :);
        for frame = 1:size(smallvid, 3)
            vframe = smallvid(:, :, frame);
            vframe = normalize(double(vframe));
            vframe(isnan(vframe)) = min(vframe(:));
            vcorrquad(mask, frame) = corr2(vframe, smalltemp);
        end
    end
        
    se = strel('disk', 4);   
    projframes = 5;
    
    for mask = 1:size(pmaskc, 3)
        relmask = pmaskc(:, :, mask);
        relmaskdil = imdilate(relmask, se);
        maskind = find(relmaskdil);
        quadloc(mask) = mode(fullquad(maskind));
        if ~all(relmask(:) == 0)
            maskind = find(relmaskdil);
            [row col] = ind2sub(size(relmask), maskind);
            r = [min(row) max(row)];
            c = [min(col) max(col)];
            
            smalltemp = ntemplate(r(1):r(2), c(1):c(2));
            smallvid = regvid(r(1):r(2), c(1):c(2), :);
            relmaskdil = relmaskdil(r(1):r(2), c(1):c(2), :);
            relmask = relmask(r(1):r(2), c(1):c(2), :);
            relind = find(relmask);
            nrelind = find(relmask == 0);
            relinddil = find(relmaskdil);
            vcorr_temp = zeros(1, size(smallvid, 3));
            vssim_temp = zeros(1, size(smallvid, 3));
            vsnr_temp = zeros(1, size(smallvid, 3));
            vdiff_temp = zeros(1, size(smallvid, 3));
            vgrad_temp = zeros(1, size(smallvid, 3));
            vIOU_temp = zeros(1, size(smallvid, 3));
            level = graythresh(smalltemp);
            BWt = imbinarize(smalltemp,level);
            tempind = find(BWt);
            parfor frame = 1:size(smallvid, 3)
                vframe = double(smallvid(:, :, frame));
                vframe(vframe < median(vframe(:))-4*mad(vframe(:), 1)) = median(vframe(:));
                vsnr_temp(frame) = mean(vframe(relind))/std(vframe(nrelind));
                vframe = (vframe-mean(vframe(:)))/std(vframe(:));
                vframe(isnan(vframe)) = min(vframe(:));
                diffframe = vframe-smalltemp; 
                vdiff_temp(frame) = mean(diffframe(relind))-mean(diffframe(relinddil));
                vcorr_temp(frame) = corr(vframe(:), smalltemp(:));
                vssim_temp(frame) = ssim(vframe, smalltemp);
                sframe = max([1 frame-projframes]);
                avgframe = mean(double(smallvid(:, :, sframe:frame)), 3);
                avgframe = (avgframe-mean(avgframe(:)))/std(avgframe(:));
                [gx,gy] = gradient(avgframe); 
                vgrad_temp(frame) = norm(sqrt(gx.^2+gy.^2),'fro');
                level = graythresh(avgframe);
                BW = imbinarize(avgframe,level);
                BWind = find(BW);
                vIOU_temp(frame) = sum(ismember(BWind, tempind))/(length(BWind)+length(tempind)-sum(ismember(BWind, tempind)));
            end
            vgrad_temp(1:projframes) = vgrad_temp(projframes+1);
            vcorr(mask, :) = vcorr_temp;
            vssim(mask, :) = vssim_temp;
            vsnr(mask, :) = vsnr_temp;
            vdiff(mask, :) = vdiff_temp;
            vgrad(mask, :) = vgrad_temp;
            vIOU(mask, :) = vIOU_temp;
            [sd1 sd2 sd3] = size(smallvid);
            vMAD(mask, :) = std(double(reshape(smallvid, [sd1*sd2 sd3])));
            vkurt(mask, :) = kurtosis(double(reshape(smallvid, [sd1*sd2 sd3])));
        end
    end
    
    cd(focusfolder)
    fT = dir(strcat(vidname(1:greenind), '*tion.mat'));
    load(fT.name)
    order = [2 1 3 4];
    for mask = 1:size(vcorr, 1)
        if ~isnan(quadloc(mask)) && quadloc(mask) ~= 0
            subplot(2, 2, order(quadloc(mask)))
            plot(normalize(log10(vcorr(mask, :)), 'center', 'median'), 'k')
            hold on
            plot(normalize(real(log10(vssim(mask, :))), 'center', 'median'), 'r')
        end
    end
    for period = 1:size(focusTable, 1)
        for q = 1:4
            if focusTable(period, q+2) == 1
                subplot(2, 2, order(q))
                plot([focusTable(period, 1) focusTable(period, 2)], [1 1], 'k', 'LineWidth', 2) 
                xline(focusTable(period, 1))
                xline(focusTable(period, 2))
            end
        end
    end
    hold off
    for q = 1:4
        subplot(2, 2, q)
        ylim([0 1])
        xlim([0 2400])
    end
    focus_layers(count).correlation = single(vcorr);
    focus_layers(count).quadrant = quadloc;
    focus_layers(count).focusTable = focusTable;
    focus_layers(count).folder = regfolder;
    focus_layers(count).name = vidname;
    focus_layers(count).plexus = uint8(pmask);
    focus_layers(count).quadcorr = single(vcorrquad);
    focus_layers(count).ssim = single(vssim);
    focus_layers(count).snr = single(vsnr);
    focus_layers(count).diff = single(vdiff);
    focus_layers(count).grad = single(vgrad);
    focus_layers(count).IOU = single(vIOU);
    focus_layers(count).MAD = single(vMAD);
    a = whos('focus_layers'); 
    a.bytes*(1e-9);
    clf
    count = count+1;
    else
    end
end