%% Writing custom functions
clear
close all

img = im2double(imread('rice.png'));
bigMask = findBigGrains(img);

% Suppose we want to remove the background with imopen and find the largest
% 10% of grains

% Start by removing the slow background drift and noisy pixels
bg = imopen(img,strel('disk',11));
imgNew = img - bg;
imgNew2 = imopen(imgNew,strel('disk',3));
imgNew3 = mat2gray(imgNew2);

% Set a threshold and identify individual grains
imgMask = im2bw(imgNew3,graythresh(imgNew3));
imgMask2 = imopen(imgMask,strel('disk',4));
imgMask3 = imerode(imgMask2,strel('square',2));

figure
imshow([img (imgMask2)])

imgSegm = bwlabel(imgMask2,8);
imgLabel = label2rgb(imgSegm);
imshow(imgLabel)

stats = regionprops(imgSegm,'area');

% Collect cell areas and take top 10%
areas = [stats.Area];
bigGrains = find(areas > prctile(areas,90));
bigGrainStats = stats(bigGrains);

indx = 0;
bigMask = zeros(size(img));
for j = 1:length(bigGrains)
    tmpMask = (imgSegm == bigGrains(j));
    bigMask = bigMask + tmpMask;
end

figure
imshowpair(img,bigMask)
% If this is a set of commands that we will run on many images in a row, it
% becomes much cleaner and easier to make this whole routine into a
% function.

%%
clear
close all

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

% Because we included the commented text, we can type help and see what is
% known about this function
help findBigGrains

% Now we can call findBigGrains function
[bigGrains,imgLabel] = findBigGrains(img,10);

figure
imshowpair(img,imgLabel)

% Aside note, make sure that functions you want to use regularly are in a
% path searched by Matlab.

