function [mosaic,pth] = compute_mosaic(movie, dominant_pan_dir, ALLOW_Y)
% Compute mosaic from movie using Dijkstra algorithm with video pyramid.
%
% INPUT
%   movie - a string (e.g., 'orig/*.jpg' or 'movie.avi')
%      Note: an avi movie should be readable by Matlab. 
%            For example, in UNIX it should most probably be
%            uncompressed.
%   dominant_pan_dir =  1 if camera panning is to the right; 
%                    = -1 if camera panning is to the left
%    /if wrong pan direction is chosen, mosaic will be degenerate (small)/
%   ALLOW_Y (default = false) - whether to search for vertical shifts 
%          (in y dimension). 
%          Caution: in the current (very unoptimal) implementation setting 
%             ALLOW_Y=true may enlarge running time considerably.
% OUTPUT
%   mosaic (YxNx3 uint8) - mosaic image
%   pth (struct) - pth.x, pth.t, pth.y_shift (1xN) are coordinates
%                  of the found optimal path. 
%
% USAGE
%   [mosaic,pth] = compute_mosaic(movie, dominant_pan_dir, ALLOW_Y=false)
% EXAMPLES
%   [mosaic,pth] = compute_mosaic('movie/*.jpg', 1);
%   [mosaic,pth] = compute_mosaic('movie.avi', -1);
%   [mosaic,pth] = compute_mosaic('jitter_movie.avi', -1, true);

% This function basically sets default parameters and calls mosaic_main()

if ~exist('ALLOW_Y','var'); ALLOW_Y = false; end

% =========== Set up default parameters ============== %

% Reference points (start and end points of the path)
%   ref_points = [x0 xend; t0 tend]
if dominant_pan_dir==1 % to the right
    % from the first to the last frame
    ref_points=[-inf inf;-inf +inf]; 
else % to the left
    % from the last to the first frame
    ref_points=[-inf inf;+inf -inf]; 
end 

% 2D mask (in xt dimensions): restrict the path
mask=[]; % all points

% Parameters:
%     /see pyramid_level_from_scratch() for more detailed 
%        description of parameter structs/
% Path and cost function
comp_params=struct(...
    'trange',[-5 -3:3 5], ... % frame shifts to explore
    'path_norm_type',10,... % combine arc weights using L_p norm
    'xwnd',10,... % size of the window in x dimension
    'ysubwnd',10); % size of subwindows in y 
% How much freedom to give on a finer pyramid level
exp_params=struct(...
    'xvel_wnd',7,...
    'yvel_wnd',0,...
    'dilation_diam',[7 7 1]);
% Dynamics of the video
dyn_params=struct(...
    'xvel',7*[-1 1],... % velocity range in x dimension (dx/dt)
    'yvel',[0 0],... % velocity range in y dimension (dy/dt)
    'max_y_shift',0); % maximal shift in y accounted for

if ALLOW_Y
    exp_params.yvel_wnd = 3;
    exp_params.dilation_diam = [7 7 7];
    dyn_params.yvel = [-3 3];
    dyn_params.max_y_shift = 10;
end

% =========== Run main function ============== %

pyr = mosaic_main(movie, ref_points,mask,...
    comp_params,exp_params,dyn_params);

% =========== Return results ============== %

mosaic = pyr{1}.mosaic;
pth = pyr{1}.pth;


% ======================== %
%    actually_get_frame    %
% ======================== %

function img=actually_get_frame(mov, frame)
% Obtain frames of the movie
% INPUT
%  mov (1x1 struct) - movie, see load_movie.m
%  frame (1xF) - frame numbers
% 
% OUTPUT
%  img (HxWxCxF uint8) - movie frames. C=1 for grayscale, 3 for color.

%disp(['loading frames ' int2str(frame) ]);
N=length(frame);
img = uint8(0);
img(mov.framesize(1), mov.framesize(2), ...
    mov.framesize(3), length(frame)) = 0;

switch mov.input_type
    case 'avi'
        if isfield(mov,'framerange');
            frame=mov.framerange(frame);
        end
        pfx = ''; if isfield(mov,'path_prefix'); pfx = mov.path_prefix; end
        tmp=aviread([pfx mov.files], frame);
        for i=1:N
            img(:,:,:,i)=tmp(i).cdata;
        end
    case 'images'
        pfx = ''; if isfield(mov,'path_prefix'); pfx = mov.path_prefix; end
        if isfield(mov,'read_func');
            rf=mov.read_func;
            for i=1:N; 
                tmp=feval(rf,[pfx mov.files{frame(i)}]); 
                img(:,:,:,i)=tmp;
            end
        else
            for i=1:N; 
                tmp=imread([pfx mov.files{frame(i)}]); 
                img(:,:,:,i)=tmp;
            end
        end
    case 'array'
        img(:,:,:,:) = mov.video_data(:,:,:,frame);
    otherwise
        error('Unknown movie type %s', mov.input_type);
end


% ====================== %
%    best_path_forest    %
% ====================== %

function forest = best_path_forest(start_nodes,stop_nodes, n_nodes, path_norm_type,graph_func,varargin)
% forest = best_path_forest(start_nodes,stop_nodes, n_nodes, path_norm_type,graph_func,GRAPH_FUNC_PARAMS);
% Computes forest of best paths, starting from any of 'start_nodes'.
% The forest is built until any of the 'stop_nodes' is reached, or till the end.
% 
% INPUT:
% start_nodes, stop_nodes (1xK,1xL) - lists of nodes. Forest is built for finding
%                the best path strarting at any of the start_nodes and ending at any of the stop_nodes.
% n_nodes (1x1) - total number of nodes (used for displaying progress)
% path_norm_type (1x1) = 1 (sum) or Inf (max) - operation to compute path cost from weights of arcs in the path
% graph_func (e.g., @graph_func_movie) - handler of the graph function
% GRAPH_FUNC_PARAMS - parameters required by the graph_func
%
% OUTPUT: 
% (see forest_struct.m for description of the fields)
% forest (1x1) - structure with fields:
%   alias (1xN)
%   parent (1xN)
%   path_cost (1xN), arc_weight (1xN)

% Initialization for Dijkstra algorithm
final=zeros(0,6);
pending=zeros(length(start_nodes),6);
pending(:,2)=start_nodes(:); pending(:,3)=-1;

% Run Dijkstra algorithm 
START_CLOCK=clock;
% fprintf('\t [%s] Started Dijkstra algorithm\n',datestr(START_CLOCK,21));
[final,pending,stop_reason,n_arcs_visited] = dijkstra_online_core(final,pending, inf,inf,stop_nodes, START_CLOCK,n_nodes, path_norm_type, graph_func,varargin{:});
%  report
% fprintf('\t [%s]',datestr(clock,21));
fprintf('\t Dijkstra algorithm stopped after %s, reason: %s\n\t\t best path computed for %d (%.1f%%) nodes, visited %d arcs per node\n',...
    seconds2str(etime(clock,START_CLOCK)), stop_reason,size(final,1),100*size(final,1)/n_nodes,round(n_arcs_visited/size(final,1)));

% Create forest structure
forest=forest_struct(final);


% ==================== %
%    decimate_movie    %
% ==================== %

function mov=decimate_movie(mov,n)
% mov=decimate_movie(mov,n)
% Take every n's frame from movie.
%
% mov - movie structure (see load_movie())

switch mov.input_type
    case 'avi'
        if isfield(mov,'framerange');
            mov.framerange=mov.framerange(1:n:end);
        else
            mov.framerange=1:n:mov.nframes;
        end
        mov.nframes=length(mov.framerange);
    case 'images'
        mov.files=mov.files(1:n:end);
        mov.nframes=length(mov.files);
    case 'array'
        mov.video_data=mov.video_data(:,:,:,1:n:end);
        mov.nframes=size(mov.video_data,4);
end

if isfield(mov,'xtslice');
    mov.xtslice=mov.xtslice(:,1:n:end,:);
end

% ========================== %
%    dijkstra_online_core    %
% ========================== %

function [final,pending,stop_reason,n_arcs_visited] = dijkstra_online_core(final,pending, STOP_ITER,STOP_TIME,STOP_NODES, ...
    START_CLOCK,N_NODES, path_norm_type, graph_func,varargin)
% Dijkstra algorithm with on-line computation of arcs and weights
%   and a few stopping criterion options.
%
% Work structures:
%   final (Nx6) [dist node parent arc_weight path_step path_length] - finished part of tree, 
%   pending (Mx6) [dist node parent arc_weight path_step path_length] - other nodes with non-infinite estimation of distance
%                  parent->node is an arc in the tree with weight arc_weight;
%                  dist = distance from the root (-1) to the node, according to the path_norm_type
%                         for 'final' the distance is globally minimal; for 'pending' - estimated
%
% Stopping criteria:
%   STOP_ITER, STOP_TIME, STOP_NODES
%   Stop: 1) after certain number of iterations; 
%      OR 2) after given time; 
%      OR 3) upon reaching one of the given nodes
%
% Progress report:
%   START_CLOCK (from 'clock' function), N_NODES - total number of nodes (to calculate percentage)
% 
% Path computation:
%   path_norm_type (integer or inf), graph_func and its parameters (varargin)

THIS_START_TIME=clock;

last_mess_len=0; 
min_tm_interval=5;
TM_last_report=THIS_START_TIME;

STOP_NODES=sort(STOP_NODES); % for faster search afterwards

