%% Watershed on tougher images
clear
close all

img = im2double(imread('QD_AFM.png'));

figure
imshow(img)

% Convert to gray image, filter, dilate.
img2 = rgb2gray(img);
img3 = mat2gray(imtophat(img2,strel('disk',5)));
img4 = imdilate(img3,strel('square',2));

figure
imshow([img2 img3 img4],[])

% Set a threshold to identify all the QDs
bw = im2bw(img4,graythresh(img3));

figure
imshowpair(img2,bw)

% Use the bwdist to calculate the distance from the background
distImg = mat2gray(bwdist(~bw));

% Then we can take the complement image, and the centers of the circle are
% now the lowest points.  They will fill up with the watershed.
distImg2 = imcomplement(distImg);
figure
imshow([distImg distImg2])

% Now we need to set the background pixels back to zero, and we can run the
% watershed
distImg3 = distImg2.*bw;
figure
imshow(distImg3)

% Now we can call the watershed
L = watershed(distImg3);
labelImg = label2rgb(L, @spring, 'c', 'shuffle');

figure
subplot(1,2,1)
imshow(labelImg,[])
subplot(1,2,2)
imshow(img)

%% Watershed, yet again
clear
close all

img = rgb2gray(imread('DAPI.bmp'));

figure
imshow(img,[])

%% Phase I:
% The very first thing we want to do is find the markers for the watershed.
%  We have our arsenal of image manipulations to go after this goal.

% First lets try a mean and median filter to see what looks cleaner
% (getting rid of the speckle)
img1 = imfilter(img,fspecial('average',9));
img2 = medfilt2(img,[9 9]);

figure
imshow([img img1 img2])
% I like the median filter a little bit better, so we'll stick with that
% for now.

% Then I'm going to use adapthisteq to try and pull out some of the dimmer
% nuclei.  This stretches the histogram across subimages in the pciture.
img3 = adapthisteq(img2);
figure
imshow([img img2 img3])

% I think the cells look pretty reasonable, so now lets open it to smooth out the bumps.
% The nuclei are all rougly circles, so I'll use the disk structured
% element.  Lets try a couple of radii.
img4 = imopen(img3,strel('disk',3));
img5 = imopen(img3,strel('disk',7));

figure
imshow([img3 img4 img5])

% Not enough, lets go a bit bigger
img6 = imopen(img3,strel('disk',13));
img7 = imopen(img3,strel('disk',17));

figure
imshow([img3 img6 img7])

% Now lets smooth that out a bit
img8 = imfilter(img7,fspecial('gaussian',11));
img9 = imopen(img8,strel('disk',11));

figure
imshow([img img8 img9])

% Lets try a regional max.  We should aim to get one marker per nuclei.  If
% there is more than one marker, then it will become two objects in our
% watershed
rMax = imregionalmax(img8);
figure
imshowpair(rMax,img3)

% That looks pretty good.  Maybe we can erode by a bit so the markers are
% all smaller
rMax2 = imerode(rMax,strel('disk',9));
imshow([rMax rMax2])

%%
% Phase II - Now lets set a gray threshold to identfiy foreground from
% background.  I'm going to use the adapted image since that brought out
% the background cells a bit.
bw = im2bw(img3,graythresh(img2)*.9);

% Use imclose to soften up the sharp edges a bit.
bw2 = imclose(bw,strel('disk',5));

% bwareaopen will get rid of any connected objects smaller than the second
% arguement.  I'll use 400 to get rid of the really small things that are
% clearly not nuclei.
bw3 = bwareaopen(bw2,400);

figure
imshow([bw bw2 bw3])

figure
imshowpair(bw3,rMax2)

bw_perim = bwperim(bw3);
figure
imshowpair(bw3, bw_perim | rMax2)

%% Phase 3
% Now lets take the complement of the original image, and set our markers
% as the minima.
imgComp = imcomplement(img3);
I_low_valleys = imimposemin(imgComp, ~bw3 | rMax);

figure
imshow(I_low_valleys)

% Now lets apply watershed
L = watershed(I_low_valleys);
Lcol = label2rgb(L,'jet','w','shuffle');
% Now we can visualize

figure
subplot(1, 2, 1)
imshow(img3)
subplot(1, 2, 2)
imshow(img3)
hold on
himage = imshow(Lcol);
himage.AlphaData = 0.2;
