%% Grains of rice (See lesson 8)
clear
close all

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


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

stats = regionprops(imgMask3,imgDouble,'Area','MeanIntensity');
stats(11:20).Area

%% 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];
histogram(areas,length(areas))
xlabel('Size of objects (pixels)')
ylabel('Number of objects')

figure 
imshow(img)

% 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
histogram([stats2.Area])

justRight2 = find(areas > 10 & areas < 300);
stats3 = stats(justRight2);
figure
histogram([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);

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 == [currentGrain]);
    maskJustRight = maskJustRight + masktmp;
end

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