% MAIN LOOP (Dijkstra)
iter=0;
TM_cur=THIS_START_TIME;
n_arcs_visited=0;
while (iter<STOP_ITER) && (etime(TM_cur,THIS_START_TIME) < STOP_TIME)
    iter=iter+1;
    
    % Find the min-dist pending node, add it to final and remove from pending
    if isempty(pending); break; end % exhausted the connected component
    [dummy,ind]=min(pending(:,1));
    new=pending(ind,:);
    % remove from pending
    pending(ind,:)=[];
    % add to final
    final(end+1,:)=new;

    % (a faster version of "if any(STOP_NODES==new(2))")
    if ismembc(new(2),STOP_NODES); break; end % found a target node - EXIT
    
    % Find all arcs adjacent to the new node
    [adj,w,steps]=feval(graph_func, new(2),final(:,2), varargin{:});
    adj=adj(:); w=w(:); steps=steps(:); 
    n_arcs_visited=n_arcs_visited+length(adj);

    % Analyze adjacent arcs
    % (faster version, instead of ismember())
    adj_in_fin=ismembc(adj,sort(final(:,2))); 
    adj=adj(~adj_in_fin); w=w(~adj_in_fin); steps=steps(~adj_in_fin); 
    % Find nodes that are already present in the list of pending nodes
    % (faster version, instead of 
    %   [adj_in_pend_1,adj_in_pend_pos_1]=ismember(adj,pending(:,2)); )
    [pnd_srt,pnd_srt_idx]=sort(pending(:,2));
    adj_in_pend_pos=ismembc2(adj,pnd_srt); 
    adj_in_pend=adj_in_pend_pos>0;
    adj_in_pend_pos(adj_in_pend)=pnd_srt_idx(adj_in_pend_pos(adj_in_pend));
    %   these we need to add:
    adj_to_add=adj(~adj_in_pend); w_to_add=w(~adj_in_pend); steps_to_add=steps(~adj_in_pend); 
    %   these we need to update:
    adj_to_upd=adj(adj_in_pend); w_to_upd=w(adj_in_pend); steps_to_upd=steps(adj_in_pend); 
    adj_to_upd_i=adj_in_pend_pos(adj_in_pend);
    
    % Compute new distances
    if isinf(path_norm_type) % Linf norm
        new_dist_to_upd=max(new(1),w_to_upd);
        new_dist_to_add=max(new(1),w_to_add);
    else % Lk norm
        new_dist_to_upd=new(1)+w_to_upd.^path_norm_type;
        new_dist_to_add=new(1)+w_to_add.^path_norm_type;
    end
    
    % Update
    %   compute and compare distances
    if ~isempty(adj_to_upd);
        old_dist=pending(adj_to_upd_i,1);
        better= new_dist_to_upd<old_dist;
        %   update distance and predecessor when needed
        pending(adj_to_upd_i(better),1)=new_dist_to_upd(better); % distance in path
        pending(adj_to_upd_i(better),3)=new(2); % parent
        pending(adj_to_upd_i(better),4)=w_to_upd(better); % arc weight
        pending(adj_to_upd_i(better),5)=steps_to_upd(better); % step in the path (e.g., in pixels in x)
        pending(adj_to_upd_i(better),6)=new(6)+steps_to_upd(better); % path length
    end
    
    % Add
    if ~isempty(adj_to_add);
        add = [new_dist_to_add adj_to_add repmat(new(2),size(adj_to_add)) w_to_add steps_to_add new(6)+steps_to_add];
        % add to pending
        pending=[pending;add];
    end
    
    % Report progress and times
    TM_cur=clock;
    if etime(TM_cur,TM_last_report)>=min_tm_interval
        fprintf(repmat('\b',1,last_mess_len));
        last_final=size(final,1);
        last_pending=size(pending,1);
        TM_iter=etime(TM_cur,START_CLOCK);
        TM_iter_str=seconds2str(TM_iter);
        TM_estim=round((N_NODES/last_final-1)*TM_iter);
        TM_estim_str=seconds2str(TM_estim);
        last_mess_len=fprintf('    [Dijkstra: %s] %d (%.1f%%) final, %d (%.1f%%) queued;   %s to 100%%',...
            TM_iter_str,last_final,100*last_final/N_NODES,...
            last_pending,100*last_pending/N_NODES,...
            TM_estim_str);
        TM_last_report=TM_cur;
    end

end % main loop

% Clean the report 
fprintf(repmat('\b',1,last_mess_len));

% Return the reason for stopping
if ismembc(final(end,2),STOP_NODES); % found a target node - EXIT
    stop_reason='OK (found target node)';
elseif isempty(pending); % the forest cannot grow further
    stop_reason='exhausted the graph';
elseif iter>=STOP_ITER; % exceeded number of iterations
    stop_reason='max iterations';
elseif etime(TM_cur,THIS_START_TIME) >= STOP_TIME; % exceeded maximum time
    stop_reason='max time';
else
    stop_reason='unknown';
end


% ======================== %
%    dynamics_from_path    %
% ======================== %

function [xvel,yvel,yshift]=dynamics_from_path(pth)
% [xvel,yvel,yshift]=dynamics_from_path(pth)
% Derive parameters of the video dynamics given mosaic path. 
% 
% INPUT
%   pth - structure with fields:
%     x,t,y_shift (1xN) - coordinates of the path
%     xstep (1x1) - shift in mosaic (in pixels) corresponding to each arc 
%                   in the path
%
% OUTPUT
%   xvel,yvel (1x2) [min max] - ranges of velocities
%   yshift (1x2) [min max] - range of values of y_shift

if ~isfield(pth,'xstep'); pth.xstep=1; end

dt=diff(pth.t);
dx=diff(pth.x)-pth.xstep; % how the point moved (~optical flow)
dy=diff(pth.y_shift); 

jumps = dt~=0;
dt=dt(jumps); dx=dx(jumps); dy=dy(jumps); 

xvel = [min(dx./dt) max(dx./dt)];
yvel = [min(dy./dt) max(dy./dt)];

yshift = [min(pth.y_shift) max(pth.y_shift)];


% =================== %
%    eval_refpoint    %
% =================== %

function ref_points=eval_refpoint(ref_points,mov,xwnd_rad)
% ref_points=eval_refpoint(ref_points,mov,xwnd_rad);
% Change all +Inf and -Inf (meaning boundary) to numeric values
%
% ref_points (2xP)
% mov - movie structure (load_movie.m)

x_min = 1+xwnd_rad;
x_max = mov.framesize(2)-xwnd_rad-1;
t_min=1; 
t_max=mov.nframes;

ref_points(1,ref_points(1,:)==-inf) = x_min;
ref_points(1,ref_points(1,:)==+inf) = x_max;
ref_points(2,ref_points(2,:)==-inf) = t_min;
ref_points(2,ref_points(2,:)==+inf) = t_max;


% ======================= %
%    expand_multipoint    %
% ======================= %

function indices = expand_multipoint(mpoint, sz)
% indices = expand_multipoint(mpoint, sz);
% Return indices of all points described by "multipoint".
% "Multipoint", e.g.: [1 nan 2], i.e., 2nd dimension can take any value.
%
% mpoint (Dx1) - Multipoint
% sz (1xD) - Dimensions for indexing
%
% indices (1xN)

exp_sz=sz;
isgiven=isfinite(mpoint);
exp_sz(isgiven)=1;

indices=zeros(1,prod(exp_sz));
for d=length(mpoint):-1:1;
    % Determine values
    if isgiven(d);
        exp_val=mpoint(d);
    else
        exp_val=1:exp_sz(d);
    end
    % Copy in correct order
    resh_sz=ones(size(exp_sz)); resh_sz(d)=exp_sz(d);
    rep_sz=exp_sz; rep_sz(d)=1;
    exp_val=repmat(reshape(exp_val,resh_sz),rep_sz);
    exp_val=exp_val(:)';
    % Add to indices
    indices = (exp_val-1) + sz(d)*indices;
end

indices=1+indices;


% =================== %
%    forest_struct    %
% =================== %

function forest=forest_struct(final)
% forest=forest_struct(final)
% Converts array to forest structure. 
% Optimized for tracing paths (path_in_forest.m). 
%
% INPUT
% final (NxK, K>=4) [dist node parent arc_weight ....]
%          - working array of 'dijkstra_online_core.m'
%
% OUTPUT
% 
% All nodes in the forest are renumbered to 1:N. All the arrays of the 'forest'
% have entries corresponding to these numbers.
% 
% forest (1x1) - structure with fields:
%   alias (1xN) - forest.alias(k) is the "real name" of the node 'k' (its number in the
%                 input 'final' array)
%   parent (1xN) - forest.parent(k) is the parent node of the node 'k' 
%                 parent = -1 if k is one of the start_nodes; 
%                 parent = NaN if the node is not in the tree
%   path_cost (1xN) - total cost of the best path from the source (one of start_nodes) 
%                     to the node (depends on the path_norm_type)
%   arc_weight (1xN) - weight of the arc from the parent to the node


% Find aliases for all nodes in the forest
[forest.alias,ind_in_f,ind_in_al]=unique(final(:,2:3));
forest.alias=forest.alias(:)';
% Rename nodes
final(:,2:3)=reshape(ind_in_al,size(final(:,2:3)));

%   init with "nan"
max_node=length(forest.alias);
forest.parent=repmat(nan,1,max_node);
forest.path_cost=repmat(nan,1,max_node);
forest.arc_weight=repmat(nan,1,max_node);
%   put values
forest.parent(final(:,2))=final(:,3);
root_symb=find(forest.alias==-1); 
forest.parent(root_symb)=root_symb(1); % remove 'nan's from the symbolic "-1" node
forest.path_cost(final(:,2))=final(:,1);
forest.arc_weight(final(:,2))=final(:,4);


