%% Morphological operators
% Morphological operators involve moving around a chunk that is smaller
% than the size of the images.  


% Creating structural elements
clear 
close all

% To create a structural element, you use the strel function.  The first
% argument is the shape of the filter, and the second is the size.  
filt1 = strel('disk',4)
filt1.Neighborhood
% Creates a structured element in the shape of a disk with radius 4 pixels

filt2 = strel('square',4)
% Obviously creates a square with side length 4

% If you want to make a line, now you need three total arguments, the
% 'line', the length, and the angle in degrees
filt3 = strel('line',7,90)
filt4 = strel('line',7,45)

% A rectangle also needs another argument.  The first is the number of
% rows, the second is the number of columns
filt5 = strel('rectangle',[2 5])

% Note that for each of these variables, the type of varible is not a
% matrix, but a 1x1 strel.  There is a specific variable that Matlab needs
% in order to work with structured elements.


%% Now that we know how to make structured elements, we can use them to morphologically change images

% First load in the image and change it to a double
img = imread('blobs.png');
imgDouble = im2double(img);

figure
imshow(imgDouble)

% First lets try a dilation.  We'll use the disk from above (filt1)
imgDil1 = imdilate(imgDouble,filt1);

figure
imshow([imgDouble, imgDil1])
title('Dilation with a disk')

% This has the effect of making everything bigger.  We can try with
% different filters

imgDil2 = imdilate(imgDouble,filt2); % Square
imgDil3 = imdilate(imgDouble,filt4); % 45 degree line

figure
imshow([imgDouble, imgDil1; imgDil2 imgDil3]);


%% We can also erode objects in exactly the same way.

imgErode1 = imerode(imgDouble,filt1);
imgErode2 = imerode(imgDouble,filt2); % Square
imgErode3 = imerode(imgDouble,filt4); % 45 degree line

figure
imshow([imgDouble, imgErode1; imgErode2 imgErode3]);


% Note, we can also dilate and erode grayscale images

img2 = imread('coins.png');
img2Double = im2double(img2);

figure
imshow(img2Double);

img2Dil1 = imdilate(img2Double,filt2);
img2Erode1 = imerode(img2Double,filt2);

figure
imshow([img2Double,img2Dil1, img2Erode1])

% You can see that on a grayscale image, dilation makes things overall
% brighter, and erosion makes things overall darker.  The shadows become
% much more pronounced.


%% 
% We can also use dilation and erosion to open and close images.
clear
close all

img = im2double(imread('rice.png'));
filt1 = strel('disk',4);

% Opening consists of an erosion, followed by a dilation.
tmpimg = imerode(img,filt1);
imgOpen = imdilate(tmpimg,filt1);

figure
imshow([img,tmpimg,imgOpen])


% We can also close by dilating, followed by eroding
tmpimg2 = imdilate(img,filt1);
imgClose = imerode(tmpimg2,filt1);

figure
imshow([img,tmpimg2,imgClose])

% Note again that closing has the effect of filling in regions between
% grains, whereas opening separates them.  

% That was the hardway, we can also just use imopen, imclose
imgOpen2 = imopen(img,filt1);
imgClose2 = imclose(img,filt1);

figure
imshow([imgOpen,imgOpen2;imgClose,imgClose2])


%% Load in the grains of rice, remove background
clear
close all

img = imread('rice.png');
imgDouble = im2double(img);

figure
imshow(imgDouble)
title('Background varies along y axis')

% Lets plot out the background along the y axis.
yplot = mean(imgDouble,2);
figure
plot(yplot)% We can see that the background is varying across the image.  To estimate
% the background, we want to remove the grains of rice, and just look at
% everything else.  


% We will use imopen (erosion, followed by dilation) to try to remove 
% grains.  If we make our strel % larger than the size of the grains of 
% rice, it will replace them with levels close to  the background.

filt1 = strel('disk',5);
filt2 = strel('disk',9);
filt3 = strel('disk',13);

% Rice grains are around 10 pixels along the widht.  We want to choose an
% element that is large enough not be contained in a grain.

bg1 = imopen(imgDouble,filt1);
bg2 = imopen(imgDouble,filt2);
bg3 = imopen(imgDouble,filt3);

% Show the regular image and the background
figure
imshow([imgDouble bg1; bg2 bg3])

% Looking at this image, it is clear that 5 is too small, b

%%

% The imopen command removes anything smaller than the size of the disk.
% Since the radius of the grains of rice is ~8 pixels, they get removed.

% Now we can subtract the background
imgCorr = imgDouble - bg2;
figure
imshow([imgDouble imgCorr])

% And adjust the contrast to display it more clearly.
imgCorr2 = imadjust(imgCorr);
figure
imshow([imgCorr imgCorr2])

%% You can also use imopen to get rid of small specks in the background

% Lets make a new filter with a radius smaller than a grain of rice
filtSmall = strel('disk',4);
imgSpeck = imopen(imgCorr2,filtSmall);

figure
subplot(2,1,1)
imshow([imgCorr2 imgSpeck])
subplot(2,1,2)
imshowpair(imgSpeck,imgDouble)
% Now lets erode by a little bit to ensure each grain is isolated
filt5 = strel('square',2);
imgErode = imerode(imgSpeck,filt5);

figure
imshow([imgDouble imgSpeck imgErode])

% Now we can set a new threshold and find connected regions with region
% props.

figure
imhist(imgErode)

thresh = graythresh(imgErode);
imgBW = imgErode > thresh;
figure
imshow([imgSpeck imgBW])

img_label = bwlabel(imgBW);
img_labelRGB = label2rgb(img_label,'spring','c','shuffle');

figure
imshow([img_labelRGB])

% Note, if you want to find the intensitites of each grain, you have to
% reference the original image, not the filtered image.  We used that just
% to segment each grain.

stats = regionprops(imgBW,imgCorr,'Area','MeanIntensity');
areas = [stats.Area];
intens = [stats.MeanIntensity];

figure
subplot(2,2,1)
hist(areas,25)
title('Grain areas')
subplot(2,2,2)
hist(intens,25)
title('Grain intensities')
subplot(2,2,[3 4])
scatter(areas,intens)
xlabel('Areas')
ylabel('Intensities')

%% Looking at the grain areas, there are still some that are too small.  
% Those are the ones touching the border.  If we are worried about
% consistent size measurements, we can remove all the cells on the border.

img_clearBorder = imclearborder(imgBW);

% Note that all the objects near the edges are now removed.
figure
imshow([imgBW img_clearBorder])

% Now we can re-run our stats and see how we look.
stats2 = regionprops(img_clearBorder,imgCorr,'Area','MeanIntensity','Eccentricity','Perimeter');
areas2 = [stats2.Area];
ecc2 = [stats2.Eccentricity];
intens2 = [stats2.MeanIntensity];
peri2 = [stats2.Perimeter];

figure
subplot(2,4,1)
hist(areas2,20)
title('Areas')
subplot(2,4,2)
hist(intens2,20)
title('Intensity')
subplot(2,4,3)
hist(ecc2,20)
title('Eccentricity')
subplot(2,4,4)
hist(peri2,20)
title('Perimeter')
subplot(2,4,5:6)
scatter(areas2,intens2)
xlabel('Areas')
ylabel('Intensity')
subplot(2,4,7:8)
scatter(areas2,peri2)
xlabel('Areas')
ylabel('Perimeter')


