%% Exercise - Segmentation in images 

% Estimate the relative amounts of grains of different sizes.
% The image has 284x367 pixels, the true width is 3 cm, height=3x284/367 =
% 2.32 cm. The true number of grains in the image is 236, including three
% grains that overlap the borders of the image. 

% read image 
I1 = imread('../../Users/sergeikatsev/Documents/grainsize.tif');
ix=3; iy=284*3/367;
imshow(I1,'Xdata',[0 ix],'Ydata',[0 iy])
title 'Original image'
% Here, ix and iy denote the coordinate axes used to calibrate the image to
% a centimeter scale. 
axis on

%%
% Reject the color information by converting the image to grayscale.
I2=rgb2gray(I1);
figure
imshow(I2,'Xdata',[0 ix],'Ydata',[0 iy])
title 'Grayscale image'
%%
% Use the function imajust to improve the contrast by adjusting the
% intensity values. This function maps the intensity values to new values
% such that 1% of the data is saturated at low and high intensities.
I3=imadjust(I2);
figure
imshow(I3,'Xdata',[0 ix],'Ydata',[0 iy])
title 'Adjusted intensity image'
%%
% To separate the grains from the background, we need to determine what the
% background is. The function imopen(im, se) determines objects in an image
% 'im' below a certain pixel size and a flat structuring element 'se' such
% as a disk with a radius of 3 pixels generated by the function 'strel'. We
% then produce a background-free image:
I4=imopen(I3,strel('disk',1));
figure
imshow(I4,'Xdata',[0 ix],'Ydata',[0 iy])
title 'No background'
%%
% Subtract the background from the original image to observe the background
% that has been eliminated:
I5=imsubtract(I3,I4);
figure
imshow(I5,'Xdata',[0 ix],'Ydata',[0 iy])
title 'Background'
%%
% The function 'im2bw' converts the background-free image I4 to a binary
% image I6 by thresholding. If the threshold is 1.0 the image is black, if
% the background is 0.0 the image is all white.
I6=im2bw(I4,0.2);
figure
imshow(I6,'Xdata',[0 ix],'Ydata',[0 iy])
title 'Binary image'
%%
% We next eliminate the objects that overlap the image border.
I7=imclearborder(I6);
figure
himage1=imshow(I6,'Xdata',[0 ix],'Ydata',[0 iy]);hold on
himage2=imshow(imsubtract(I6,I7),'Xdata',[0 ix],'Ydata',[0 iy]);
set(himage2,'AlphaData',0.4);
title 'Image border', hold off
%%
% Now we can trace the boundaries using 'bwboundaries' in a binary image
% where non-zero pixels belong to an object whereas zero pixels are
% background. The option 'noholes' supresses the tracing of the holes.
% Function 'label2grb' converts the label matrix L resulting from
% bwboundaries to an RGB image. We use the colormap 'jet'.
[B, L]=bwboundaries(I7,'noholes');
figure
imshow(label2rgb(L,@jet,'w','shuffle'),'Xdata',[0 ix],'Ydata',[0 iy])
title 'Define objects'
%%
% The function 'bwlabeln' obtains the number of connected objects in the binary image.
[labeled, numObjects]=bwlabeln(I7,8);
numObjects
% Compare the produced number to the true number. See if reducing the disk
% size in 'strel' from disk=1 to disk=5 can help separate connected grains. 
%%
% Now determine the areas of the grains
graindata=regionprops(labeled,'basic');
grainareas=[graindata(:) .Area];
objectareas=3^2*grainareas*367^(-2);
%%
% Now we can find the maximum, minimum, and mean areas for all grains in
% the image, in cm^2.
max_area = max(objectareas)
min_area = min(objectareas)
mean_area = mean(objectareas)
%%
% Finally, we plot the histogram of all the grain areas.

clf
v= 0.0005 : 0.0005 : 0.15;
hist (objectareas,v)
xlabel ('Grain Size in Millimeters^2')
ylabel ('Number of grains')
axis([0 0.1 0 30])