function stain = meshstaining(im1, p_embryo, t, verbose, varargin)


method = 'cgreen';
persistent zt;

if nargin < 4
    verbose = 0;
else 
    % 2007/06: varargin parsing
    if ~isempty(varargin)
        i = 1;
        while i <= length(varargin),
            switch(varargin{i}),
                case 'method'
                    i = i+1;
                    method = varargin{i};
            end
            i = i+1;
        end
    end
end

% im1g = rgb2gray(im1); never used

switch method
    case 'gray'
        if ~isgray(im1),
            im1s = rgb2gray(im1);
        else
            im1s = im1;
        end
    case 'cgreen'
        im1s = im1(:,:,2); % staining only
    case 'cred'
        im1s = im1(:,:,3);
    case 'corrdic'
        if isempty(zt),
            dthresh = 0.005;
            
            x = repmat(0:255, [256 1]);
            y = x';
            z = (x .* 256) ./ (y + 1);
            zg = del2(z);
            zt = zg < dthresh;            
        end
        imrgb = double(im1);        
        [m, n, r] = size(imrgb);
    
        im_r = im2col(imrgb(:,:,1),[1 1], 'distinct');
        im_b = im2col(imrgb(:,:,3),[1 1], 'distinct');

        % Fixed 2009/10 to account for rare cases (virtually none in BDGP set) were im_r=0
        % Code used for TIs was: im_bmin = 255 - sum(zt(:,im_r));
        im_bmin = 255 - sum(zt(:,im_r+1));
        im_btadj = im_b - im_bmin;
        im_badj = im_btadj <= 0;
        im_b(im_badj) = im_bmin(im_badj);
        
        im1s = (im_r .* 256) ./ (im_b + 1);
        im1s = col2im(uint8(im1s), [1 1], [m n]);
    otherwise
        error('meshstaining: Method %s not known', method);
end

p_embryo_f = fix(p_embryo);

Ntri = size(t,1);

stain = zeros(Ntri,1);

if verbose > 0
    fprintf('Starting to process %d triangles:\n', Ntri);
    hold on;
    colormap(gray);
    axis ij;
end

for i = 1:Ntri
    if verbose > 0
        fprintf('Processing triangle %d = ', i);
    end
    
    roic = p_embryo_f(t(i,:),1)';
    roir = p_embryo_f(t(i,:),2)';
    
    % was: im_roi = roipoly(im1s, roic, roir);
    im_roi = roipoly(im1s, roic, roir);
    stain(i) = median(im1s(im_roi == 1));
    
    if verbose > 0
        fprintf('%d\n', stain(i));
        fill(roic, roir, stain(i));
    end
end

if verbose > 0
    hold off;
end;
