function [newshift] = DijkstraPeakSearch(peakshift, nomshifts, transKey, distthresh, fillthresh, buffer, maxshift, doesplot)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Identifies the shortest path through a matrix of X and Y nodes to use
% for the registration
% Inputs:
%     - peakshift: struct containing the 2D coordinates of the nodes that
%     are potential registration coordinates
%     - nomshifts: Struct containing the registration coordinates generated
%     by NoRMCorre without memory (using the peak correlation value)
%     - transKey: Cell array that contains the conversion from pixel
%     location in the correlation image to the shift size to be used by
%     NoRMCorre for each patch.  
%     - distthresh: Cutoff threshold for movement distances
%     - fillthresh: Cutoff threshold for filling in missing nodes
%     - buffer: 40 frames, used as temporal window
%     - maxshift: Maximum pixel shift
%     - doesplot: Boolean to plot the data from figure 1
% Output:
%     - newshift: Struct containing the new registration coordinates
%     calculated using Dijkstra's algorithm. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

shifts_nr = cat(ndims(peakshift(1).shifts)+1,peakshift(:).shifts);
peaks_nr = cat(ndims(peakshift(1).PeakVal)+1,peakshift(:).PeakVal);
threshold = cat(ndims(peakshift(1).threshold)+1,peakshift(:).threshold);
[d1, d2, ~, ~, d5] = size(shifts_nr); 

% Threshold for total movement from 0 needed to perform path finding
jthresh = 20; 

% Threshold for the standard deviation of the base norm corre path to recalculate
stdthresh = 20;

% A buffer for frames identified for path finding
time = 41; 

% Weighting of the heuristic
alpha = 0.5; 

if doesplot == 1
    f = figure;
    f.Position(1) = f.Position(1) - f.Position(3)/2;
    f.Position(3) = f.Position(3)*2;
    tg = uitabgroup;
    panelcount = 1;
end

