function [w, orig, im, bim, S2] = segment_antirrhinum_anthocyanin_leaves(I, w, binim, segsectors)

% function [w, orig, im, bim, S2] = antirrhinum_anthocyanin_leaves(I, w)

% 

% A function that takes an image and a handle to a waitbar and returns the

% waitbar handle, the original image, the image with the sectors overlayed,

% the binary image of the organ, and the binary image of the sectors.

%

% This script has been written for Sandra Bensmihen (JIC) for use with the 

% expression of a cell autonomous marker (anthocyanin) in Antirrhinum

% majus leaves.

%

% Dr. A. I. Hanna (2006)


if nargin<1

    error('Must at least give an image');

end

if nargin<2

    w = waitbar(0, 'Processing image...');

end



% Make sure that the image is in a descend format.

I = double(I);

I = I./max(I(:));

orig = I;



% Calculate the binary image

bim = calcBinIm(orig);
S2 = [];
im = [];
if segsectors==0
    return;
end

waitbar(1/10, w);

% Label the binary image

L = bwlabel(bim);

stats = regionprops(L, 'BoundingBox');

BB = round(stats.BoundingBox);



% Extract the leaf from the full RGB image

for i=1:3

    I2(:,:,i) = I(:,:,i).*bim;

    I3(:,:,i) = imcrop(I2(:,:,i), BB);

end

waitbar(2/10, w);

I2 = I3;

clear I3;

H = fspecial('gaussian', 5, 5);

A = normalise(I2(:,:,1) - I2(:,:,2).^5);

A2 = normalise(I2(:,:,1) - I2(:,:,2));

thehist = load('leafhist2');

thehist = thehist.h;

waitbar(3/10, w);

A2 = histeq(A2, thehist);

A = conv2(A, H, 'same');

A2 = conv2(A2, H, 'same');

waitbar(4/10, w);

h = imhist(A);

for i=length(h):-1:1

    if (h(i) >= 300)

        thr = i/256;

        break

    end

end

waitbar(5/10, w);

if(i==1)

    thr = 0.8;

end

S = A>thr;

p = S.*A;

q = imadjust(p, [.5 1], [], 1);

S = (q> .4) + (A> .8) + (A2>.6);

waitbar(6/10, w);

v = normalise(A.*A2);

h = imhist(v);

lh = log(h)./max(log(h));

for i=length(lh):-1:1

    if(lh(i) >= 0.46)

        thr = i/256;

        break

    end

end

waitbar(7/10, w);

v = A2.^3;

h = imhist(v);

h = h./max(h);

dh = diff(h);

for i=length(dh):-1:1

    if(dh(i) < -1.403e-3)

        thr = i/256;

        break;

    end

end

%S = (S>0).* (normalise(A.*A2)>thr);

S = (normalise(A.*A2)>thr);

waitbar(8/10, w);

S = normalise((histeq(A2).*histeq(A)))>.9;

v = normalise(A.*A2);

S = A2.^3 > thr;

S = (A2.^3>.5) & (imextendedmax(v, .2));

S = (A2>.7) & (I2(:,:,2) < .7);

S = imclose(S, strel('disk', 5));

waitbar(9/10, w);

S = bwfill(S, 'holes');

S2 = zeros(size(orig(:,:,1)));

S2(BB(2):BB(2) + BB(4), BB(1):BB(1)+BB(3)) = S;

waitbar(10/10, w);

im = I;

return;

%%%%%

%

%

%%%%%

function n = normalise(a)

if ndims(a)==3           % Assume colour image

    hsv = rgb2hsv(a);

    v = hsv(:,:,3);

    v = v - min(min(v)); % Just normalise value component

    v = v/max(max(v));

    hsv(:,:,3) = v;

    n = hsv2rgb(hsv);

else                     % Assume 1D or 2D greyscale values

    if ~isa(a,'double'), a = double(a); end

    n = a - min(min(a));

    n = n/max(max(n));

end

%%%%%

%

%

%%%%%

function bim = calcBinIm(orig)

d = orig(:,:,1) -orig(:,:,3);

[counts, X] = imhist(d);

ind = intersect(find(X>0), find(X<.3));

X = X(ind);

counts = counts(ind);

ind = find(counts>1000);

if isempty(ind)

    thr = 0.1;

else

    thr = X(ind(1));

end
bim = d> thr;
bim = imclearborder(bim);
bim = bwareaopen(bim, 50000);
bim = imclose(bim, strel('disk', 10));
bim = bwfill(logical(bim), 'holes');