%% Setup - run this first

clear all; close all; clc; %#ok<CLALL> 

% MATLAB has a few built-in images that you can load to test with
rgb_pears = imread('pears.png');
rgb_chips = imread('coloredChips.png');

% Different formats can be read by the same function, but the data type
% that results can vary. 
I_rice = imread('rice.png');
I_coins = imread('coins.png'); % Note that imread takes a variety of different extensions, but the data type isn't necessarily the same. 



%
%% Segmentation

% Segmentation is the process of separating a whole into its parts. The
% simplest example of this is background/foreground, or black and white.

imAndBW = @(x) montage({x,imbinarize(x)});

imAndBW(rgb_pears);
imAndBW(rgb_chips); % Notice the yellow and orange chips are missing. Also notice that background is bright.
imAndBW(I_rice); % Notice the missed grains towards the bottom. 
imAndBW(I_coins);

% Generally you will want a threshold that is different from image to
% image. Otsu's method assumes a single threshold but infers it from image
% statistics, and is the method behind the imbinarize function. 

imAndOtsu = @(x) montage({x,imbinarize(x,graythresh(x))})

imAndOtsu(rgb_pears); % Notice that a single threshold was applied to all channels, not per-channel.
imAndOtsu(rgb_chips); % Same
imAndOtsu(I_rice);
imAndOtsu(I_coins);

% Compare im2bw and otsu method
imBWAndOtsu = @(x) montage({x,imbinarize(x),imbinarize(x,graythresh(x))}, 'Size', [1 3])
imBWAndOtsu(rgb_pears); % Notice that a single threshold was applied to all channels, not per-channel.
imBWAndOtsu(rgb_chips); % Same
imBWAndOtsu(I_rice);
imBWAndOtsu(I_coins);


%
%% Adaptive thresholding

% It is possible to derive a locally adaptive threshold and apply that to
% the image, to compensate for spatially-varying illumination or other
% variation. This has interesting effects if applied to multichannel
% images.

imAndAdapt = @(x) montage({x,imbinarize(x,adaptthresh(x))});

imAndAdapt(rgb_pears); % 
imAndAdapt(rgb_chips); % Same
imAndAdapt(I_rice); % Notice that more grains were captured, but more noise was as well.
imAndAdapt(I_coins); % The foreground of the coins is less well-defined.


%

% If you have an external boundary of objects in your image, it is
% straightforward to fill in the holes that remain.

coin_adaptive_threshed = imbinarize(I_coins,adaptthresh(I_coins));
filled_coin_adaptive_threshed = imfill(coin_adaptive_threshed,'holes');
montage({coin_adaptive_threshed,filled_coin_adaptive_threshed});

% Show an overlay. Magenta is the second image, green is the first, and
% white is the overlay of both.
imshowpair(coin_adaptive_threshed,filled_coin_adaptive_threshed);


% It is important to clean up your image as well as you can before trying
% to locate features and content!
bounds_coin = bwboundaries(coin_adaptive_threshed);
bounds_filled = bwboundaries(filled_coin_adaptive_threshed);

% Unfilled image
figure; colormap('gray');
imagesc(coin_adaptive_threshed); axis equal tight off; hold on;
% Plot the boundaries.
for k = 1:length(bounds_coin)
   boundary = bounds_coin{k};
   plot(boundary(:,2), boundary(:,1), 'r', 'LineWidth', 2)
end
title('Boundaries of unfilled')

% Filled image
figure; colormap('gray')
imagesc(filled_coin_adaptive_threshed);  axis equal tight off; hold on;
% Plot the boundaries.
for k = 1:length(bounds_filled)
   boundary = bounds_filled{k};
   plot(boundary(:,2), boundary(:,1), 'r', 'LineWidth', 2)
end
title('Boundaries of filled')

%
%% Getting information about each region

% There are many built-in functions that can assist in quantifying your
% image regions. This is complementary to the ROI methods shown in the next
% notebook.
regprops_basic=regionprops(filled_coin_adaptive_threshed);
regprops_all=regionprops(filled_coin_adaptive_threshed,'all');

% Each element of the regprops corresponds to a region within the image.
% bwconncomp yields pixel labels in an object-oriented way.
connectedComp = bwconncomp(filled_coin_adaptive_threshed);

% bwlabel gives back an image with each region labeled with a different
% gray value.
connLabels = bwlabel(filled_coin_adaptive_threshed);
figure;
imagesc(connLabels); title('Indexed regions'); axis equal tight off; colorbar;


rgb_labels = label2rgb(connLabels);
figure;
imagesc(rgb_labels); title('Labeled regions'); axis equal tight off;
for k = 1 : numel(regprops_basic)
    text(regprops_basic(k).Centroid(1),regprops_basic(k).Centroid(2),num2str(k))   
end

% Notice that the connected components include a very small region (region
% 8) which would be an outlier in our dataset.


% MATLAB variables can be saved in .mat files for later examination.
save('./coin_regionprops.mat','regprops_all')
regprops_all(3).Centroid
regprops_all(3).EquivDiameter