function [w, orig, im, binim, snew] = amelia_petal_sectors(im, w, binim)
%#function checknargin
%clear all;
%w = waitbar(0, 'Testing...');
%im = imread('C:\Documents and Settings\aih\My Documents\post_doc\amelia\Images\StagedTestSet\S4\P10F15P1I.JPG');
%im = imread('C:\Documents and Settings\aih\My Documents\post_doc\amelia\Images\StagedTestSet\S4\P7F22P1I.JPG');
%im = imread('C:\Documents and Settings\aih\My Documents\post_doc\amelia\Images\StagedTestSet\S1\P12F8P1I.JPG');
%im = imread('C:\Documents and Settings\aih\My Documents\post_doc\amelia\Images\StagedTestSet\S4\P12F18P1I.JPG');
%im = imread('C:\Documents and Settings\aih\My Documents\post_doc\amelia\Images\StagedTestSet\S4\P10F17P1I.JPG');
%im = imread('C:\Documents and Settings\aih\My Documents\post_doc\amelia\Images\StagedTestSet\S3\P12F12P1I.JPG');
%im = imread('C:\Documents and Settings\aih\My Documents\post_doc\amelia\Images\P7\P21F11-2P1I_small_cluster.JPG');
%im = imread('C:\Documents and Settings\aih\My Documents\post_doc\amelia\Images\P7\F33F3082_cells.JPG');
%im = imresize(im, .5);
if nargin<2
    w = waitbar(0, 'Segmenting...');
    sz = size(im);
end
im = double(im);
im = im./max(im(:));

r = im(:,:,1);
[h, c] = imhist(r);
ind = intersect(find(c>.5), find(c<.8));
h = h(ind); c = c(ind);
[val, ind] = min(h);
thr = c(ind);

%im = uint8(im.*255);
%im = double(im);
%im = im./max(im(:));
orig = im;
im = imadjust(im, stretchlim(im));
if (nargin<3)
    binim = calc_binim(im, thr);
end
if(sum(binim(:)) == 0)
    fprintf('Error with the binary image\n');
    orig = zeros(size(im));
    snew = zeros(size(im));
    return;
end
waitbar(1/50, w);
imn = rgb2ntsc(im);
a = imn(:,:,2);
a = normim(a);
b = imn(:,:,3);
b = normim(b);
k = a.*b;
k = normim(k);
bwim = k > .3;
waitbar(2/50, w);
bwim = bwareaopen(bwim, 30);
y = k.*binim;
y = double(y);
y = y./max(y(:));
y2 = imadjust(y, [.4 1], [], .4);
h = imhist(y);
waitbar(3/50, w);
f = filter(ones(10, 1)/10, [1], h);
d = rot90(filter(ones(10, 1)/10, [1], diff(f)), 2);
d = d./max(d(:));
for i=2:length(d)-1
    if(d(i) < -0.1)
        thr = (255 - i)/255;
        break;
    end
end
s =(y>thr);
s = (s>0);
s = s.*(y2>.3);
L = bwlabel(s);
stats = regionprops(L, 'BoundingBox', 'Area');
snew = zeros(size(s));
for i=1:max(L(:))
    waitbar((i+4)/(2*length(stats)), w);
    u = (L==i);
    bbox = stats(i).BoundingBox;
    bbox(1:2) = floor(bbox(1:2));
    ind = find(bbox<1); bbox(ind) = 1;
    
    bbox(3:4) = ceil(bbox(3:4)) + 1;
    su = u(bbox(2):bbox(2)+bbox(4), bbox(1):bbox(1)+bbox(3));
    sy = y(bbox(2):bbox(2)+bbox(4), bbox(1):bbox(1)+bbox(3));
    k = (imadjust(sy.*su, [0 .8],[]));
    sn = (imadjust(k, stretchlim(k))>.7);
    sn = bwfill(sn, 'holes');
    snew(bbox(2):bbox(2)+bbox(4), bbox(1):bbox(1)+bbox(3)) = snew(bbox(2):bbox(2)+bbox(4), bbox(1):bbox(1)+bbox(3)) + sn;
