function [svid, MoveXt, MoveYt, offset, useorder, tCOM] = moveGanglia(filtgang, gangCOM, neighGang, gangMask, fluo, movement, t)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Generates a video of the ganglia and internodal strands moving without a
%background
% Inputs:
%     - filtgang: cell array of the masks that go into each ganglia, which
%     have been translated to their desired location in the video
%     - gangCOM: The COM of each ganglia
%     - neighgang: Cell array containing the indexes of the neighboring ganglia for each
%     individual ganglia
%     - gangMask: MxNxnGanglia matrix that contains the background mask for
%     each ganglia
%     - fluo: nNeuronxT matrix of the fluorescence traces to be used for
%     the neurons
%     - movement: structure from genMovement_fit
%     - t: number of frames
% Output:
%     - svid: Simulated video of the ganglia and internodal strand movement
%     - MoveXt: The base movement in the X dimension over the entire image
%     - MoveYt: The base movement in the Y dimension over the entire image
%     - offset: The calculated offset for the beginning of the movement
%     calculated based on the ganglia's center of mass
%     - useorder: the order in which the traces from fluo were randomly
%     assigned to the neurons in the video
%     - tCOM: nGangliaxframex2 matrix of the X and Y movement applied to
%     each of the different ganglia
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

MoveXminor = movement.MoveXminor;
MoveYminor = movement.MoveYminor;
MoveXmajor = movement.MoveXmajor; 
MoveYmajor = movement.MoveYmajor;
propangle = movement.propangle;
propangleortho = movement.propangleortho;

%Final movement trajectory
MoveXt = MoveXmajor+MoveXminor;
MoveYt = MoveYmajor+MoveYminor;

%Smooth the edges of the neurons so that they appear more natural
[d1 d2 d3] = size(cat(3, filtgang{:}));
[Xm Ym] = meshgrid(1:d1, 1:d2);

%Determine the maximum new x value obtained when the axis are rotated along the
%propagation angle
maxvals = Xm*cos(propangle)-Ym*sin(propangle);
minval = min(maxvals(:));
maxval = max(maxvals(:))-minval;

maxvalso = Xm*cos(propangleortho)-Ym*sin(propangleortho);
minvalo = min(maxvalso(:));
maxvalo = max(maxvalso(:))-minvalo;

fluo = fluo+1.15;

%Determine where ganglia connect to each other to draw the internodal
%strands and set which fluorescence traces are used
[loc1 loc2 loc3] = findConnectPoints(filtgang, neighGang);
[ntrace nframe] = size(fluo);
useorder = 1:d3;
count = 1;
tCOM = [];
blankvideo = single(zeros(d1, d2, length(t)));

gaussmu = 70;
gaussstd = 210;
gaussdiss = normpdf(maxvals, gaussmu, gaussstd);
gaussscale = 0.25;
gaussdiss = gaussdiss/max(gaussdiss(:));
gaussdiss = gaussdiss*gaussscale+1;
    
toff = 20;

