function [results, rdist] = esimilar3(emesh, idx, query, varargin)
% [results] = esimilar3(emesh, idx, query, varargin)
% e.g.
% esimilar3(emesh, allStage4_idx, 'insitu20707.jpe')
%
% emesh: emesh struct
% idx: use only data with those indices
% query: query, either filename or TI
% 
% Additional parameters:
% images: Path for image directories
% range: Display results in that range (e.g. [10:20])
% db: db connection to MySQL database for symbol queries
% distance: distance function for Matlab pdist or 
%           'no_pdist' for internal correlation distance calculation (may
%           be faster and less memory consuming). 'no_pdist' doesn't
%           throw an error if unstained/ubiquitous patterns are included. 
% verbose:  0=no display
%           1=text results
%           2=images (default)
%

% persistent edist
%if isempty(edist),
%    edist = squareform(pdist(emesh.stain(:,idx)', 'correlation'));
%end

fpath = [];
range = [];
conn = [];
dist_fx = 'correlation';
verbose = 2;

if ~isempty(varargin)
    i = 1;
    while i <= length(varargin),
        switch(varargin{i}),
            case 'images'
                i=i+1;
                fpath = varargin{i};
            case 'range'
                i=i+1;
                range = varargin{i};
            case 'db'
                i=i+1;
                conn = varargin{i};
            case 'distance'
                i=i+1;
                dist_fx = varargin{i};
            case 'verbose'
                verbose = 1;
            otherwise
                error('Argument %s not known', varargin{i});
        end
        i = i+1;
    end
end


if isempty(range)
    range = 1:10;
end

estain = emesh.stain(:,idx);

filename = '';
if ischar(query)
    filename = query;
    
    ef_idx = strmatch(filename, emesh.files);
    
    if isempty(ef_idx),
        error('Error: %s doesn''t have a mesh associated in the database', filename);
    end
    
    qidx = find(idx == ef_idx);
    equery = estain(:, qidx);
    
elseif isvector(query)
    equery = double(query);
    estain = [equery, estain];
    qidx = 1;
else
    error('Query must be a filename or vector of staining intensities');
end

if strcmp(dist_fx, 'no_pdist'),
    dp = size(estain,1);
    n = size(estain,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);
    edist = edist';
    [val, esim] = sortrows(edist);
else
    edist = squareform(pdist(estain', dist_fx));
    [val, esim] = sortrows(edist(:,qidx));
end

results = idx(esim(min(range)+1:max(range)+1));
rdist = val(min(range)+1:max(range)+1);
[rho, pval] = corr(estain(:, esim(min(range):max(range)+1)));
pval = pval(1, 2:max(range)+1);

if verbose > 1,
    mpanels = (max(range) - min(range) + 2);
    npanels = 2;

    figure;
    if qidx ~= 1.
    subplot(mpanels,npanels,1);
    displayimage(conn, fpath, filename);
    end
    subplot(mpanels,npanels,2);
    displaymesh(equery, emesh.p_scale, emesh.t);
    for i = range
        enum = esim(i+1);
        subplot(mpanels, npanels, (i - min(range) + 1)*2 +1);
        displayimage(conn, fpath, emesh.files{idx(enum)});
        subplot(mpanels, npanels, (i - min(range) + 1)*2 +2);
        displaymesh(estain(:,enum), emesh.p_scale, emesh.t);
    end
end

if verbose > 0,
    fprintf('%s\n', emesh.files{results});
    fprintf('%e\n', pval);
end

end

function displayimage(conn, fpath, filename)
    if isempty(fpath)
        return
    end
    
    sym=[];
    fullname = myinsitu2fn(filename);
    if ~isempty(conn)
        sym = pg_query(conn, 0, my_insitus('symbol'), fullname);
    end
        
    fullname_path = strcat(fpath,'/',fullname);

    im = imread(fullname_path);
    im2 = imresize(im, 0.5);
    imshow(im2);
    axis off;
    if isempty(sym),
        suptitle(filename);
    else
        title([sym{1}, ' (', filename, ')']);
    end
end