% ===================== %
%    gaussian_window    %
% ===================== %

function [win,sigma] = gaussian_window(rect_length,win_length)
% [win,sigma] = gaussian_window(rect_length,win_length = automatic);
% Gaussian window function for given window size (not standard deviation).
% 
% Standard deviation of the gaussian (sigma) and size of the window
%    are chosen automatically.
%
% If win_length is not defined, or empty, or <=0, it is computed automatically.

THROW_OUT = .1; % Values < THROW_OUT*max_value are thrown out

sigma = rect_length/2/sqrt(2);
if ~exist('win_length','var'); win_length = []; end
if isempty(win_length) || (win_length<=0);
    N = max(3, round(sigma*sqrt(2*log(1/THROW_OUT)))*2+1);
else
    N = win_length;
end
x = 1:N;
middle = 1+(N-1)/2;
win = exp(-.5*((x-middle)./sigma).^2);
win = win./sum(win);


% =============== %
%    get_frame    %
% =============== %

function img=get_frame(mov, frame,flag)
% frame = get_frame(mov, frame[, flag=0])
%
% Returns video frame(s) with given index
%
% IN:
%       mov -  is structure returned by function 'load_movie'
%       frame - 1xK list of frame numbers to get
%       flag - 1x1, some options: 
%           0(default) - use cache
%           1 - reset cache
%           2 - skip cache use
%           3 - test mode
%
%
% OUTPUT:
%       img - HxWxCxT    one or more images

if ~exist('flag','var'); flag = 0; end
if isequal(mov.input_type,'array') 
    % Bypass cache when movie is in memory
    flag = 2;
end

global GET_FRAME_CACHE;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Load into cache and obtain indices    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if flag==0, % Use cache
    ind = load_frame(mov, frame);
elseif flag==1, % Reset cache
    init_movie_cache(mov);
    ind = load_frame(mov, frame);
elseif flag==2, % Directly read the frames --- bypass cache
    img=actually_get_frame(mov, frame);
    return;
elseif flag==3                       % test this function
    N=mov.nframes;
    F_MAX = min(N,20); % how many frames try to read
    disp 'Testing the cache...'
    % Make sure cache is initialized, without forcing initialization
    get_frame(mov, max(1,fix(N*rand(1,5))));
    % Get size of the cache
    K = size(GET_FRAME_CACHE.data, 2);
    for i=1:30
        if rand<.8 % try less than in the cache
            F = min(K,F_MAX);
        else
            F = F_MAX;
        end
        n_frames = max(1,round(rand*F));
        ind = max(1,round(N*rand(1,n_frames)));
        % Call load_frame and take directly from cache
        ind_in_cache = load_frame(mov,ind);
        fr_1 = GET_FRAME_CACHE.data(:,:,:,ind_in_cache);
        % Call get_frame with cacheing
        fr_2 = get_frame(mov, ind);
        % Call get_frame without cacheing
        fr_3 = get_frame(mov,ind,2);
        % Check
        if ~isequal(fr_1,fr_2,fr_3)
            GET_FRAME_CACHE
            error 'Got bad results'
        end
    end
    disp '    okay.'
    disp 'timing...'
    tic
    for i=1:min(N-4, 40); img=get_frame(mov, i:i+4); end;
    t1=toc;
    disp '------------------'
    tic
    for i=1:min(N-4, 40); img=get_frame(mov, i:i+4,2); end;
    t2=toc;
    fprintf('Speedup is %.2f times (%.1f vs %.1f)\n', t2/t1, t1,t2);
    return;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Return frames from cache     %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if isempty(ind); img = []; return; end

% Allocate images
img_sz=mov.framesize; 
if length(img_sz)<3; img_sz(3)=1; end; img_sz(4)=length(frame);
img=uint8(0); 
img(img_sz(1),img_sz(2),img_sz(3),img_sz(4))=0;

% Copy from cache
img = GET_FRAME_CACHE.data(:, :, :, ind);


% ====================== %
%    graph_func_movie    %
% ====================== %

function [adj_nodes,weights,steps]=graph_func_movie(node,nodes_exclude, mov,mask3D,xwnd, xrange,trange,yrange, ysubwnd)
% [adj_nodes,weights,steps]=graph_func_movie(node,nodes_exclude, mov,mask3D,xwnd, xrange,trange,yrange, ysubwnd)
% Computes adjacent arcs and their weights for the given movie. 
% /Helper function for dijkstra_online_core()./
% Arcs measure how consistent are pairs of space-time(-yshift) 
%     points to be put into mosaic. 
% Arcs will be only in the given 3D mask. Points too close to 
%     the boundary will be ignored.
% 
% INPUT:
%   node (1x1)
%   nodes_exclude (1xK) - disregard arcs to these nodes (e.g., the ones 
%                   that are already in the tree in dijkstra algorithm)
%   mov - movie structure (see load_movie.m)
%   mask3D (XxTxYshifts) - 3-dimensional mask of nodes which should be 
%                   in the graph. 
%
%   xwnd (1x1) - window size in X-dimension
%   xrange, trange, yrange (1xA,1xB,1xC) - increments to x, t, y-shift 
%                   values to generate arcs
%   ysubwnd (1x1) - size of comparison subwindow in Y-dimension: 
%                   each column is divided into subwindows, cost = maximum 
%                        over average differences in these subwindows;
%                   E.g.:
%                   ysubwnd==1: strict match over Y
%                   ysubwnd==Inf: average match over Y
%
% OUTPUT
%   adj_nodes, weights (1xN) - nodes, adjacent to the given one, and 
%              weights of corresponding arcs
%   steps (1xN) - distance in mosaic (in pixels) represented by each arc
%              (equal to 1).

% Empty output - in case there are no arcs 
adj_nodes=[]; weights=[]; steps=[]; 

max_y_shift=size(mask3D,3)-1;
ywnd=mov.framesize(1)-max_y_shift;
ysubwnd=max(1,min(ywnd,round(ysubwnd)));
ysuboff=floor(rem(ywnd,ysubwnd)/2);

xstep=1;
sz=size(mask3D); sz=sz(1:2); % [X T]

% x,t,y_shift of the given node
node_x= 1+mod(node-1,sz(1));
node_t= 1+mod(floor((node-1)/sz(1)),sz(2));
node_ysh= 1+floor((node-1)/prod(sz));

% Check if the node is too close to boundary or not in mask
if (node_x-xwnd<1) || (node_x+xwnd>sz(1)) || ~mask3D(node_x,node_t,node_ysh);
    return;
end
    
% Determine if shifts inside frame are required - we can save on computing them 
i_t_zero=find(trange==0);
IN_FRAME_SHIFT=~isempty(i_t_zero);
if IN_FRAME_SHIFT; trange(i_t_zero)=[]; end % save on computation

% --------------------------------
% ---- Determine to-points -------
% x and t for points in mask and in range
poss_ts=node_t+trange; 
poss_ts=poss_ts(poss_ts>=1 & poss_ts<=sz(2));
xtymask = mask3D(:,poss_ts,:);

if ~isempty(xrange);
    xmask=false; xmask(sz(1),1)=0;
    poss_xs=node_x+xrange+xstep; % we shift arcs, and xrange is for corresponding points
    poss_xs=poss_xs(poss_xs-xwnd-xstep>=1 & poss_xs+xwnd-xstep<=sz(1));
    xmask(poss_xs)=1;
    xmask=xmask(:,ones(1,size(xtymask,2)),ones(1,size(xtymask,3)));
    % Combine x- and t-masks
    xtymask=xmask & xtymask;
end
if ~isempty(yrange);
    ymask=false; ymask(1,1,max_y_shift+1)=0;
    poss_ysh=node_ysh+yrange;
    poss_ysh=poss_ysh(poss_ysh>=1 & poss_ysh<=max_y_shift+1);
    ymask(poss_ysh)=1;
    ymask=ymask(ones(1,size(xtymask,1)),ones(1,size(xtymask,2)), :);
    % Combine it with xt-mask
    xtymask=ymask & xtymask;
end

inds=find(xtymask);
if isempty(inds); return; end
xs= 1+mod(inds-1,size(xtymask,1));
ts= 1+mod(floor((inds-1)./size(xtymask,1)),size(xtymask,2));
ysh= 1+floor((inds-1)./(size(xtymask,1)*size(xtymask,2)));
ts=poss_ts(ts);
%   make t first
[ts,inds]=sort(ts);
xs=xs(inds); ysh=ysh(inds);
xs=xs(:); ts=ts(:); ysh=ysh(:);

% Transform coordinates into node indices
adj_nodes = xs+sz(1)* (ts-1+ sz(2)*(ysh-1));

% Remove unnecessary nodes
if ~isempty(nodes_exclude)
    %  "ismembc" is a C routine used in "ismember" function
    %  calling it directly is faster:
    nodes_exclude=sort(nodes_exclude);
    to_keep=~ismembc(adj_nodes,nodes_exclude);
    %  instead of:
    %     to_keep=~ismember(adj_nodes,nodes_exclude);
    adj_nodes=adj_nodes(to_keep);
    xs=xs(to_keep); ts=ts(to_keep); ysh=ysh(to_keep); 
end

if isempty(adj_nodes); return; end

%
% Now   xs,ts,ysh  (Kx1) are coordinates of 
% all to-points for the given node
% !!! ts are sorted!!!!
%_____________________________________


