function [moveMask] = genMovMask(maskfolder, vidfolder, ogfolder, vidname)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create a video of the moving masks to demonstrate how the fluorescence is
% calculated
% Inputs:
%     - maskfolder: string of the folder containing the masks to be used for the COM
%     locations
%     - vidfolder: String of the folder containing the registered videos that the COM
%     points are draw on
%     - ogfolder: String of the folder containing the unregistered video to
%     be used for showing how the COM movement tracks the movement of the
%     video
%     - vidname: string Name of the video to use
% Outputs:
%     - moveMask: the video of the moving masks
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sigma = 10; 
nbackground = 40;
greenind = strfind(vidname, 'okada');
if isempty(greenind)
    greenind = strfind(vidname, '_');
end

% Find the data in the folders and load it
cd(maskfolder)
maskname = dir(strcat(vidname(1:greenind), '*masks.mat'));
if isempty(maskname)
    F = [];
    B = [];
    xTraj =[];
    yTraj =[];
else
    
    load(maskname(1).name)

    assocname = dir(strcat(vidname(1:greenind), '*associated.mat'));
    load(assocname(1).name, 'maxmove');
    maskmax = maxmove;
    
    cd(vidfolder)
    vidn = dir(strcat(vidname(1:greenind), '*.h5'));
    video = h5read(vidn(1).name, '/video3');
    assocname = dir(strcat(vidname(1:greenind), '*associated.mat'));
    load(assocname(1).name);
    vidmax = maxmove;
    
    savequad = {};

    %Go through each iterate of the registration and identify which patch each
    %pixel corresponds to. This is then used to construct the full movement
    %trajectory of a pixel
    for ind = 1:length(NCparamStruct.overlap)
        
        options_nr =NoRMCorreSetParms('d1',d1(ind+1),'d2',d2(ind+1),'bin_width',50, ...
            'grid_size',NCparamStruct.gridsize(ind, :),'mot_uf', 1,'correct_bidir',false, ...
            'overlap_pre',NCparamStruct.overlap(ind),'overlap_post', NCparamStruct.overlap(ind)/2,'max_shift',....
            NCparamStruct.maxshift(ind), 'upd_template', false, 'shifts_method', 'cubic');
        
        
        relquad = reconstructMovement(options_nr, d1(ind+1), d2(ind+1), 1);
        if ind ~= length(NCparamStruct.overlap)
            sumove = sum(vidmax((ind+1):3, :), 1);
        else
            sumove = [0 0 0 0];
        end
        tind = sum(vidmax(1:ind, :), 1);
        trackval = zeros(d1(1), d2(1), 2);
        trackval(tind(3)+1:end-tind(4), tind(1)+1:end-tind(2), :) = relquad;
        savequad{ind} = trackval;
        relquad = relquad(sumove(3)+1:end-sumove(4), sumove(1)+1:end-sumove(2), :);
        
    end
    
    %Convert the masks and video to the correct field of view
    [ssmask, ssvid, nmaxmove] = convergeMaxMove(masks, video, maskmax, vidmax);
    
    if size(ssmask, 1) ~= size(ssvid, 1)
        [ssmask, ssvid, nmaxmove] = convergeMaxMove(permute(masks, [2 1 3]), video, maskmax, vidmax);
    end
    
    
    [v1 v2 v3] = size(video);
    [X Y] = findCOM(double(ssmask));
    neuronnum = 1:length(X);
    X = ceil(X);
    Y = ceil(Y);
    bind = [];
    for fn = 1:length(X)
        if isnan(X(fn))
            bind = [bind fn];
        end
    end
    
    useind = 1:length(X);
    useind(bind) = [];
    for fn = 1:length(useind)
        nind = find(ssmask(:, :, useind(fn)));
        [r c] = ind2sub(size(ssmask(:, :, fn)), nind);
        rloc{fn} = r+nmaxmove(3);
        cloc{fn} = c+nmaxmove(1);
    end
    
    X(bind) = [];
    Y(bind) = [];
    neuronnum(bind) = [];
    
    %Generate the COM trajectories by identifying the relevant patch for each
    %COM for each iteration of the registration
    xTraj = repmat(X+nmaxmove(1), 1, length(shifts2{2}));
    yTraj = repmat(Y+nmaxmove(3), 1, length(shifts2{2}));
    
    
    for ind = fliplr(1:length(NCparamStruct.overlap))
        relshift = peakshift{ind+1};
        shifts_nr = cat(ndims(relshift(1).shifts)+1,relshift(:).shifts);
        relquad = savequad{ind};
        dim1 = d1(ind+1);
        dim2 = d2(ind+1);
        
        parfor frame = 1:size(xTraj, 2)
            shifts_up = imresize(relshift(frame).shifts,[dim1 dim2]);
            shiftsx = relquad(:, :, 2);
            shiftsx(find(shiftsx)) = shifts_up(:, :, :, 2);
            shiftsy = relquad(:, :, 1);
            shiftsy(find(shiftsy)) = shifts_up(:, :, :, 1);
            Yc = floor(yTraj(:, frame));
            Yc(Yc < 1) = 1;
            Xc = floor(xTraj(:, frame));
            Xc(Xc < 1) = 1;
            trackind = sub2ind(size(shiftsy), Yc, Xc);
            xTraj(:, frame) = xTraj(:, frame)-shiftsx(trackind);
            yTraj(:, frame) = yTraj(:, frame)-shiftsy(trackind);
        end
        
    end
    xTraj = xTraj-repmat(xTraj(:, 1)-(X+nmaxmove(1)), 1, size(xTraj, 2));
    yTraj = yTraj-repmat(yTraj(:, 1)-(Y+nmaxmove(3)), 1, size(yTraj, 2));

    %Draw the video
    cd(ogfolder)
    vname = dir(strcat(vidname(1:(greenind-2)), '*.h5'));
    vname(1).name
    vidname

    video = h5read(vname(1).name, '/video');

    [d1 d2 d3] = size(video);
    moveMask = single(zeros(size(video)));
    clear video;
    
    parfor frame = 1:d3
        newframe = moveMask(:, :, frame);
        for neuron = 1:length(neuronnum)
            rpos = rloc{neuron};
            cpos = cloc{neuron};
            
            rpos = rpos+round(yTraj(neuron, frame))-(Y(neuron)+nmaxmove(3));
            cpos = cpos+round(xTraj(neuron, frame))-(X(neuron)+nmaxmove(1));
            if all(rpos < d1) && all(rpos >= 1) && all(cpos >= 1) && all(cpos < d2)
                newind = sub2ind(size(newframe), rpos, cpos);
            else
                remind = find(rpos >= d1 | rpos < 1 | cpos < 1 | cpos >= d2);
                rpos(remind) = [];
                cpos(remind) = [];
                newind = sub2ind(size(newframe), rpos, cpos);
            end
            newframe(newind) = neuronnum(neuron);
        end
       moveMask(:, :, frame) = newframe;
    end
end