%Apply the movement first to the ganglia
for fn = 1:size(gangCOM, 1)
    fn
    relGang = single(filtgang{fn}); 
    backgroundgang = single(gangMask(:, :, fn));
    [Xn Yn] = findCOM(relGang);
    rotx = nanmean(Xn)*cos(propangle)-nanmean(Yn)*sin(propangle);
    offset(fn) = double(ceil(toff*(1-((rotx-minval)/maxval))));
    if offset(fn) > 0
        moveXmajorg = padarray(MoveXmajor, [0 offset(fn)], 0, 'pre');
        moveYmajorg = padarray(MoveYmajor, [0 offset(fn)], 0, 'pre');
    else
        moveXmajorg = MoveXmajor;
        moveYmajorg = MoveYmajor;
    end
    gaussmod = gaussdiss(ceil(nanmean(Xn)), ceil(nanmean(Yn)));
    moveXmajorg = moveXmajorg*gaussmod;
    moveYmajorg = moveYmajorg*gaussmod;
    
    rotx = nanmean(Xn)*cos(propangleortho)-nanmean(Yn)*sin(propangleortho);
    offsetminor(fn) = double(ceil(toff*(1-((rotx-minvalo)/maxvalo))));
    if offsetminor(fn) > 0
        moveXminorg = padarray(MoveXminor, [0 offsetminor(fn)], 0, 'pre');
        moveYminorg = padarray(MoveYminor, [0 offsetminor(fn)], 0, 'pre');
    else
        moveXminorg = MoveXminor;
        moveYminorg = MoveYminor;
    end
    
    moveX = moveXmajorg(1:length(t))+moveXminorg(1:length(t));
    moveY = moveYmajorg(1:length(t))+moveYminorg(1:length(t));
    
    img = single(zeros(d1, d2, length(t)));
    count2 = count;
    tic
    for nn = 1:size(relGang, 3)
        linGang = reshape(relGang(:, :, nn), [d1*d2 1]);
        linFluo = linGang*fluo(useorder(count2), 1:length(t));
        reshGang = reshape(linFluo, size(img));
        img = img+reshGang;
        count2 = count2+1;
    end
    clear linFluo linGang reshGang
    toc
    tic
    blankGanglia = single(zeros(d1, d2, length(t)));
    parfor frame = 1:length(t)        
        indgang = img(:, :, frame);  
        tCOM(fn, frame, :) = [moveX(frame) moveY(frame)];
        medint = max(indgang(indgang ~=0));
        
        if all(tCOM(fn, frame, :) == 0)
            newpos = indgang;
            bgn = backgroundgang;
            newposind = find(newpos);
            bgind = find(bgn);
            newind = intersect(newposind, bgind);
            lowind = find(newpos(newind) < 0.1*medint);
            newpos(newind(lowind)) = 0.1*medint;
        else
            newpos = indgang;
            bgn = backgroundgang;
            newposind = find(newpos);
            bgind = find(bgn);
            newind = intersect(newposind, bgind);
            lowind = find(newpos(newind) < 0.1*medint);
            newpos(newind(lowind)) = 0.1*medint;
            newpos = imtranslate(newpos, squeeze(tCOM(fn, frame, :))');
        end

        blankGanglia(:, :, frame) = newpos;    
    end
    blankvideo = blankvideo+blankGanglia;
    toc
    count = count+size(relGang, 3);
end

blankimage = zeros(d1, d2);
svid = blankvideo;

%Randomly select the brightness and width of the internodal strands
for fn = 1:length(neighGang)
    neighbors = neighGang{fn};
    for nn = 1:length(neighbors)
        exsize(fn, nn) = randi(1)+1;
        brightness(fn, nn) = (0.45*rand)+0.3;
    end
end

%Insert the movement of the internodal strands
parfor frame = 1:length(t)
    (frame/length(t))*100
    tic
    for fn = 1:length(neighGang)
        neighbors = neighGang{fn};
        for nn = 1:length(neighbors)
            if neighbors(nn) > fn
                sloc1 = loc1{fn}(:, nn);
                sloc2 = loc2{fn}(:, nn);
                sloc3 = loc3{fn}(:, nn);
                tCOM1 = squeeze(tCOM(fn, frame, :));
                tCOM2 = squeeze(tCOM(neighbors(nn), frame, :));
                [index] = drawStrands(blankimage, sloc1, sloc2, tCOM1, tCOM2, sloc3);
                tframe = blankvideo(:, :, frame);
                vframe = blankimage;
                neurind = find(tframe(:));
                
                vframe(index) = quantile(tframe(neurind), brightness(fn, nn));
                vframe(index(ismember(index, neurind))) = vframe(index(ismember(index, neurind)))*0.3;
                
                SE = strel("disk", exsize(fn, nn));
                vframe = imdilate(vframe, SE);
                index = find(vframe);
                subind = randsample(1:length(index), ceil(length(index)*0.2));
                vframe(index(subind)) = 0;
                vframe = imgaussfilt(vframe, 3);
                
                
                
                svid(:, :, frame) = vframe+svid(:, :, frame);
            end
        end
    end
    toc
end
end