% ====== Load all images at once 
% ts_unique=poss_ts(find(any(xtmask,1))); - may be not good because of nodes_exclude
% ts_unique=unique(ts);
% A faster version for t_unique (uses the fact that "ts" is sorted):
ts_unique=[];last_t=-1; for k=1:length(ts); if ts(k)~=last_t;last_t=ts(k); ts_unique(end+1,1)=last_t; end; end
ts_unique=unique([ts_unique; node_t]);
% Load images
imgs=get_frame(mov,ts_unique);
% Delete non-needed x range
min_x=max(1,min(node_x,min(xs)-abs(xstep))-xwnd);
max_x=min(sz(1),max(node_x,max(xs)+abs(xstep))+xwnd);
imgs=imgs(:,min_x:max_x,:,:);
imgs=double(imgs);
% if size(imgs,3)>1; imgs(:,:,2)=2*imgs(:,:,2); end
% Simulate x-values (NOTE: adj_nodes are already computed on correct xs)
xs_eff=xs-min_x+1;
node_x_eff=node_x-min_x+1;

% From-window - shifts from the current node
wnd_1=[-(xwnd-1) 0];
wnd_2=[1 xwnd];
wnd_1=sort(sign(xstep)*wnd_1);
wnd_2=sort(sign(xstep)*wnd_2);

% Cut out the from-image
im_from=imgs(node_ysh:node_ysh+end-max_y_shift-1,:,:,find(ts_unique==node_t));
im_from_1=im_from(:,node_x_eff+wnd_1(1):node_x_eff+wnd_1(2),:);
im_from_2=im_from(:,node_x_eff+wnd_2(1):node_x_eff+wnd_2(2),:);

% To-window
wnd_1=wnd_1-xstep; wnd_2=wnd_2-xstep;

% --------------------------------------------
% ----- Iterate over all adjacent nodes ------
adj_i=1; 
n_adj=length(xs);
weights=zeros(n_adj,1);
while adj_i<=n_adj
    t=ts(adj_i);
    % Extract one frame
    im_to=imgs(:,:,:,find(ts_unique==t));
    % Do for all adjacent nodes with current t (ts should be sorted!!!)
    while (adj_i<=n_adj) && (t==ts(adj_i));
        x=xs_eff(adj_i); y_shift=ysh(adj_i);
        % Cut out the to-image
        im_to_1=im_to(y_shift:y_shift+end-max_y_shift-1,x+wnd_1(1):x+wnd_1(2),:);
        im_to_2=im_to(y_shift:y_shift+end-max_y_shift-1,x+wnd_2(1):x+wnd_2(2),:);
        
        % Compare
        diff_1=(im_to_1-im_from_1).^2;
        diff_2=(im_to_2-im_from_2).^2;
        % Accumulate 
        diff_1=sum(sum(diff_1,2),3);
        diff_1=[0;cumsum(diff_1,1)];
        diff_1=diff_1(1+ysuboff+ysubwnd:ysubwnd:end)-diff_1(1+ysuboff:ysubwnd:end-ysubwnd);
        diff_1=max(diff_1);
        diff_2=sum(sum(diff_2,2),3);
        diff_2=[0;cumsum(diff_2,1)];
        diff_2=diff_2(1+ysuboff+ysubwnd:ysubwnd:end)-diff_2(1+ysuboff:ysubwnd:end-ysubwnd);
        diff_2=max(diff_2);
        
        % Compute the weight
        weights(adj_i)=min([diff_1,diff_2]);
        
        % Go to the next adjacent node
        adj_i=adj_i+1;
    end
end % for all adjacent nodes
% --------------------------------------------


% Add "trivial" shift to the neighboring point
if IN_FRAME_SHIFT % add only one arc - to the neighboring point in the frame
    neighb_x=node_x+xstep;
    if (neighb_x>=1) && (neighb_x<=sz(1)) && mask3D(neighb_x,node_t,node_ysh)
        neighb_ind=neighb_x+sz(1)* (node_t-1+ sz(2)*(node_ysh-1));
        adj_nodes=[adj_nodes;neighb_ind];
        weights=[weights;0];
    end
end

% Normalize weights (make them "per pixel and color channel")
n_pixels=xwnd*ysubwnd;
if length(mov.framesize)>2; n_pixels=n_pixels*mov.framesize(3); end
weights=(1/n_pixels).*weights;

% Return trivial 'steps' vector
steps=ones(size(adj_nodes));


% ====================== %
%    init_movie_cache    %
% ====================== %

function init_movie_cache(mov)
% init_movie_cache(mov)
% Initialize movie cache (used by get_frame.m and load_frame.m).

global GET_FRAME_CACHE
MAX_AGE = 9e99; % use instead of inf
MEM_LIMIT = 50*2^20; % 50 Mb

K=min(mov.nframes, fix(MEM_LIMIT/prod(mov.framesize))); % initial cache size
% fprintf('\nInitializing movie cache to %d frames\n\n',K);
GET_FRAME_CACHE.signature=mov.signature;
GET_FRAME_CACHE.data=...
    repmat(uint8(0), [mov.framesize K]);      % the data we have
GET_FRAME_CACHE.indices=zeros(1,K);  % the indices of the data
GET_FRAME_CACHE.age=repmat(MAX_AGE,[1 K]);      % age (unused) of each entry
% age is the number of calls without access to the frame
% age is 0 when the frame have just been accessed
% GET_FRAME_CACHE.stat=[0 0];          % [saved loaded]
%    % Statistics for every frame:
%    GET_FRAME_CACHE.loaded_from_disk=zeros(1,mov.nframes);
%    GET_FRAME_CACHE.loaded_from_memory=zeros(1,mov.nframes);


% ================ %
%    load_frame    %
% ================ %

function frame_in_cache = load_frame(mov, frame)
%
% frame_in_cache = load_frame(mov, frame)
% Load video frame(s) to memory cache.
% For use with get_frame.m
%
% Cache is linked to a specific movie. If frames from a different
%   movie are requested, cache is automatically reset. 
%
% 'frame' array can contain multiple instances of the same frame. 
%   Each frame will be loaded once.
%
% IN:
%   mov - 1x1 structure returned by function 'load_movie'
%   frame - 1xF list of frame numbers to get
%
% OUTPUT:
%   frame_in_cache (1xF) - indices of the loaded frames in cache
%           (requested frames are GET_FRAME_CACHE.data(:,:,:,frame_in_cache))

global GET_FRAME_CACHE

% Reset cache in the following cases:
%    no global variable GET_FRAME_CACHE
%    for new movie (check signature)
if isempty(GET_FRAME_CACHE) ...
        || ~isequal(GET_FRAME_CACHE.signature, mov.signature)
    init_movie_cache(mov);
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     Load all required images into the cache     %
%     replacing old (least accessed) ones         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% ====== Deal with non-unique frames in 'frame' ===== %
% find unique frames (~4 times faster than unique())
ts = sort(frame);
m = diff(ts)~=0; m(1,end+1) = true;
ts = ts(m);

% ========================================== % 
%  See which frames we actually need to get  % 
% ========================================== % 
%   SLOW VARIANT (approximately equivalent):
%        [have_nums have_ind have_loc]=intersect(frame, GET_FRAME_CACHE.indices);
%        [get_num get_ind]=setdiff(frame, GET_FRAME_CACHE.indices);
%   FAST VARIANT:
[indices_sorted,srt_idx]=sort(GET_FRAME_CACHE.indices);
have_loc=ismembc2(ts,indices_sorted); % numbers
have_ind=have_loc>0;
have_loc=srt_idx(have_loc(have_ind)); % index to the unsorted indices
have_nums=ts(have_ind);
get_ind=~have_ind;
get_num=ts(get_ind);
% ===== Now: =====
% have_ind (1xF logical) - ts(have_ind) can be taken from cache
% have_nums (1xK) - frame numbers, corresponding to have_ind
% have_loc (1xK) - indices in cache, which can be used
% get_ind (1xF logical) - ts(get_ind) need to be read from disk
% get_num (1xK) - frame numbers, corresponding to get_ind

% ========================================== % 
% Check if we need to enlarge cache (in case too many frames are
% requested at once)
enlarge_cache = max(0, numel(get_num) - ...
    (numel(GET_FRAME_CACHE.indices)-numel(have_nums)));
if enlarge_cache>0;
    sz = size(GET_FRAME_CACHE.data);
    GET_FRAME_CACHE.data(sz(1), sz(2), sz(3), sz(4)+enlarge_cache) = 0;
    GET_FRAME_CACHE.indices(sz(4)+enlarge_cache) = 0;
    GET_FRAME_CACHE.age(sz(4)+1:sz(4)+enlarge_cache) = MAX_AGE;
end

% ========================================== % 
%      Update age of the current cache       %
% ========================================== % 
GET_FRAME_CACHE.age = 1 + GET_FRAME_CACHE.age;
GET_FRAME_CACHE.age(have_loc) = 0;
% /now all frames that in cache which are not requested have age>0/

% Where to find requested frames in cache:
ts_cache = zeros(1,numel(ts));
ts_cache(have_ind) = have_loc;

% ============================== % 
%      Read missing frames       %
% ============================== % 
if ~isempty(get_num)
    % Find the least used ones that should be replaced:
    % (with maximal age)
    [Y I]=sort(GET_FRAME_CACHE.age);
    replace_these=I(end-numel(get_num)+1:end);

    % Read new data from disk and fill in the cache
    GET_FRAME_CACHE.data(:,:,:,replace_these)=actually_get_frame(mov, get_num);
    GET_FRAME_CACHE.indices(replace_these)=get_num;
    GET_FRAME_CACHE.age(replace_these)=0;
    
    % Where to find loaded frames
    ts_cache(get_ind) = replace_these;
