function [tidx, tmesh, tmesh_slice] = temporal_meshsort(emesh, idx, verbose)
%
% temporal_meshsort(emesh, st4idx_X50(strmatch('slp1', st4sym_X50)), 0)
% temporal_meshsort(emesh, st4idx_X50(strmatch('Alhambra', st4sym_X50)), 2);
% temporal_meshsort(emesh, st4idx_X50(strmatch('CG10479', st4sym_X50)), 2)

    % fixed parameters
    DIFF_THRESH = 5;
    % CLUS_DIST = 'euclidian';
    % correlation works better with experiments with different staining
    % intensities, otherwise euclidian is actually slightly better
    CLUS_DIST = 'correlation';
    % Threshold for comparing just 2 images with correlation distance
    % EVerything smaller than that is considered identical
    % Most APC cluster entries have differences less than 0.4 within
    THRESH_SAME = 0.4;
    num_tri = 0;
    
    % less than 3 embryos, APC doesn't work
    if length(idx) < 3,
        if length(idx) == 1,
            % one embryo -> just return it
            tidx = [idx, 1];
            tmesh = emesh.stain(:,idx);
            tmesh_slice = 1;
        else
            % 2 embryos
            
            % calculate the distance
            apc_data = emesh.stain(:, idx);
            apc_dist = squareform(pdist(apc_data', 'correlation'));
            
            if apc_dist < THRESH_SAME,
                % Embryos are the same
                tidx = [ idx(1), 1; idx(2), 1 ];
                tmesh = min(apc_data,[], 2);
                tmesh_slice = 1;
            else 
                % Embryos are presumably different
                % -> sort them naively
                stain_intensity = sum(abs(255 - emesh.stain(:, idx)), 1);
                if stain_intensity(1) < stain_intensity(2),
                    tidx = [ idx(1), 1; idx(2), 2 ];
                else
                    tidx = [ idx(1), 2; idx(2), 1 ];
                end
                tmesh = apc_data(:, tidx(:,2));
                tmesh_slice = [1,2];
            end
        end
        return;
    end
    

    N = length(idx);
    s = zeros(N*N-N,3); 
    
    apc_data = emesh.stain(:, idx);
    if strcmp(CLUS_DIST, 'euclidian'),
        apc_dist = squareform(pdist(normalize(apc_data')));
    else
        apc_dist = squareform(pdist(apc_data', CLUS_DIST));
    end
    
    j=1;
    for i=1:N
        for k=[1:i-1,i+1:N]
            s(j,1)=i; s(j,2)=k; s(j,3)=-apc_dist(i,k);
            j=j+1;
        end;
    end

    if verbose > 0, fprintf('Clustering ... '); end

    apc_p = median(s(:,3));
    apc_idx = apcluster(s, apc_p);
    
    apc_uidx = unique(apc_idx);
    apc_uidx_len = length(apc_uidx);

    if verbose > 0, 
        fprintf('found %d clusters\n', apc_uidx_len);
    end

    % Make sure there is more than one cluster
    if max(isnan(apc_idx)) || apc_uidx_len == 1,
        % if yes, just combine
        tidx = [ idx(1:end)', repmat(1, [length(idx), 1])];
        tmesh = min(apc_data,[], 2);
        tmesh_slice = 1;
        return;
    end

    apc_med = zeros(size(emesh.stain, 1), apc_uidx_len);
    for i=1:apc_uidx_len,
        % try to represent consensus either as median or min
        % apc_med(:, i) = median(apc_data(:, apc_idx == apc_uidx(i)),2);
        % NOTE: change also consensus for 2 embryos above!!!!
        apc_med(:, i) = min(apc_data(:, apc_idx == apc_uidx(i)),[], 2);
    end
    
    % Variants are:
    % - Euclidean distance
    % - normalized Eucledian distance
    % - correlation distance
    % UG = squareform(pdist(apc_med'));
    UG = squareform(pdist(normalize(apc_med')));
    % UG = squareform(pdist(apc_med', 'correlation'));
    UG = tril(UG);
    UG = sparse(UG);
    [ST,pred] = graphminspantree(UG);

    if verbose > 1,
        fprintf('Showing shortest path\n');
        view(biograph(ST,[],'ShowArrows','off','ShowWeights','on'));
    end
        
    
    tidx = zeros(length(idx), 2);
    tidx(:, 1) = idx;
    
    if num_tri == 1,
        % code for comparing/sorting # of stained triangeles
        
        if verbose > 0, fprintf('Discretizing meshes ... '); end
    
        med_dmesh = TImrf(emesh, apc_med, 'verbose', 0, 'nocorona');

        if verbose > 0, fprintf('done\nSorting them\n'); end


        med_sum = [sum(med_dmesh,1); apc_uidx'; 1:apc_uidx_len];
        med_ssum = sortrows(med_sum', 1);
        %sum(abs(255 - apc_med) .* med_dmesh) ./ sum(med_dmesh > 0)

        med_diff = diff(med_ssum(:,1)) ./ med_ssum(1:end-1,1) .* 100;

        med_da = find(med_diff < DIFF_THRESH);
        
        med_idx = med_ssum(:, 3);
    
        tidx = idx(med_ssum(:, 2));
    else
        % code for graph theoretic approach for sorting
        
        % find enpoints
        G_ep = sum(full(ST + ST') > 0) == 1; % <- shows endpoints in graph
        
        stain_intensity = sum(abs(255 - apc_med), 1);
        
        if verbose > 0, fprintf('Found %d endpoints in graph\n', sum(G_ep)); end

        if sum(G_ep) > 2,
            
            if sum(G_ep) > 3,
                fprintf('More than 3 endpoints found, cannot resolve\n');
            else
                
                fprintf('3 endpoints found\n');
                G_epidx = find(G_ep > 0);
                G_epp = zeros(3,2);
                G_epp(1,:) = [G_epidx(1), G_epidx(2)];
                G_epp(2,:) = [G_epidx(1), G_epidx(3)];
                G_epp(3,:) = [G_epidx(2), G_epidx(3)];
                
                path_len = zeros(3, 1);
                for i = 1:3,
                    [dist, path] = graphshortestpath(ST, G_epp(i,1), G_epp(i,2), 'directed',false);
                    path_len(i) = length(path);
                end
                
                [val, min_path] = min(path_len);
                
                G_a = setdiff(reshape(G_epp(1:3 ~= min_path ,:),4,1)', G_epp(min_path,:));

                if stain_intensity(G_a) < stain_intensity(G_epp(min_path,1)),
                    [val, G_b] = max([stain_intensity(G_epp(min_path,1)), stain_intensity(G_epp(min_path,2))],[], 2);
                    G_end = G_epp(min_path, G_b);
                    G_start = G_a;
                else
                    [val, G_b] = min([stain_intensity(G_epp(min_path,1)), stain_intensity(G_epp(min_path,2))],[], 2);
                    G_start = G_epp(min_path, G_b);
                    G_end = G_a;
                end
                
                G_out = G_epp(min_path, 1:2 ~= G_b);
                
                [dist, path] = graphshortestpath(ST, G_start, G_end, 'directed',false);
                [dist, path_out] = graphshortestpath(ST, G_out, G_end, 'directed',false);
                
                path_is = find(ismember(path,path_out) > 0);

                apc_staged = apc_uidx(path);

                for i=1:length(apc_staged),
                    tidx(apc_idx == apc_staged(i),2) = i;
                end
                
                tidx(apc_idx == apc_uidx(G_out), 2) = path(path_is(1));
                
                med_idx = [ path(1:path_is(1)), G_out, path(path_is(1)+1:end) ];
                tmesh_slice = [1:path_is(1), path_is(1), path_is(1)+1:length(path)];
            end
        else
            G_epnum = find(G_ep > 0);  
            [val, G_start] = min(stain_intensity(G_ep),[], 2);
            [val, G_end] = max(stain_intensity(G_ep), [], 2);
            G_start = G_epnum(G_start);
            G_end = G_epnum(G_end);
            
            % find the path through the graph
            [dist, path] = graphshortestpath(ST, G_start, G_end, 'directed',false);
            
            apc_staged = apc_uidx(path);
            
            for i=1:length(apc_staged),
                tidx(apc_idx == apc_staged(i),2) = i;
            end
           
            med_idx = path;
            tmesh_slice = 1:length(path);
        end
          
    end
    
    
    tmesh = apc_med(:, med_idx);

    if verbose > 0,
        
        figure;
        for i=1:apc_uidx_len,
            subplot(1, apc_uidx_len, i);
            displaymesh(apc_med(:, med_idx(i)), emesh.p_scale, emesh.t);
            ti_clus = sprintf('Cluster %d', med_idx(i));
            title(ti_clus);
        end
        
        if verbose > 1,
            figure;
            
            fig_pan = histc(apc_idx, sort(apc_uidx));
            max_pan = max(fig_pan);
            curr_pan = 1;
            
            for i=1:apc_uidx_len,
                cmesh = find(apc_idx == apc_uidx(i));
                for j=1:fig_pan(i),
                    subplot(apc_uidx_len, max_pan, curr_pan);
                    displaymesh(apc_data(:, cmesh(j)), emesh.p_scale, emesh.t);
                    curr_pan = curr_pan + 1;
                end
                curr_pan = curr_pan + (max_pan - fig_pan(i));
            end
            
        end
        
        
    end
    
end