function [elabel] = exp_label2(emesh, fn, thresh, varargin)
% [elabel] = exp_label(emesh, fn, dmesh, dmesh_idx)
% similar to matlab bwlabel: find connnected regions with expression
% uses the discretized meshes for determining the expression probabilities
% returns linear array with numbered labels (like bwlabel)
%
% fn = filename or index (in emesh) of the embryo
% thresh = set threshold, default 1 if empty
% variable args:
% 'dmesh', dmesh, dmesh_idx
%       dmesh = discretized mesh
%       dmesh_idx = discretized mesh indices
% 'min_size', size
%       size = minimum number of triangles to be labeled


dmesh = [];
dmesh_idx = [];
min_size = 0;


if ~isempty(varargin),
    var=1;
    while var <= length(varargin),
        switch varargin{var}
            case 'dmesh'
                dmesh = varargin{var+1};
                dmesh_idx = varargin{var+2};
                var = var+2;
            case 'min_size'
                var = var+1;
                min_size = varargin{var};
            otherwise
                error('Unknown argument: %s', varargin{var});
        end
        var = var + 1;
    end
end




if isempty(thresh),
    thresh = 1;
end

% datapoints
dp = size(emesh.stain,1);

if ischar(fn),
    eidx = strmatch(fn, emesh.files);
else
    eidx = fn;
end


ed = [];

if isempty(dmesh),
    try,
        ed = emesh.dstain(:, eidx);
    catch
        error('emesh.dstain not found, specify it in command line');
    end
else
    ed = dmesh(:, find(dmesh_idx == eidx));
end

% threshold them
et = find(ed >= thresh);
et_all = et;

elabel = zeros(dp,1);

while ~isempty(et)
  
  dig_conn = et(1);
  dug_conn = [];
  cconn = [];
  while ~isempty(dig_conn)
      conn = [dig_conn(1), emesh.tn{dig_conn(1)}];
      cconn = unique([ cconn, intersect(conn, et_all) ]);
      dug_conn = [dug_conn, dig_conn];
      dig_conn = setdiff(cconn, dug_conn);
  end
      
  %conn = [et(1), emesh.tn{et(1)}];
  %cconn = intersect(conn, et_all); % here some pruning would be cool
  
  clab = intersect(cconn, find(elabel > 0));
%  clab = unique(elabel(cconn));
  if isempty(clab),
      % no previous labels found
      elabel(cconn) = max(unique(elabel)) + 1;
  else
      % previous labels found
      clab_u = unique(elabel(cconn));
      elabel(cconn) = max(clab_u);
      if length(clab_u) > 2,
        % 2-several previous labels found -> make all the same
        for i=2:length(clab_u),
          elabel(elabel == clab_u(i)) = max(clab_u);
        end
      end
  end
  et = setdiff(et, cconn);
end

% prune for min_size
if min_size > 0,
    labels = unique(elabel);
    for i=1:length(labels),
        if length(find(elabel == labels(i))) < min_size,
            elabel(elabel == labels(i)) = 0;
        end
    end
    
end

% renumber labels in consecutive numbers
labels = unique(elabel);
for i=2:length(labels),
    elabel(elabel == labels(i)) = i-1;
end