end
s = snew;
s = bwfill(s, 'holes');
f = fspecial('gaussian', 2, 5);
L = bwlabel(s);
stats = regionprops(L, 'BoundingBox', 'Area');
if(mean([stats.Area])>1000)
    uthr = .35;
elseif(mean([stats.Area])>500)
    uthr = .2;
else
    uthr = .1;
end
uthr_orig = uthr;
snew = zeros(size(s));
waitbar(4/50, w);
G_CHANNEL = orig(:,:,2);
for i=1:length(stats)
    waitbar((2*(i+4))/(2*length(stats)), w);
    u = (L==i);
    bbox = stats(i).BoundingBox;
    bbox(1:2) = floor(bbox(1:2));
    ind = find(bbox<1); bbox(ind) = 1;
    
    bbox(3:4) = ceil(bbox(3:4)) + 1;
    su = u(bbox(2):bbox(2)+bbox(4), bbox(1):bbox(1)+bbox(3));
    sy = y(bbox(2):bbox(2)+bbox(4), bbox(1):bbox(1)+bbox(3));
    eq = su.*sy;
    eq = imadjust(eq, stretchlim(eq));
    eqa = histeq(eq);
    eqa = imadjust(eqa, stretchlim(eqa));
    if(((size(su,1)*size(su,2))>1000) & (uthr < .3))
        uthr = .2;
    else
        uthr = uthr_orig;
    end
    if((size(su,1)*size(su,2)) <= 30000)
        sn = imextendedmax(eqa, .3) +  imextendedmax(eqa, .2) + imextendedmax(eq, uthr);
    else
        sn = su.*(y2(bbox(2):bbox(2)+bbox(4), bbox(1):bbox(1)+bbox(3)) > .3);
    end
    sn = sn>0;
    if (sum(sn(:)) == 0)
        sn = eq>0;
    end
    sn = bwfill(sn, 'holes');
    sn = sn + (G_CHANNEL(bbox(2):bbox(2)+bbox(4), bbox(1):bbox(1)+bbox(3))<.2);
    sn = sn>0;
    snew(bbox(2):bbox(2)+bbox(4), bbox(1):bbox(1)+bbox(3)) = snew(bbox(2):bbox(2)+bbox(4), bbox(1):bbox(1)+bbox(3)) + sn;
end
%snew = snew + (orig(:,:,2) < .2);
%snew = (snew > 0);
snew = snew.*binim;
snew = bwfill(snew, 'holes');
L = bwlabel(snew);
stats = regionprops(L);
if (mean([stats.Area]) > 400)
    snew = bwareaopen(snew, 10);
end
im = add_spots(orig, snew, [0 0 1]);

function [binim] = calc_binim(im, thr);
%p = imadjust(im(:,:,1), stretchlim(im(:,:,1)));
%binim = normim(p) > .4;
p = imadjust(im(:,:,1), stretchlim(im(:,:,1))) - imadjust(im(:,:,3), stretchlim(im(:,:,3)));
binim = im(:,:,1)>thr;
binim = binim | normim(p) > .45;
binim = bwareaopen(binim, 5000);
binim = imclearborder(binim);
%L = bwlabel(binim);
%S = regionprops(L,'Area');
%[val, ind] = max([S.Area]);
%ind = find(L==ind);
%binim = zeros(size(binim));
%binim(ind) = 1;

dim = im(:,:,1) - im(:,:,3);

[counts, X] = imhist(dim);
ind = intersect(find(X>.2), find(X<.6));
counts = counts(ind); X = X(ind);
[val, ind] = min(counts);
thr = X(ind);

binim = binim | dim>thr;
binim = bwareaopen(binim, 5000);
binim = imclearborder(binim);


binim = imclose(binim, strel('disk', 40));
binim = bwfill(binim,'holes');
binim = imdilate(binim, strel('disk', 2));
return;

function [im] = add_spots(im, spots, col)

if (nargin == 2)
    col = [1 0 0];
end
ind = find(spots == 1);
imr = im(:,:,1);
img = im(:,:,2);
imb = im(:,:,3);
imr(ind)=col(1);
img(ind)=col(2);
imb(ind)=col(3);
im(:,:,1) = imr;
im(:,:,2) = img;
im(:,:,3) = imb;