end

% ===== Return position in cache for 'frame' ====== %
ts_map = 0; ts_map(1, max(ts)) = 0;
ts_map(ts) = ts_cache;
frame_in_cache = ts_map(frame);


% ================ %
%    load_movie    %
% ================ %

function mov=load_movie(basename, range)
%
%INPUT:
%   basename - banana.avi or C:a/adsf/*.jpg
%   range - which images to load. Only valid for list of images.
% 
% OUTPUT:
%   mov - structure
%      .input_type = 'avi' or 'images' (sequence of images) or 'array' (matlab array)
%      .framesize   1x3 size of the frame, [Y X n_colors]
%      .nframes
%      .signature - unique identifier of the movie (used in cacheing)
%    For mov.input_type = 'avi':
%      .files - name of avi file
%      .path_prefix /OPTIONAL/ - if present, will be added before file name
%      .framerange (1xR) /OPTIONAL/ - frames that are considered (as 'range' argument of this function)
%    For mov.input_type = 'images':
%      .files - cell array of frames of the video
%      .path_prefix /OPTIONAL/ - if present, will be added before file name
%    For mov.input_type = 'array':
%      .video_data - Y x X x n_colors x T - video volume

REMEMBER_PATH = true;

if ischar(basename); % path to the video is given
    if findstr(basename, '.avi')
        mov.files=basename;
        mov.input_type='avi';
        img=aviread(basename, 1);
        img=img.cdata;
        mov.framesize = [size(img) repmat(1,1,3-length(size(img)))];
        info=aviinfo(basename);
        if ~exist('range','var');
            mov.nframes = info.NumFrames;
        else
            mov.framerange=range(:)';
            mov.nframes = length(mov.framerange);
        end

        % Compute unique movie signature
        mov.signature=cat(2, mov.framesize,mov.nframes,double(cd),...
            double(basename));
        
        % 
        if REMEMBER_PATH; mov.path_prefix = [cd '/']; end
        return
    else
        % We expect  C:/..../*.jpg
        D=dir(basename);
        
        % We need to add the prefix of the path...
        ind=find(basename=='/' | basename=='\');
        if ~isempty(ind), basename=basename(1:ind(end));
        else basename=''; end
        
        mov.files={D(:).name};
        mov.files=sort(mov.files);
        for i=1:length(mov.files), mov.files{i}=[basename mov.files{i}]; end
        
        % clip to the desired range
        if exist('range','var'), mov.files=mov.files(range); end
        
        mov.input_type='images';
        img=imread(mov.files{1});
        mov.framesize = [size(img) repmat(1,1,3-length(size(img)))]; 
        mov.nframes = length(mov.files);
        
        % Determine read-function for this type of images
        % ALL IMAGE FILES MUST BE IN THE SAME FORMAT !!!
        imfi=imfinfo(mov.files{1});
        imfo = imformats(imfi.Format);
        mov.read_func=imfo.read;
        
        % Compute unique movie signature
        mov.signature=cat(2, mov.framesize,mov.nframes,double(cd),...
            double(basename));

        if REMEMBER_PATH; mov.path_prefix = [cd '/']; end

        return;
    end
else % array itself is given
    if exist('range','var');
        mov.video_data = basename(:,:,:,range);
    else
        mov.video_data = basename;
    end
    mov.input_type = 'array';
    sz = size(mov.video_data); sz = [sz ones(1,4-numel(sz))];
    mov.framesize = sz(1:3);
    mov.nframes = sz(4);
    
    % Movie signature: make it unique for sure
    mov.signature = clock; 
end


% ===================== %
%    mkdir_overwrite    %
% ===================== %

function new_dir = mkdir_overwrite(new_dir)

CONFIRMATION = 0;

if exist(new_dir,'dir') % need to remove
    % Ask confirmation
    if CONFIRMATION
        conf = input(['Remove     "' new_dir '"    ? (enter = yes)'],'s');
        if ~isempty(conf)
            warning('User aborted directory removal');
            new_dir = [];
            return;
        end
    end
    % Remove
    if ~rmdir(new_dir,'s'); % remove the directory with its contents
        warning('Directory could not be removed');
        new_dir = [];
        return;
    end
end

% Make the directory
if ~mkdir(new_dir); % error creating directory
    new_dir = [];
end

% ============================ %
%    mosaic_image_from_path    %
% ============================ %

function img = mosaic_image_from_path(pth, mov, BCK_COLOR,SHOW_ARC_COST)
% img = mosaic_image_from_path(pth, mov, BCK_COLOR=black/0/,SHOW_ARC_COST=1);
% Create mosaic image given path in XT. 
% Path may include shifts in Y.
% Optionally, arc weights ("errors") are depicted at the bottom of the mosaic image
%   as different intensities of red
% 
% INPUT
%   pth - structure with fields:
%     x,t, y_shift (1xN)
%   mov - movie structure (from load_movie.m)
%   BCK_COLOR - background color; can be either 1x1 (intensity) or 1x3 (RGB)
%   SHOW_ARC_COST (1x1) logical - whether to include dislay of path arc weights into the image
% OUTPUT
%   img - mosaic image (uint8)

% ----------------------------
%    Deal with parameters
% ----------------------------
if ~exist('BCK_COLOR','var'); BCK_COLOR=uint8([0 0 0]); end
if ~exist('SHOW_ARC_COST','var'); SHOW_ARC_COST=true; end
if length(BCK_COLOR)==1; BCK_COLOR=BCK_COLOR([1 1 1]); end
BCK_COLOR=reshape(BCK_COLOR,1,1,3);
BCK_COLOR=uint8(BCK_COLOR);
if mov.framesize(3)==1; % grayscale
    BCK_COLOR=rgb2gray(BCK_COLOR);
end
COST_COLOR=[255 0 0]; % red
if norm(double(COST_COLOR(:))-double(BCK_COLOR(:)))<5; % too similar color
    COST_COLOR=[0 255 0]; % green
end

% ----------------------------
%    Prepare data
% ----------------------------
x=pth.x; t=pth.t; 
if isfield(pth,'y_shift');
    y_shift=pth.y_shift;
else
    y_shift=zeros(size(x));
end
y_shift=y_shift-min(y_shift); % ensure it starts from zero
xstep=1;
n_points = length(x);
Y = mov.framesize(1);
mos_sz = mov.framesize;
mos_sz(1)=mos_sz(1)+max(y_shift);
mos_sz(2)=1+xstep*(n_points-1); % each point except the first corresponds to xstep points
img = repmat(BCK_COLOR,mos_sz(1:2));

% ---------------------------------
%   Put columns from video frames
%        into mosaic
% ---------------------------------
tt = unique(t);
for tti = 1:length(tt); this_t = tt(tti);
    this_inds = find(t==this_t);
    orig_frame = get_frame(mov,this_t);
    for k=1:length(this_inds);
        this_i=this_inds(k);
        this_ysh=y_shift(this_i);
        this_x=x(this_i);
        this_col=1+xstep*(this_i-1);
        if this_i==1;
            img(end-this_ysh-Y+1:end-this_ysh, this_col,:)...
                = orig_frame(:,this_x,:);
        else
            img(end-this_ysh-Y+1:end-this_ysh, this_col-xstep+sign(xstep):sign(xstep):this_col,:)...
                = orig_frame(:,this_x-abs(xstep)+1:1:this_x,:);
        end
    end;
end

% ---------------------------------
% Add arc cost (weight) if needed
% ---------------------------------
if SHOW_ARC_COST && isfield(pth,'arc_weight');
    cost_sz=mos_sz; cost_sz(1)=max(10,ceil(.05*mos_sz(1)));
    cost_im=repmat(BCK_COLOR,ceil(cost_sz(1:2).*[.1 1]));
    %   normalize cost
    cost=pth.arc_weight; 
    if max(cost)~=min(cost); 
        cost=(cost-min(cost))./(max(cost)-min(cost));
    else
        cost(:)=0;
    end
    cost=[0 cost];
    %   put into image
    cost=repmat(cost(:)',[1 1 length(COST_COLOR)]).*repmat(reshape(COST_COLOR,1,1,[]),[1,length(cost),1]);
    cost=repmat(cost,[ceil(cost_sz(1)*.9) 1]);
    cost_im=uint8([cost_im;round(cost)]);
    %    add to the mosaic
    img=[img;cost_im];
end


% ================= %
%    mosaic_main    %
% ================= %

function pyr = mosaic_main(movie,...
    ref_points,mask,comp_params,exp_params,dyn_params)
% Compute mosaic from movie using Dijkstra algorithm with video pyramid.
%
% INPUT
%   movie - a string (e.g., 'orig/*.jpg' or 'movie.avi')
%   ref_points (2xP, P>=2) - path will go through these points.
%          ref_points [x0 x1 ... xP; t0 t1 ... tP]
%          Can contain +-Inf and NaN values, treated specially 
%               (see mov_best_path_thru_points() for description).
%   mask (XxT logical or empty) - mask in xt-coordinates, where
%          the path is constrained to be located. Empty is treated
%          as no constraints (full mask of true's).
%   comp_params, exp_params, dyn_params (struct) - algrorithm parameters.
%          See pyramid_level_from_scratch() for description.
% OUTPUT
%   pyr (cell 1xL) - intermediate results for each level of the pyramid
%      (see pyramid_level_from_scratch() for description.)
%
% USAGE
% pyr = mosaic_main(movie, dominant_pan_dir, ...
%             ref_points,mask,comp_params,exp_params,dyn_params)


% Load original movie
movies=load_movie(movie);

% =========== Prepare movie pyramid if needed ============== %
fprintf('*** Creating pyramid of subsampled movies ***\n');
movies=movie_pyramid_create(movies);
fprintf('\n');

% ================================================================ %
% ===                 Run the pyramid iterations                == %
% ================================================================ %
pyr=pyramid_init(ref_points,movies,mask,comp_params,exp_params,dyn_params);
pyr=pyramid_run(pyr);


% ========================== %
%    mov_best_path_params    %
% ========================== %

function params=mov_best_path_params(comp_params,dyn_params,exp_params)
% params - see mov_best_path_thru_points.m

% Bounds for x and y ranges
xbnd = dyn_params.xvel + exp_params.xvel_wnd*[-1 1];
ybnd = dyn_params.yvel + exp_params.yvel_wnd*[-1 1];
xbnd(1)=floor(xbnd(1)); xbnd(2)=ceil(xbnd(2));
ybnd(1)=floor(ybnd(1)); ybnd(2)=ceil(ybnd(2));

% Parameters
params=comp_params;
params.xrange=xbnd(1):xbnd(2);
params.yrange=ybnd(1):ybnd(2);


% =============================== %
%    mov_best_path_thru_points    %
% =============================== %

function [pth,forest]=mov_best_path_thru_points(ref_points, mov,mask3D, params)
% [pth,forest]=mov_best_path_thru_points(ref_points, mov,mask3D, params)
% Compute best path for building mosaic from the given movie. 
% The path should pass through the given reference points.
% The path is restricted to go only in the given mask. 
% 
% INPUT
%  ref_points (2xN) [x;t] - path should pass through these x-t points.
%       Can contain +Inf/-Inf values (meaning boundaries);
%                   NaN - meaning any value (inside boundaries).
%       E.g.: [nan nan;-inf inf] means "from the first to the last frame".
%              /Shifts in y are determined automatically./
%  mov - movie structure (as in load_movie.m)
%  mask3D (X x T x Yshifts) logical - determines which points 
%       are allowed in the path. Can be used to remove unwanted 
%       regions or to constrain search when refining a known path.
%
%  params - struct with fields:
%    path_norm_type (1x1) - k (sum w^k, k integer) or Inf (max w) 
%          - operation used to compute cost of the path 
%            from weights 'w' of constituting arcs
%    xwnd (1x1) - window size in x-dimension; 
%    ysubwnd (1x1) - size of subwindow in y-dimension 
%       (cost function requires averages over these subwindows to match),
%       Inf means the whole column; the smaller ysubwnd, 
%                                   the more strict is matching.
%    xrange,trange,yrange (1xK,1xL,1xM) - allowed shifts in each dimension 
%           (xrange and yrange can be empty, meaning no restrictions).
%
% OUTPUT
%  pth - structure with fields
%    x, y, y_shift (1xP) - coordinates of the path nodes
%    arc_weight (1 x P-1) - weights of arcs in the path
% forest (1 x N-1 cell) - best-path forest for each segment, 
%    structures with fields: parent, arc_weight, path_cost

n_segments=size(ref_points,2)-1;
m_tmp=mask3D(1+params.xwnd:end-params.xwnd, :,:); 
n_nodes=sum(sum(sum(m_tmp)));
% % Uncomment for additional info
% fprintf('\t Computing best path through %d reference points, ',n_segments+1);
% fprintf(' mask contains %.1f%% points\n',100*n_nodes/numel(m_tmp));
% xr=length(params.xrange); tr=length(params.trange); yr=length(params.yrange);
% fprintf('\t\t');
% if xr*tr*yr>0;
%     fprintf('search range %dx%dx%d (expected %d arcs per node); ',yr,xr,tr,xr*tr*yr);
% end
% fprintf('window size (x) = %d; path norm = %d\n',params.xwnd,params.path_norm_type)

% Remove '+-Inf's from ref_points
ref_points=eval_refpoint(ref_points,mov,params.xwnd);
% Parameters for computing forest of best path (best_path_forest.m)
forest_params={n_nodes, params.path_norm_type, ...
        @graph_func_movie, mov,mask3D,...
        params.xwnd, params.xrange,params.trange,params.yrange, params.ysubwnd};

% ====================================
% == Find best path for every segment
pth=[];
% expand all 'NaN's
from_nodes=expand_multipoint([ref_points(:,1);nan],[size(mask3D) 1]); 
forest = cell(1, n_segments);
for r=1:n_segments;
    if n_segments>1; fprintf('\t *** Segment %d/%d ***\n',r,n_segments); end
    
    % expand all 'NaN's
    to_nodes=expand_multipoint([ref_points(:,r+1);nan],[size(mask3D) 1]); 
    
    % Find best path forest
    forest{r} = best_path_forest(from_nodes,to_nodes, forest_params{:});
    
    % Extract best path (closest to the to_nodes)
    [best_to_node,to_dist]...
        =nearest_node_in_forest(forest{r},to_nodes,size(mask3D));
    if to_dist>0;
        fprintf('\t *** Warning: have not found path to the given nodes, taking nearest path\n');
    end
    pth_segm=path_in_forest(forest{r},best_to_node);
    
    % Add computed segment to previous ones
    if isempty(pth); % 1st segment
        pth=pth_segm;
        pth=rmfield(pth,'path_cost');
        fprintf('\t *** Added %d nodes to the path ***\n',...
            length(pth_segm.nodes));
    else % add a new segment
        pth.nodes=[pth.nodes pth_segm.nodes(2:end)];
        pth.arc_weight=[pth.arc_weight pth_segm.arc_weight];
        fprintf('\t *** Added %d nodes to the path ***\n',...
            length(pth_segm.nodes)-1);
    end
    
    % Prepare for the next segment
    from_nodes=pth.nodes(end); % for the next segment
end

% == Create path structure
[pth.x,pth.t, y_shift] = ind2sub(size(mask3D),pth.nodes);
pth.y_shift=y_shift-1;



% **************************************************
function [node,min_dist]=nearest_node_in_forest(forest,to_nodes,dims)
% list all nodes in the forest
if isfield(forest,'alias')
    nodes_in_forest=forest.alias; 
    nodes_in_forest(nodes_in_forest==-1)=[];
else
    nodes_in_forest=find(~isnan(forest.parent));
end

[x,t,ysh]=ind2sub(dims,nodes_in_forest);
[x_to,t_to,ysh_to]=ind2sub(dims,to_nodes);

[k,dist]=dsearchn([x(:) t(:) ysh(:)], [x_to(:) t_to(:) ysh_to(:)]);
[min_dist,kk]=min(dist);

node=nodes_in_forest(k(kk));


% ========================== %
%    movie_pyramid_create    %
% ========================== %

function movies=movie_pyramid_create(mov,min_frame_dim)
% movies=movie_pyramid_create(mov,min_frame_dim=30)
%
% INPUT
%   mov (1x1) - movies struct (load_movie.m)
%   min_frame_dim (1x1) - frame dimensions in the smalles level of the pyramid will be not less than this value
% OUTPUT
%   movies (cell 1xL) - movies for pyramid levels

if ~exist('min_frame_dim', 'var'); min_frame_dim = 30; end

% Create pyramid
movies=spatial_pyramid_create(mov,'.',min_frame_dim);
movies=movie_pyramid_from_spatial(movies);

% Report
fprintf('Pyramid levels:');
for level=1:length(movies);
    fprintf(' %d) %dx%dx%d',level,movies{level}.framesize(1:2),movies{level}.nframes);
end
fprintf('\n');

% *******************************************************
function movs=spatial_pyramid_create(orig_mov,target_dir, MIN_SIZE)
%
% movs=spatial_pyramid_create(orig_mov,target_dir='.', MIN_SIZE)
% E.g.: spat_pyr=spatial_pyramid_create(load_movie('orig/*.jpg'),'.');
% Create pyramid of movies, frames at each level are half the size (in x and y)
%  of the size at the previous level.
% 
% INPUT
% orig_mov - movie structure (e.g., from load_movie.m)
% target_dir (string), e.g., 'my_subdir/'
% MIN_SIZE (1x1) - pyramid creation stop before image size goes below MIN_SIZE

fprintf('Original video: %03dx%03d, %d frames\n',...
    orig_mov.framesize(1:2),orig_mov.nframes);

if ~exist('target_dir','var'); target_dir='.'; end

target_dir=standardize_path(target_dir);
if ~exist(target_dir,'dir')
    if ~mkdir(target_dir)
        error(['Error creating output directory ' target_dir]);
    end
end

n_levels = 1+floor(log2(min(orig_mov.framesize(1:2))/MIN_SIZE));

filt=gaussian_window(3);

% Create output directories
if iscell(orig_mov.files)
    orig_pth=fileparts(orig_mov.files{1});
else
    orig_pth=fileparts(orig_mov.files);
end
orig_pth=standardize_path(orig_pth);
sz=orig_mov.framesize(1:2);
for l=2:n_levels;
    sz=ceil(sz./2);
    level_dir{l}=[orig_pth(1:end-1) sprintf('_%03dx%03d/',sz)];
    if isempty(mkdir_overwrite(level_dir{l})); error('Failed to create directory %s',level_dir{l}); end
end

% Resize and save images
START_CLOCK=clock;
T=orig_mov.nframes;
xt_sz=[orig_mov.framesize(2) T];
if length(orig_mov.framesize)>2; xt_sz=[xt_sz orig_mov.framesize(3)]; end
xtslice=repmat(uint8(0),xt_sz);
last_mess_len=0; 
for t=1:T;
    im=get_frame(orig_mov,t);
    xtslice(:,t,:)=reshape(uint8(median(double(im),1)),xt_sz([1 3]));
    for l=2:n_levels;
        im=imfilter(im,filt);
        im=im(1:2:end,1:2:end,:);
        out_fname=[target_dir level_dir{l} sprintf('%04d.jpg',t)];
        imwrite(im,out_fname,'Quality',100);
    end
    
    fprintf(repmat('\b',1,last_mess_len)); 
    last_mess_len=fprintf('Frame %d/%d, finishing in %s',t,T,seconds2str((T/t-1)*etime(clock,START_CLOCK)));
end
fprintf(repmat('\b',1,last_mess_len)); 
fprintf('Movie pyramid: finished in %s\n',seconds2str(etime(clock,START_CLOCK)));

movs=cell(1,n_levels);
movs{1}=orig_mov;
movs{1}.xtslice=xtslice;
for l=2:n_levels;
    movs{l}=load_movie([target_dir level_dir{l} '*.jpg']);
    movs{l}.xtslice=imresize(xtslice,[movs{l}.framesize(2) movs{l}.nframes]);
end

% *******************************************************
function pyr=movie_pyramid_from_spatial(spatial_pyr)
% pyr=movie_pyramid_from_spatial(spatial_pyr)
% Returns pyramid of movies, for given spatial movie pyramid
% (i.e. full-length movies subsampled in width and height only).
%
% INPUT
% Cell array of movie structures (see load_movie.m).
%
% OUTPUT
% Cell array of movie structures (see load_movie.m),
% constituting the pyramid

level=1;
x_level=1;
t_level=1;
pyr{1}=spatial_pyr{1};

% Equal subsampling in X and T
while x_level<length(spatial_pyr);
    x_level=x_level+1;
    t_level=t_level+1;
    level=level+1;
    pyr{level}=decimate_movie(spatial_pyr{x_level},2^(t_level-1));
end


% ==================== %
%    path_in_forest    %
% ==================== %

function pth = path_in_forest(forest,end_node)
%
% pth = path_in_forest(forest,end_node)
% Finds path from root to the given node.
%
% forest (1x1) - structure (see forest_struct.m)
%    It may or may not contain field "alias". 
%
% pth - struct with fields:
%    nodes (1xN) - nodes of the path
%    arc_weight (1xN-1) - weight of the arc
%    path_cost (1xN) - cost for each sub-path

DIRECT=~isfield(forest,'alias');
if ~DIRECT
    alias=forest.alias;
    end_node=find(alias==end_node);
end

P=forest.parent; % /struct fields are not accelerated by Matlab - important for the loop/
if isnan(P(end_node)); % No path
    pth=[];
    return;
end
n_nodes = sum(P>0); % maximal number of nodes in path

% Track path from end_node till the root (marked by -1)
curr_node = end_node;
if DIRECT; 
    reached_root= (curr_node<0);
else
    reached_root= (alias(curr_node)<0); 
end
pth_nodes = zeros(1,n_nodes);
cnt = 0;
while (~reached_root)&&(cnt<n_nodes); % avoid loops in bad trees
    cnt = cnt+1;
    pth_nodes(cnt) = curr_node;
    curr_node = P(curr_node);
    if DIRECT; 
        reached_root= (curr_node<0);
    else
        reached_root= (alias(curr_node)<0); 
    end
end

if ~DIRECT;
    curr_node=alias(curr_node);
end

% Check sanity and return results
if curr_node~=-1;
    error('Bad forest: a path did not lead back to root');
end

pth_nodes=pth_nodes(cnt:-1:1);
pth.arc_weight=forest.arc_weight(pth_nodes); 
pth.arc_weight(1)=[]; % remove dummy weight to the source
pth.path_cost=forest.path_cost(pth_nodes);

if ~DIRECT; pth_nodes=alias(pth_nodes); end
pth.nodes=pth_nodes;


% ================== %
%    pyramid_init    %
% ================== %

function pyr=pyramid_init(ref_points,movies,mask,comp_params,exp_params,dyn_params)
% pyr=pyramid_init(ref_points,movies,mask,comp_params,exp_params,dyn_params)
% Create array of level structures.
%
% INPUT
%  movies (cell 1xL) - movie structure for each level
%  ref_points, mask, comp_params, exp_params, dyn_params 
%         - all are given at the 1st (finest) level
%
% OUTPUT
%  pyr (cell 1xL) - structure with fields
%   ref_points (2xP)
%   mov (1x1) - movie structure
%   mask (XxT) logical
%   comp_params
%   exp_params
%   dyn_params
%   /see pyramid_level_from_scratch.m for description of parameter structs/

if length(exp_params.dilation_diam)==1; 
    exp_params.dilation_diam=exp_params.dilation_diam([1 1 1]); 
end

sz_1=[movies{1}.framesize(2) movies{1}.nframes movies{1}.framesize(1)];
n_levels=length(movies);
pyr = cell(1,n_levels);
for k=1:n_levels;
    sz=[movies{k}.framesize(2) movies{k}.nframes movies{k}.framesize(1)];
    % Just copy
    pyr{k}.mov=movies{k};
    if isempty(mask); pyr{k}.mask=mask; else
        pyr{k}.mask=imresize(mask,sz(1:2)); 
    end
    pyr{k}.comp_params=comp_params;
    pyr{k}.exp_params=scale_exp_params(exp_params,sz_1,sz);
    % Reference points
    pyr{k}.ref_points=scale_ref_points(ref_points,sz_1(1:2),sz(1:2));
    % Dynamics parameters
    pyr{k}.dyn_params=scale_dyn_params(dyn_params,sz_1,sz);
end

% *****************************8
function ref_points=scale_ref_points(ref_points,sz_from,sz_to)
% ref_points (2xP) can include NaN or +-Inf - they will not be scaled
% sz_from,sz_to (1x2)

fact=sz_to./sz_from;

% x
fin=isfinite(ref_points(1,:));
ref_points(1,fin) = min(sz_to(1),max(1,round(   fact(1)*ref_points(1,fin)   )));
% t
fin=isfinite(ref_points(2,:));
ref_points(2,fin) = min(sz_to(2),max(1,round(   fact(2)*ref_points(2,fin)   )));

% *****************************8
function dyn_params=scale_dyn_params(dyn_params,sz_from,sz_to)
% sz_from,sz_to (1x3)

fact=sz_to./sz_from;

dyn_params.xvel = fact(1)/fact(2) * dyn_params.xvel;
dyn_params.yvel = fact(3)/fact(2) * dyn_params.yvel;
dyn_params.max_y_shift = fact(3) * dyn_params.max_y_shift;

% *****************************8
function exp_params=scale_exp_params(exp_params,sz_from,sz_to)
% sz_from,sz_to (1x3)

MIN_WND=0;
fact=sz_to./sz_from;

exp_params.xvel_wnd = max(MIN_WND,ceil(fact(1)*exp_params.xvel_wnd));
exp_params.yvel_wnd = max(MIN_WND,ceil(fact(3)*exp_params.yvel_wnd));


% =============================== %
%    pyramid_level_from_coarse    %
% =============================== %

function lev=pyramid_level_from_coarse(lev,lev_coarse)
% lev=pyramid_level_from_coarse(lev)
% 
% lev_coarse - struct (see pyramid_level_from_scratch)
%  Fields needed:
%   mov (1x1) - movie structure (load_movie.m)
%   pth (1x1) - struct with fields x,t,y_shift
% 
% lev
%  Fields needed: 
%   ref_points
%   mov
%   mask
%   comp_params
%   exp_params
%
% lev
%  Fields will be added:
%   dyn_params
%   pth
%   forest
%   mosaic
%   mask3D

% Determine dynamics from the coarse path
% Assign lev.dyn_params
[xvel,yvel,yshift]=dynamics_from_path(lev_coarse.pth);
% fprintf('\t Dynamics of the coarse level path: %.1f <= dx/dt <= %.1f; %.1f <= dy/dt <= %.1f; %.1f <= y <= %.1f;\n',xvel,yvel,yshift);
sz=[lev.mov.framesize(2), lev.mov.nframes, lev.mov.framesize(1)];
sz_coarse=[lev_coarse.mov.framesize(2), lev_coarse.mov.nframes, ...
    lev_coarse.mov.framesize(1)];
fact=sz./sz_coarse;
lev.dyn_params.xvel=fact(1)/fact(2)*(xvel+[-1 1]) - [-1 1];
lev.dyn_params.yvel=fact(3)/fact(2)*(yvel+[-1 1]) - [-1 1];
lev.dyn_params.max_y_shift...
    =round((yshift(2)-yshift(1))*fact(3)+lev.exp_params.yvel_wnd);
% fprintf('\t Parameters chosen for the fine (current) level: x velocity %.1f:%.1f, y velocity %.1f:%.1f, max y shift %.1f\n',...
%     lev.dyn_params.xvel,lev.dyn_params.yvel,lev.dyn_params.max_y_shift);

% Create 3D mask
%    path from coarse
pth_fc=lev_coarse.pth;
pth_fc.x=fact(1)*pth_fc.x; pth_fc.x=min(sz(1),max(1,round(pth_fc.x)));
pth_fc.t=fact(2)*pth_fc.t; pth_fc.t=min(sz(2),max(1,round(pth_fc.t)));
pth_fc.y_shift=fact(3)*(pth_fc.y_shift-yshift(1));
pth_fc.y_shift=round(pth_fc.y_shift+lev.exp_params.yvel_wnd/2);
pth_fc.y_shift=min(lev.dyn_params.max_y_shift,max(1,pth_fc.y_shift));
%    interpolate between points
no_gaps=no_jumps_path([pth_fc.x(:),pth_fc.t(:),pth_fc.y_shift(:)+1]);
no_gaps=round(no_gaps);
%    put interpolated path to the mask
mask3D=repmat(false,[sz(1:2) lev.dyn_params.max_y_shift+1]);
mask3D(sub2ind(size(mask3D),no_gaps(:,1),no_gaps(:,2),no_gaps(:,3)))=1;
%    add boundaries from the coarse level at the t's of the relevant reference points
xwnd_fc=ceil(fact(1)*(lev_coarse.comp_params.xwnd+1)-1); % size of the coarse-level window at the fine level
rp=eval_refpoint(lev.ref_points,lev.mov,lev.comp_params.xwnd);
%       leave only near-border values
rel_rp=isnan(rp(:,1)) | rp(:,1)<=xwnd_fc | rp(:,1)>=sz(1)-xwnd_fc+1;
rp=rp(:,rel_rp);
%       consider separately left and right border
left_border=isnan(rp(1,:)) | rp(1,:)<=xwnd_fc;
brd_ts=rp(2,left_border); % t-coords of the points
if any(isnan(brd_ts)); brd_ts=1:sz(2); end %  all points
mask3D(1:max(1,xwnd_fc-lev.exp_params.dilation_diam(1)+1),unique(brd_ts))=1; % compensate for dilation
right_border=isnan(rp(1,:)) | rp(1,:)>=sz(1)-xwnd_fc+1;
brd_ts=rp(2,right_border); % t-coords of the points
if any(isnan(brd_ts)); brd_ts=1:sz(2); end %  all points
mask3D(min(end,end-xwnd_fc+lev.exp_params.dilation_diam(1)):end,unique(brd_ts))=1; % compensate for dilation
%  ***** Dilate *****
mask3D=imdilate(mask3D,ones(lev.exp_params.dilation_diam));
%    put in to the lev structure
lev.mask3D = mask3D;


% Create parameters
params=mov_best_path_params(lev.comp_params,lev.dyn_params,lev.exp_params);
% Compute best path
[lev.pth,lev.forest]=mov_best_path_thru_points(lev.ref_points, ...
    lev.mov,lev.mask3D, params);
% Generate mosaic image
lev.mosaic = mosaic_image_from_path(lev.pth, lev.mov);


% ********************************************
function coords_new = no_jumps_path(coords)

jumps = find(any(abs(diff(coords))>1,2));
coords_new=coords(1:jumps(1),:);

for jm = 1:length(jumps);
    i_from = jumps(jm);
    pnt_from = coords(i_from,:);
    pnt_to = coords(i_from+1,:);
    n_steps = max(abs(pnt_to-pnt_from))-1;
    coeff = (1:n_steps)./(n_steps+1); coeff=coeff(:);
    extra_coords=pnt_from(ones(n_steps,1),:).*(1-coeff(:,ones(length(pnt_from),1))) + pnt_to(ones(n_steps,1),:).*coeff(:,ones(length(pnt_from),1));
    coords_new=[coords_new;extra_coords];
    if jm<length(jumps);
        coords_new=[coords_new;coords(jumps(jm)+1:jumps(jm+1),:)];
    else
        coords_new=[coords_new;coords(jumps(jm)+1:end,:)];
    end
end % for all jumps


% ================================ %
%    pyramid_level_from_scratch    %
% ================================ %

function lev=pyramid_level_from_scratch(lev)
% lev=pyramid_level_from_scratch(lev)
% 
% lev - struct with fields
%   ref_points (2xP) 
%   mov (1x1) - movie structure (see load_movie.m)
%   mask (XxT) logical - mask in xt-slice
%   comp_params - computation parameters (see description 
%                        in mov_best_path_thru_points())
%     trange, path_norm_type,  xwnd, ysubwnd
%   dyn_params - dynamics parameters
%     xvel, yvel (1x2) [min max] - minimal and maximal shift 
%                            in x and y between neighboring frames
%     max_y_shift (1x1)  - maximal total (cumulative) shift in y 
%                             between two frames in the movie
%   exp_params - expansion parameters
%     xvel_wnd, yvel_wnd (1x1) - window to search around 
%                                the estimated velocity
%     dilation_diam (1x3) [x t y] - dilate the path mask when 
%                                refining from coarse path
%
% This function adds (or replaces) the following fields in lev:
%    pth (1x1) - struct with fields x,t,y_shift,xstep
%    forest (1 x P-1) - forest for each segment of the path 
%                       (see  forest_struct())
%    mosaic (YxMx3 or YxM) - mosaic image
%    mask3D (XxTxYshifts) logical

% Create 3D mask
mask=lev.mask;
if isempty(mask);
    mask=repmat(true,[lev.mov.framesize(2) lev.mov.nframes]);
end
lev.mask3D=repmat(mask,[1 1 ceil(lev.dyn_params.max_y_shift)+1]);

% Create parameters
params=mov_best_path_params(lev.comp_params,lev.dyn_params,lev.exp_params);
% Compute best path
[lev.pth,lev.forest]=mov_best_path_thru_points(lev.ref_points, ...
    lev.mov,lev.mask3D, params);
% Generate mosaic image
lev.mosaic = mosaic_image_from_path(lev.pth, lev.mov);


% ================= %
%    pyramid_run    %
% ================= %

function pyr=pyramid_run(pyr,level_from,level_to)
% pyr=pyramid_run(pyr,level_from=COARSEST,level_to=FINEST)
% Best path mosaic: iterate through pyramid levels.
%
% pyr (cell 1xL) - structure (see pyramid_init.m for description)

INTERMEDIATE_SAVE=false; % save results after each level
INTERMEDIATE_DISPLAY=false; % display results after each level

n_levels=length(pyr);
if ~exist('level_from','var'); level_from=n_levels; end
if ~exist('level_to','var'); level_to=1; end

PYR_CLOCK=clock;
fprintf('\n');
fprintf('[%s] Started computation of best path mosaic in %d levels\n\n',datestr(clock,21),n_levels);
for lev_no=level_from:-1:level_to;
    LEVEL_CLOCK=clock;
    fprintf('******** LEVEL %d ********\n',lev_no);
    fprintf('Video %dx%dx%d (YxXxT)',pyr{lev_no}.mov.framesize(1:2),pyr{lev_no}.mov.nframes);
    if (lev_no==n_levels) || ~isfield(pyr{lev_no+1},'pth');
        fprintf('; computing from scratch\n');
        pyr{lev_no}=pyramid_level_from_scratch(pyr{lev_no});
    else
        fprintf('; using estimation from level %d\n',lev_no+1);
        pyr{lev_no}=pyramid_level_from_coarse(pyr{lev_no},pyr{lev_no+1});
    end
    pyr{lev_no}.running_time=etime(clock,LEVEL_CLOCK);
    fprintf('Level %d completed in %s\n\n',lev_no,seconds2str(pyr{lev_no}.running_time));
    % Save and display results
    if INTERMEDIATE_SAVE;
        imwrite(pyr{lev_no}.mosaic,...
            sprintf('%d%02d%02d_%02d%02d%02d_mosaic%02d.jpg',...
            round(clock),lev_no),'Quality',98);
        save(sprintf('%d%02d%02d_%02d%02d%02d_level%02d.mat',...
            round(clock),lev_no),'pyr','lev_no');
    end
    if INTERMEDIATE_DISPLAY;
        figure; 
        imshow(pyr{lev_no}.mosaic);
        title(sprintf('Level %d',lev_no)); 
        drawnow
    end
end
fprintf('[%s] %d levels in %s\n',datestr(round(clock),21),...
    n_levels,seconds2str(etime(clock,PYR_CLOCK)));


% ================= %
%    seconds2str    %
% ================= %

function timestr=seconds2str(secs)
% timestr=seconds2str(secs)
% Converts seconds into time string:
%  5d 12:30:15 = 5 days 12 hours 30 minutes 15 seconds
%  4:15:05 = (0 days) 4 hours 15 minutes 5 seconds


% Days-hours-minutes-seconds vector
dhms=[fix(secs/86400) rem(fix(secs/3600),24) rem(fix(secs/60),60) rem(secs,60)];
dhms=round(dhms);

% String
if dhms(1)>0; % include days
    timestr=sprintf('%dd %d:%02d:%02d',dhms);
else % don't include days
    timestr=sprintf('%d:%02d:%02d',dhms(2:end));
end


% ====================== %
%    standardize_path    %
% ====================== %

function output_path = standardize_path(input_path)
% output_path = standardize_path(input_path)
%
% Makes path of the form: 'c:/dir/subdir/subsubdir/'
%  (with '/' as separator, and ending on '/')

% Replace '\' by '/'
output_path = strrep(input_path,'\','/');

% Remove all '/' at the end
if isempty(output_path); output_path = './'; return; end
while output_path(end)=='/'
    output_path(end) = [];
    if isempty(output_path); output_path = './'; return; end
end

% Add exactly one '/' at the end
output_path(end+1) = '/';