%% Segmenting - Look at image with rice grains
clear
close all

img = imread('rice.png');

figure
imshow(img)

% Convert it to a double
imgDouble = im2double(img);

% Looking at the image, there are grains of rice that are clearly visible
% by eye.  We might want to know how many grains of rice there are.  The
% easiest thing to do is set a simple threshold.
figure
imhist(imgDouble)

% There are three main peaks.  The lowest peak (around 0.2) is clearly
% background, but we might ask if the grains of rice are located around
% 0.4, of if they are around 0.7.  
%
% To test this hypothesis, we can make a mask to see if they make sense.

thresh1 = 0.35;
thresh2 = 0.65;

imgMask1 = imgDouble > thresh1;

figure
imshow([imgDouble,imgMask1])
title(['Threshold with ' num2str(thresh1)])

% Now try the second threshold
imgMask2 = imgDouble > thresh2;
figure
imshow([imgDouble,imgMask2])
title(['Threshold with ' num2str(thresh2)])

% This is close, but it's missing the grains on the bottom since the
% background is different.
%
% We can also let Matlab choose a threshold for us

thresh3 = graythresh(imgDouble);
imgMask3 = imgDouble > thresh3;

figure
imshow([imgDouble imgMask1; imgMask2 imgMask3])

% The matlab selected threshold looks prety good, but you do pick up some
% background spots and miss some on the bottom.

%% Identifying connected regions
% Regionprops is an awesome command that will identify all the connected
% regions in a logical image and bin them together. This is a great way to
% count things in an image.

stats = regionprops(imgMask3,'Area'); % Returns the area in pixels

% stats is a structure with elements Area
stats(1).Area 
stats(2).Area
stats(1:10).Area

% Can also grab other features
stats = regionprops(imgMask3,imgDouble,'Area','MeanIntensity');

stats(5).Area
stats(5).MeanIntensity

%% Removing elements that are definitely wrong
% There are 151 objects that regionprops found.  But this is higher than
% the number of actual grains of rice.  Some of the areas are just 1 pixel.
%  These are almost certainly background pixels.
%
% We can look at the distribution of sizes
% First we put the structure into an array.

areas = [stats.Area];
hist(areas,length(areas))
xlabel('Size of objects (pixels)')
ylabel('Number of objects')

% You can see that a large number of objects have size < 10.  We know those
% are not grains of rice.  Lets get rid of those.

justRight = find(areas > 10);

% This returns a list of indices.  For example.  Objects 9, 10, 11, 18 etc
% were all too small.  We want to get rid of those

stats2 = stats(justRight);
figure
hist([stats2.Area])

justRight2 = find(areas > 10 & areas < 300);
stats3 = stats(justRight2);
figure
hist([stats3.Area])

% We can save all the variables that we just generated so we don't have to
% start from scratch.
save rice_stats.mat

% If we clear everything...
clear
close all

% ... we can use the load command to bring back all of those variables
load rice_stats.mat

% We can also select only certain variables to save
save rice_stats_trunc stats stats2 

% Then we can display everything just like we did before.  This is very
% convenient if the data took a long time to process.
figure
imshow(imgMask3);

%% Start here for new lesson
clear
close all

% ... we can use the load command to bring back all of those variables
load rice_stats.mat

% Then we can display everything just like we did before.  This is very
% convenient if the data took a long time to process.
figure
imshow(imgMask3);


%% Labeling black and white images
% Starting from our mask before, we can have Matlab label each of the
% connected elements
img_label = bwlabel(imgMask3);

figure
imshow(img_label)

% Each element gets its own number, and the background is 0. So the first
% grain has values of 1, the second 2, etc...  Convientent to view the
% results

img_labelRGB = label2rgb(img_label);
figure
imshow(img_labelRGB)

% If we want to look at just the 5th item.
maskSingleGrain = (img_label == 85);
figure
imshow(maskSingleGrain)

% Use a for loop to identify all the grains that are just right
numJustRight = length(justRight);

% Declare a new mask which is all zeros, exactly the size of our old mask
maskJustRight = zeros(size(imgMask3));
for j = 1:numJustRight
    currentGrain = justRight(j);
    masktmp = (img_label == [1 2]);
    maskJustRight = maskJustRight + masktmp;
end

% Overlay two images with the command imshowpair (first image green, second
% is pink)
figure
imshowpair(imgDouble, maskJustRight)

%% Looking at our results with hist and scatter
% We'll start with stats2, where we removed those things that were
% definitely too small.

areas2 = [stats.Area];
intens2 = [stats.MeanIntensity];

figure
subplot(2,1,1)
hist(areas2)
title('Rice areas')
subplot(2,1,2)
hist(intens2)
title('Rice intensities')

% A second argument in the hist command will let us set the number of bins
figure
subplot(2,1,1)
hist(areas2,20)
title('Histogram with 20 bins')
subplot(2,1,2)
hist(areas2,length(areas2))
title('Histogram with 100% bins')

% Scatter will let us compare two quantities togther
figure
scatter(areas2,intens2)
xlabel('Rice Areas')
ylabel('Rice Intensitiy')

% There are two populations in the intensity.  Lets see if we can look at
% those below 0.62.

dim_grains = find(intens2 < 0.62)

% Comes back with a list of 20 grains.  We'll use the trick of a for loop
% adding into a new mask to generate the dim grains.
nDim = length(dim_grains);
dimMask = zeros(size(img));
for j = 1:nDim
    tmpMask = (img_label == dim_grains(j));
    dimMask = dimMask + tmpMask;
end

% Display all the dim grains
figure
imshow([dimMask imgDouble])