function [similar score dscore] = TIbin2(emesh, eg_sel, dom_labels, top_no, varargin)
% Cleaned up version of find_restricted
% Finds meshes in eg_sel with shared boundaries to each label (~=0) in dom_labels
% [similar score dscore] = TIbin2(emesh, stain_lib, dom_labels, top_no, ...)
%
% emesh: mesh storage and parameters (t,tn)
% eg_sel: the collection of meshes that will be searched
% dom_labels: labelled domains for each datapoint (0 is ignored)
% top_no: Number of results to be returned
%
% returns a cell array "similar" with top_no candidates, "score" with
% top scores, and, if specified dscore with all the detailed scores
% 
% dscore is a cell array containing 3D arrays (staining probs, in/out,
% embryo_no) if 'details' option is given
%
% Variable arguments recognized:
% 'method', '<method>': <method> is mean, local, localmax or
%                       naivemrf (default)
%                       The default naivemrf corresponds to the manuscript
% 'scoring', '[descending|ascending]': sorts scores ascending or descending
% 'details': fills dscore variable
% 'correlation':    'positive' same as pos_correlation
%                   'negative' opposite as positive
%                   'both' default
% 'full': Evaluate fully connected triangles (2 shared corner points) only
%           (not default but recommended)
% For naivemrf only:
% 'mrf_smooth': The smoothness factor for scoring inside domain
% 'mrf_border': The difference that qualifies as edge
% both default to 10
% 'mrf_normalize': Normalize the MRF scores to number of triangles
% (default on)

    % defaults
    method = 'naivemrf';
    scoring = -1;
    detail_score = 0;
    pat_corr = 0;
    mrf_smooth = 10; mrf_border = 10; mrf_normalize =1;
    t_conn = emesh.t_conn;

     
     
    if ~isempty(varargin),
        var=1;
        while var <= length(varargin),
            switch varargin{var}
                case 'method'
                    var = var + 1;
                    choices = {'mean', 'local', 'localmax','naivemrf'};
                    if ~strmatch(varargin{var}, choices),
                        error('Method %s unknown', varargin{var});
                    else
                        % fprintf('Method: %s\n', varargin{var});
                        method = varargin{var};
                    end
                case 'scoring'
                    var = var + 1;
                    if strmatch(varargin{var}, 'ascending'),
                        scoring = 1;
                    else
                        scoring = -1;
                    end
                case 'full'
                    t_conn = fully_connected(emesh);
                case 'details'
                    detail_score = 1;
                case 'correlation'
                    var = var+1;
                    if strmatch(varargin{var}, 'positive'),
                        pat_corr = 1;
                    elseif strmatch(varargin{var}, 'negative'),
                        pat_corr = -1;
                    elseif strmatch(varargin{var}, 'both'),
                        pat_corr = 0;
                    else
                        error('Correlation parameter %s not known', varargin{var});
                    end
                case 'mrf_smooth'
                    var = var + 1;
                    mrf_smooth = varargin{var};
                case 'mrf_border'
                    var = var + 1;
                    mrf_border = varargin{var};
                otherwise
                    error('Unknown argument: %s', varargin{var});
            end
            var = var + 1;
        end
    end
    
  
    dom = unique(dom_labels);
    if ~isempty(find(dom == 0)),
        dom = dom(find(dom ~= 0));
    end
    dom_no = length(dom);

    fprintf('Finding similar patterns to domains (%s): ', method);
    for cdom=1:dom_no
        fprintf('%d, ', cdom);
        
        cdm = logical(zeros(311,2));
        % cdm(:,1) contains the triangles inside the domain
        cdm(find(dom_labels == dom(cdom)), 1) = 1;
        [c_t, c_v] = find(t_conn(cdm(:,1),:) == 1);
        % cdm(:,2) contains the corona triangles
        cdm(unique(c_v), 2) = 1;
        cdm(cdm(:,1), 2) = 0; % eliminate the ones inside domain
        cdm_dscore = []; % just for safety
        
        switch method,

            case 'mean'
                cdm_inside = mean(eg_sel(cdm(:,1), :),1) .* (1 - std(eg_sel(cdm(:,1), :),1));
                cdm_corona = mean(eg_sel(cdm(:,2), :),1) .* (1 - std(eg_sel(cdm(:,2), :),1));
                cdm_diff = cdm_inside - cdm_corona;

                cdm_diff = [cdm_diff; 1:length(cdm_diff)]';
                cdm_diff = sortrows(cdm_diff, scoring); % sort descending

            case 'dist'
                cdm_dm = eg_sel(cdm(:,2),:)';
                cdm_inside = mean(eg_sel(cdm(:,1), :), 1);
                cdm_dm = [ repmat(cdm_inside', [1 size(cdm_dm, 2)]); cdm_dm];
                cdm_dist = squareform(pdist(cdm_dm, 'seuclidean'));
                cdm_dist = cdm_dist(length(cdm_dist)/2+1:end, 1:length(cdm_dist)/2);
                cdm_diff = diag(cdm_dist);

                cdm_diff = [cdm_diff'; 1:length(cdm_diff)]';
                cdm_diff = sortrows(cdm_diff, scoring); % sort ascending
            
            case 'roieuclidean'    
                roi = cdm(:,1) + cdm(:,2);
                roi = logical(roi);
                equery = zeros(length(cdm(:,1)),1);
                equery(cdm(:,1)) = 1;
                equery = equery(roi);
                estain = eg_sel(roi,:);
                
                n = size(eg_sel,2);
                edist = sum(abs(estain'  - repmat(equery', [n 1])) .^ 2, 2);
                cdm_diff = edist;
                
                cdm_diff = [cdm_diff'; 1:length(cdm_diff)]';
                cdm_diff = sortrows(cdm_diff, -scoring); % sort ascending
                cdm_score = [cdm_diff(:,1)];

            case 'mebyme'
                roi = cdm(:,1) + cdm(:,2);
                roi = logical(roi);
                equery = zeros(length(cdm(:,1)),1);
                equery(cdm(:,1)) = 1;
                equery = equery(roi);
                estain = eg_sel(roi,:);

                dp = length(equery);
                n = size(eg_sel,2);
                ex = (estain - repmat(mean(estain),[dp,1])) ./ repmat(std(estain),[dp,1]);
                ey = (equery - repmat(mean(equery),[dp,1])) ./ repmat(std(equery),[dp,1]);
                edist = 1 - abs(sum(ex .* repmat(ey,[1,n])) ./ dp);
                cdm_diff = edist';

                cdm_diff = [cdm_diff'; 1:length(cdm_diff)]';
                cdm_diff = sortrows(cdm_diff, scoring); % sort ascending
                cdm_score = [cdm_diff(:,1)];

            case 'local'
                avgfx = @(cdm_dm, lt_dom) sum(cdm_dm, 2) ./ sum(lt_dom, 2);
                [cdm_diff cdm_score cdm_dscore] = localdiff(avgfx);
                % TODO: localmin (min of neighboring triangles)
                %       localone (only one neighboring triangle)

            case 'localmax'
                % identical to local - uses max instead of average
                maxfx = @(cdm_dm, lt_dom) max(cdm_dm, [], 2);
                [cdm_diff cdm_score cdm_dscore] = localdiff(maxfx);
                
            case 'overlap'
                % take sum/product of all stain inside domain
                % ? uniformity -> std, t-test (narrow should be first)
                cdm_inside = mean(eg_sel(cdm(:,1), :),1) ./ std(eg_sel(cdm(:,1), :),1);
                cdm_diff = cdm_inside;
                
                cdm_diff = [cdm_diff; 1:length(cdm_diff)]';
                cdm_diff = sortrows(cdm_diff, scoring); % sort descending

                cdm_score = [cdm_diff(:,1)];

            case 'naivemrf'
                % Very naive MRF from MRF book. The smaller the number the
                % better
                %[cdm_diff cdm_score] = naive_mrf(10, 10);
                [cdm_diff cdm_score] = naive_mrf(mrf_smooth, mrf_border);
                
        end
        
        similar{cdom} = selection(cdm_diff(1:top_no,2));
        score{cdom} = cdm_score(1:top_no, :);
        if detail_score == 1 && ~isempty(cdm_dscore),
            dscore{cdom} = cdm_dscore(:,:,cdm_diff(1:top_no,2));
        else
            dscore = [];
        end
        
    end
    fprintf('done\n');
    
    
    function [cdm_diff cdm_score cdm_dscore] = localdiff(diff_fx)
        % connectivity matrix
        lt_conn = t_conn;
        % set all rows which are not neighboring to inside triangles to 0
        lt_conn(:, ~cdm(:,2)) = 0;
        % only look at rows with triangles inside domain for
        % connectivity to the outside
        lt_dom = lt_conn(cdm(:,1),:);
        eg_dom = eg_sel(cdm(:,1),:); 
        % 11/13: eliminate the triangles deep inside the domain
        % that don't have triangles on the outside to eliminate
        % divisions by 0
        lt_dom = lt_dom(sum(lt_dom,2) > 0,:);
        eg_dom = eg_dom(sum(lt_dom,2) > 0,:);
        % initialize
        cdm_no = size(eg_sel,2);
        cdm_diff = zeros(cdm_no, 1);
        % foreach embryo
        cdm_dscore = [];
%        cdm_h = zeros(cdm_no, 1); cdm_sig = zeros(cdm_no, 1);
        for embryo = 1:cdm_no,
            eg_mult = repmat(eg_sel(:, embryo)', [size(lt_dom, 1) 1]);
            cdm_dm = lt_dom .* eg_mult;
            cdm_dm = diff_fx(cdm_dm, lt_dom);
            cdm_diff(embryo) = sum(eg_dom(:, embryo)  - cdm_dm);
            % playing around with returning all the details
            if detail_score == 1,
                cdm_dscore(:,:,embryo) = [eg_dom(:,embryo), sum(cdm_dm,2)];
            end
            % playing around with t-test inside to outside - doesn't
            % discriminate
%           [cdm_h(embryo), cdm_sig(embryo)] = ttest2(eg_sel(cdm(:,1),embryo), eg_sel(cdm(:,2), embryo));
        end

        cdm_diff = [cdm_diff'; 1:length(cdm_diff)]';
        cdm_diff = sortrows(cdm_diff, scoring); % sort ascending

        cdm_mean_ins = mean(eg_sel(cdm(:,1), cdm_diff(:,2)),1); 
        cdm_mean_out = mean(eg_sel(cdm(:,2), cdm_diff(:,2)),1);
        cdm_score = [cdm_diff(:,1),...
            repmat(length(find(cdm(:,1) == 1)), [cdm_no 1]),...
            repmat(size(lt_dom, 1), [cdm_no 1]),...
            repmat(length(find(cdm(:,2) == 1)), [cdm_no 1]),...
            cdm_mean_ins', cdm_mean_out',...
%            cdm_h(cdm_diff(:,2)), cdm_sig(cdm_diff(:,2)), ...
            ];
        % score: 
        % sum score, # inside domain, # inside with outside
        % triangles, # outside domain touching domain
        % mean inside domain, mean outside domain
      
    end

    function [cdm_diff cdm_score] = naive_mrf(thresh_smooth, thresh_border)
        % create connectivity matrix with inside neighbors
        cdm_i = repmat(cdm(:,1), [1,311]);
        t_conn_i = t_conn;
        t_conn_i(:, ~cdm(:,1)) = 0;
        t_conn_i = t_conn_i .* ~eye(311);
        cdm_i = cdm_i .* t_conn_i;
        cdm_i_no = sum(sum(cdm_i,1));
        cdm_i_lin = logical(reshape(cdm_i, 1, []));

        % create connectivity matrix with outside neighbors
        cdm_c = repmat(cdm(:,1), [1,311]);
        t_conn_c = t_conn;
        t_conn_c(:, ~cdm(:,2)) = 0;
        t_conn_c = t_conn_c .* ~eye(311);
        cdm_c = cdm_c .* t_conn_c;
        cdm_c_no = sum(sum(cdm_c,1));

        % initialize
        cdm_no = size(eg_sel,2);
        cdm_diff = zeros(cdm_no, 1);
        % foreach embryo
        cdm_dscore = [];
        

        
        for embryo = 1:cdm_no,
            % matrix of each embryo repeated
            % ... horizontally
            eg_mc = repmat(eg_sel(:, embryo), [1,311]);
            % ... vertically
            eg_mr = repmat(eg_sel(:, embryo)', [311,1]);
            
            % use repeated eg_mc/eg_mr matrices to calculate
            % differences between neighboring triangles
            % -> multiply with connectivity matrix
            % ... inside domain
            ki = eg_mc .* cdm_i - eg_mr .* cdm_i;
            % ... corona
            kc = eg_mc .* cdm_c - eg_mr .* cdm_c;
            
            if thresh_smooth > 0,
                % find out how many inside triangle pairs have
                % differences smaller than thresh_smooth
                ki = abs(ki(cdm_i_lin)) < thresh_smooth;
                % find out how many corona triangle pairs have
                % differences larger than thresh_border
                if pat_corr > 0,
                    kc = -kc - thresh_border;
                elseif pat_corr < 0
                    kc = kc - thresh_border;
                else
                    kc = abs(kc) - thresh_border;
                end
                kc = kc > 0;
                
                if mrf_normalize,
                    % try to normalize for # of triangles to make
                    % scores comparable (best = 2 (all inside/corona
                    % triangles below threshold), worst = 0)
                    cdm_diff(embryo) = sum(sum(ki,2) / cdm_i_no * sum(kc,2) / cdm_c_no);
                    
                    % multiplication makes a bit more sense?
                    % cdm_diff(embryo) = sum(ki) / cdm_i_no * sum(kc) / cdm_c_no;
                else
                    cdm_diff(embryo) = sum(ki) + sum(kc); % TODO
                end
            else
                ki = ki * ki;
                kc = kc * kc;
                cdm_diff(embryo) = sum(sum(ki,1)) + sum(sum(kc,1));
            end
        end

        
        cdm_diff = [cdm_diff'; 1:length(cdm_diff)]';
        cdm_diff = sortrows(cdm_diff, scoring); % sort ascending

        cdm_score = [cdm_diff(:,1)];
    end
    
end



function [t_conn] = fully_connected(emesh)

    for i=1:length(emesh.t)
        tnfr{i} = [];   
        tnf{i} = []; 
        for j=1:size(emesh.t,2)       
            [prow, pcol] = find(emesh.t==emesh.t(i,j));
            tnfr{i} = [tnfr{i}, prow(prow ~= i)'];
        end
        for j=1:length(emesh.tn{i})
           if length(find(tnfr{i} == emesh.tn{i}(j))) > 1,
            tnf{i} = [tnf{i}, emesh.tn{i}(j)];
          end
        end
    end

    t_conn=zeros(311,311);
    for i=1:length(tnf), 
        for j=1:length(tnf{i}), 
            t_conn(i,tnf{i}(j)) = 1; 
            t_conn(i,i) = 1; 
        end;
    end

end