for fn = 1:d1
    for nn = 1:d2
        Nc = transKey{fn, nn, 1, 2};
        Nccorr = find(Nc == 0);
        Nr = transKey{fn, nn, 1, 1};
        Nrcorr = find(Nr == 0);
        xval = squeeze(shifts_nr(fn, nn, :, 2, :));
        yval = squeeze(shifts_nr(fn, nn, :, 1, :));
        
        fx = xval(:, 1);
        fxNaN = find(isnan(fx), 1, 'first')-1;
        if isempty(fxNaN)
            fxNaN = size(xval, 1);
        end
        fy = yval(:, 1);
        fyNaN = find(isnan(fy), 1, 'first')-1;
        if isempty(fyNaN)
            fyNaN = size(yval, 1);
        end
        Ncdiff = [0 diff(Nc)];
        Nrdiff = [0 diff(Nr)];
        intind = find(Nrdiff(fy(1:fyNaN)) > 0 & Ncdiff(fx(1:fxNaN)) > 0, 1, 'first');
        reorder = 1:size(xval, 1);
        reorder(intind) = [];
        reorder = [intind reorder];
        yval = yval(reorder, :);
        xval = xval(reorder, :);
        
        
        
        
        threshval = squeeze(threshold(fn, nn, :));
        sposx = find(Nc == round(nomshifts(1).shifts(fn, nn, :, 2)));
        sposy = find(Nr == round(nomshifts(1).shifts(fn, nn, :, 1)));
        
        xvalstd = movstd(diff(xval(1, :)), time);
        yvalstd = movstd(diff(yval(1, :)), time);
        maxstd = max(cat(1, xvalstd, yvalstd));
        stdthresh = max([median(maxstd)+3*mad(maxstd, 1) stdthresh]);
        ind = find(maxstd' > stdthresh);

        ftrajmov = [sqrt(diff(xval(1, :)).^2+diff(yval(1, :)).^2) 0];
        shifttrajmov = max(cat(1, abs(xval(1, :)-length(Nc)/2), abs(yval(1, :)-length(Nr)/2)));
        jind = find(shifttrajmov' > distthresh | ftrajmov' > jthresh);
        ind = unique(cat(1, ind, jind));
        blank = ones(size(threshval));
        blank(ind) = 0;
        starts = strfind([1 blank'], [1 0]);
        ends = strfind([blank' 1], [0 1]);
        for chunk = 1:length(starts)
            nstart = max([starts(chunk)-buffer 1]);
            nend = min([ends(chunk)+buffer length(blank)]);              
            blank(nstart:starts(chunk)) = 0;
            blank(ends(chunk):nend) = 0;
        end
        blank(end) = 1;
        bind = find(blank);
        
        
        xval(2:end, bind) = NaN;
        yval(2:end, bind) = NaN;
        
        pastpos = [yval(1, 1) xval(1, 1)];
        newshift(1).shifts(fn, nn, 1, :) = [Nr(pastpos(1))/2 Nc(pastpos(2))/2];
        newshift(1).shifts_up(fn, nn, 1, :) = [Nr(pastpos(1))/2 Nc(pastpos(2))/2];
        xpath = xval;
        ypath = yval;
          
           
        [xpath ypath estx esty] = fillNodes(xpath', ypath', fillthresh);
        remind = find(xpath > length(Nc) | ypath > length(Nr) | xpath < 1 | ypath < 1);
        xpath(remind) = NaN;
        ypath(remind) = NaN;

        valind = find(~isnan(xpath));
        movement = sqrt(Nc(round(xpath(valind))).^2+Nr(round(ypath(valind))).^2);
        exmov = find(movement > maxshift*2);
        xpath(valind(exmov)) = NaN;
        ypath(valind(exmov)) = NaN;
        xblank = xpath;
        yblank = ypath;
        for frame = 1:d5
            ind = find(~isnan(xblank(frame, :)));
            xpath(frame, :) = NaN;
            xpath(frame, 1:length(ind)) = xblank(frame, ind);
            ypath(frame, :) = NaN;
            ypath(frame, 1:length(ind)) = yblank(frame, ind);
        end
        
         xest = xval(1, :);
           yest = yval(1, :);
           starts = strfind([1 blank'], [1 0]);
           ends = strfind([blank' 1], [0 1]);
           for chunk = 1:length(starts)
               s = starts(chunk);
               e = ends(chunk);
               
               yest(1+(s:e)) = yest(s+1)+cumsum(esty(s:e));
               xest(1+(s:e)) = xest(s+1)+cumsum(estx(s:e));
               sly = (yval(1, e)-yest(e))/length(s:e);
               slx = (xval(1, e)-xest(e))/length(s:e);
               yest(1+(s:e)) = yest(s:e)+sly*(1:length(s:e));
               xest(1+(s:e)) = xest(s:e)+slx*(1:length(s:e));
           end
        path = Astar2D_estimate(xpath, ypath, xest, yest, alpha);
        xplot = xpath;
        yplot = ypath; 
        logx = round(xplot(path));
        logy = round(yplot(path));
        logy(logy < 1) = 1;
        logy(logy > length(Nr)) = length(Nr);
        logx(logx < 1) = 1;
        logx(logx > length(Nc)) = length(Nc);
        
        
        for tt = 2:d5
            newshift(tt).shifts(fn, nn, 1, :) = [Nr(logy(tt))/2 Nc(logx(tt))/2];
            newshift(tt).shifts_up(fn, nn, 1, :) = [Nr(logy(tt))/2 Nc(logx(tt))/2];
        end
        if doesplot == 1
            thistab = uitab(tg);
            %Plot the data for the view used in figure 1
            axes('Parent',thistab);
            subplot(2, 1, 1)
            plot(xpath-Nccorr, 'b.')
            hold on
            plot(xval'-Nccorr, 'k.');
            ylabel('Column shift')

            plot(xval(intind, :)-Nccorr, 'c', 'LineWidth', 1.2)
            plot(xest-Nccorr, 'm', 'LineWidth', 1.2)
            plot(xplot(path)-Nccorr, 'r', 'LineWidth', 1.2)
            
            hold off
            subplot(2, 1, 2)
            plot(ypath-Nrcorr, 'b.')
            hold on
            plot(yval'-Nrcorr, 'k.');
            xlabel('Frame')
            ylabel('Row shift')
            plot(yval(intind, :)-Nrcorr, 'c', 'LineWidth', 1.2)
            plot(yest-Nrcorr, 'm', 'LineWidth', 1.2)
            plot(yplot(path)-Nrcorr, 'r', 'LineWidth', 1.2)
            hold off
            thistab.Title = num2str(panelcount);
            panelcount = panelcount +1;
        end
    end
end
end