%{
Code by:
Arjit Soni
B.Tech Student
Department of Civil Engineering 
and 
Department of Engineering Science 
Indian Institute of Technology Hyderabad 
es13b1005@iith.ac.in

Under Guidance of:
Dr. B. Umashankar
Associate Professor
Department of Civil Engineering
Indian Institute of Technology Hyderabad
buma@iith.ac.in
%}
%% Read me
%This code should be in same folder where images are.
%select only 10-12 pixels in both kerosene and water region.

%% Clear workspace
clear all               %to clear workspace
close all               %to close all windows
clc                     %to clear command window

%% Image input dynamically image should be in same folder as code 
in = input('Type input image name:','s');
a=sscanf(in,'%s');      %scan name of the image
I = imread(a);          %to read image
img = im2double(I);     %to convert uint8 to double
whos img                %to show image type

%% to make rgb layer
red = img(:,:,1);       %to extract red layer from image
green = img(:,:,2);     %to extract green layer from image
blue = img(:,:,3);      %to extract blue layer from image

%% to select pixel colour
figure('Name','To select pixels (double click to finish)','NumberTitle','off');   %to name a figure
imshow(img);
pix = impixel(img);     %to select pixel colour
disp('     Red      Green      Blue');
disp(pix);
mi = min(pix);          %minimum value of pixels
ma = max(pix);          %maximum value of pixels
disp('minimum values');
disp(mi);
disp('maximum values');
disp(ma);

%% taking minimum and maximum of pixel values
outr = red>mi(1) & red<ma(1);       %extracting red layer using min and max pixel values
outg = green>mi(2) & green<ma(2);   %extracting green layer using min and max pixel values
outb = blue>mi(3) & blue<ma(3);     %extracting blue layer using min and max pixel values
out = outr & outg & outb;           %using AND funcion on above extracted layers

%% displaying final images
figure('Name','Normal and RGB  images','NumberTitle','off') %to name a figure
subplot(2,2,1), imshow(img);        %to show normal image
title('original');
subplot(2,2,2), imshow(outr);       %to show red layer of image
title('red');
subplot(2,2,3), imshow(outg);       %to show green layer of image
title('green');
subplot(2,2,4), imshow(outb);       %to show blue layer of image
title('blue');
figure('Name','processed image','NumberTitle','off')
imshow(out);

%% Morphological Operations on Binary Image
holes = imfill(out,'holes');
erode = bwmorph(holes,'erode');
dilate = bwmorph(erode,'dilate');
figure('Name','After Morphological Operations','NumberTitle','off') %to name a figure
subplot(2,2,1), imshow(out);        %to show normal image
title('original');
subplot(2,2,2), imshow(holes);
title('after filling holes');
subplot(2,2,3), imshow(erode);
title('after erode');
subplot(2,2,4), imshow(dilate);
title('after dilate');

outf = dilate;         %initiating the final binary image

%% finding number of pixel in given area
figure('Name','select water region','NumberTitle','off')%to name a figure
imshow(outf);
freew = imfreehand;                 %for drawing the region using imfreehand tool
maskw = createMask(freew);          %creating a mask using imfreehand function
outw = maskw.*outf;                 %extracting values inside mask on final image
numw = sum(outw(:));                %using sum to find number of white pixels inside water

figure('Name','select kerosene region','NumberTitle','off')%to name a figure
imshow(outf);
freek = imfreehand;                 % for drawing the region using imfreehand tool
maskk = createMask(freek);          %creating a mask using imfreehand function
outk = maskk.*outf;                 %extracting values inside mask on final image
numk = sum(outk(:));                %using sum to find number of white pixels inside kerosene

%% Finding Free Swell Index
disp('number of white pixels inside water region');
disp(numw);
disp('number of white pixels inside kerosene region');
disp(numk);
disp('');
ratio = numw/numk;
ratio = round((ratio*100))/100;     %rounding off the percentage increase to two digits
FSI = (ratio-1)*100;
if FSI < 0              %using 'if' to make sure FSI value is not negative
    FSI = 0;
end
disp('the Free Swell Index is')
disp(FSI);

%% END
%close all               %to close all windows
%clc                     %to clear command window
%disp('the Free Swell Index is ')
%disp(FSI);