clear all
close all

% Try basic erosion.  Start by defining a structuring element

se=[0 1 0;1 1 1;0 1 0]

% Let's load one of the standard matlab images

load penny

% Threshold this

Pt=P>170;

% Display

imshow(P,[])
figure
imshow(Pt,[])

% Perform a simple erosion

Pte=imerode(Pt,se);

% Display

figure
imshow(Pte,[])

% Try a few other structuring elements

se1=[1 1 1;1 1 1;1 1 1];
se2=[1 1 1 1 1 1 1];

% Try dilation

Ptd=imdilate(Pt,se);

% Display

figure
imshow(Ptd,[]);

% Let's look at closing and opening

Pte=imerode(Pt,se);
Pto=imdilate(Pte,se);

% Display

figure
subplot(1,3,1)
imshow(Pt,[])
subplot(1,3,2)
imshow(Pte,[])
subplot(1,3,3)
imshow(Pto,[])
 
Ptd=imdilate(Pt,se);
Ptc=imerode(Ptd,se);

% Display

figure
subplot(1,3,1)
imshow(Pt,[])
subplot(1,3,2)
imshow(Ptd,[])
subplot(1,3,3)
imshow(Ptc,[])

% Geodesic dilation

clear all
close all

% Start by defining a structuring element

se=[0 1 0;1 1 1;0 1 0]

% Define a synthetic input image

f=zeros(64);
f(32,32)=1;

% Define a mask

g=zeros(64);
g(31:33,32)=1;

% Perform the geodesic dilation of f under mask g

fdil_g=min(imdilate(f,se),g);

% Display

figure
subplot(1,4,1)
imshow(f,[])
subplot(1,4,2)
imshow(g,[])
subplot(1,4,3)
imshow(imdilate(f,se),[])
subplot(1,4,4)
imshow(fdil_g,[])

% Morphological reconstruction

% Here is an image of blood cells:

clear all
close all

i=imread('blood.tif');
i=rgb2gray(i);
figure
imshow(i,[])

% Notice the cell roughly in the middle, it's a white blood cell with
% it's cell nucleus stained.  We want to get a hold of ONLY this one cell
% Hmmmmmmm, let's try a thresholding:

it1=i<100;
figure
imshow(it1,[])

% No luck, this gives us just the nucleus.  Let's change threshold

it2=i<180;
figure
imshow(it2,[])

% This gives us roughly what we want, but notice that now we got many other
% cells too.

% Reconstruction to our aid.  Start by thresholding to make a really strict
% marker for the nucleus (we also use an opening to remove some noise)

it1=i<100;
it1=imopen(it1,ones(5));
figure
imshow(it1,[])

% Now lets make a very "loose" threshold

it2=i<180;
figure
imshow(it2,[])

% Now reconstruct the marker under the mask using dilation.

ir=imreconstruct(it1,it2);
figure
imshow(ir,[])

% Lets finish this off by a closing to remove a few interior holes

ir=imclose(ir,ones(3));
figure
imshow(ir,[])

% This is just a nice trick to put a red border around our cell

irb=bwperim(ir,8);
i(find(irb))=255;
figure
imshow(i,[])
map=gray(256);
map(256,:)=[1 0 0];
colormap(map)

% Simple gray level erosion example

clear all
close all

% Load image of Marilyn Monroe

im=imread('marilyn.tif');
im=im(:,:,1);
figure
imshow(im,[])

ime=imerode(im,ones(3));
figure
imshow(ime,[])

% Often it is very instructive to plot single lines

L=135
figure
plot(im(L,:),'r')
hold on
plot(ime(L,:),'b')
figure(1)
l1=line([0 256],[L L])
set(l1,'color','r')
figure(2)
l2=line([0 256],[L L])
set(l2,'color','b')

% Line orientation example

clear all
close all

% Load image of carbon fibres

im=imread('carbonfibreepoxy.png');
figure
imshow(im,[])

ii=255-double(im);
figure
imshow(ii,[0 255])

index=1;
for m=0:5:180
    dummy=imerode(ii,strel('line',11,m));
    distr(index)=sum(sum(dummy));
    index=index+1;
end

figure
plot(distr/max(distr));
labs = 0:5:180; 
set(gca,'XTick',0:(length(labs)-1),'XTickLabel', sprintf('%d|', labs)); 
grid on
figure
polar((0:5:180)*pi/180,distr)

% Simple gray level dilation example

clear all
close all

% Load image of Marilyn Monroe

im=imread('marilyn.tif');
im=im(:,:,1);
figure
imshow(im,[])

imd=imdilate(im,ones(3));
figure
imshow(imd,[])

% Often it is very instructive to plot single lines

L=135
figure
plot(im(L,:),'r')
hold on
plot(imd(L,:),'b')
figure(1)
l1=line([0 256],[L L])
set(l1,'color','r')
figure(2)
l2=line([0 256],[L L])
set(l2,'color','b')


% Simple gray level closing example

% Not everyone liked Marilyns birthmark.  Can i tbe reomoved?

clear all
close all

% Load image of Marilyn Monroe

im=imread('marilyn.tif');
im=im(:,:,1);
figure
imshow(im,[])

% Let's perform a morphological closing of this image
% Start by dilating, then erode

imd=imdilate(im,ones(3));
figure
imshow(imd,[])

imc=imerode(imd,ones(3));
figure
imshow(imc,[])

% Often it is very instructive to plot single lines

L=174
figure
plot(im(L,:),'r')
hold on
plot(imc(L,:),'b')
figure(1)
l1=line([0 256],[L L])
set(l1,'color','r')
figure(3)
l2=line([0 256],[L L])

% Morphological gradients

clear all
close all

% Load image

i=imread('rice.tif');
figure
imshow(i,[])

% What does sobel produce

s=fspecial('sobel');
is1=imfilter(i,s);
is2=imfilter(i,s');
is=double(is1)+double(is2);
figure
imshow(is1,[min(min(is1)) max(max(is1))])
figure
imshow(is2,[min(min(is2)) max(max(is2))])
figure
imshow(is,[min(min(is)) max(max(is))])

% Let's look at some morphological gradients
% The classical one is the Beucher gradient,
% it is simply the dilation of the input image
% minus the erosion of the input image

id=imdilate(i,[0 1 0;1 1 1;0 1 0]);
ie=imerode(i,[0 1 0;1 1 1;0 1 0]);

ibeuch=imsubtract(id,ie);
figure
imshow(ibeuch,[])

% Another possibility consists in looking at the
% difference between the dilated image
% and the image itself

idil=imsubtract(id,i);
figure
imshow(idil,[])

% A final possibility consists in looking at the
% difference between the image itself
% and the eroded image

iero=imsubtract(i,ie);
figure
imshow(iero,[])

% Lets look at some line plots here

figure

plot(id(100,:),'r')
hold on
plot(ie(100,:),'b')
plot(ibeuch(100,:),'g')

% The white top hat transform

clear all
close all

% Load image and invert it

i=double(imread('francis_gray.jpg'));
i=i(:,:,1);
i=max(max(i))-i;
figure
imshow(i,[])

% Start by opening this image using a fairly large structuring element
% We will use matlabs strel function in order to define a circular
% structuring element

se=strel('disk',11,0)
ie=imerode(i,se);
io=imdilate(ie,se);

% Let's take a look at the opened image

figure
imshow(io,[])

% Line plots are also very interesting here

figure

plot(i(100,:),'r')
hold on
plot(ie(100,:),'b')
plot(io(100,:),'g')

% Now let's look at the wth

iwth=i-io;

figure
imshow(iwth,[]);

% Try a threshold now

figure
imshow(iwth>30)

% Will a GLOBAL application of Otsu work now

% Initialise L

L=256;

% Generate a normalized histogram and display

ih=histc(iwth(1:prod(size(iwth))),0:255);
ihn=ih/sum(ih);
figure
plot(ihn)
sum(ihn)

ETA=[];

for k=1:L
    k
    omega_0=0;
    omega_1=0;
    mu_0=0;
    mu_1=0;
    mu_T=0;
    sigma_0=0;
    sigma_1=0;
    sigma_T=0;
    for n=1:k
        omega_0=omega_0+ihn(n);
    end
    omega_1=1-omega_0;
    if (omega_0>0)
        for n=1:k
            mu_0=mu_0+n*ihn(n)/omega_0;
        end
    end
    if (omega_1>0)
        for n=k+1:L
            mu_1=mu_1+n*ihn(n)/omega_1;
        end
    end
    for n=1:L
        mu_T=mu_T+n*ihn(n);
    end
    if (omega_0>0)
        for n=1:k
            sigma_0=sigma_0+(n-mu_0)^2*ihn(n)/omega_0;
        end
    end
    if (omega_1>0)
        for n=k+1:L
            sigma_1=sigma_1+(n-mu_1)^2*ihn(n)/omega_1;
        end
    end
    for n=1:L
        sigma_T=sigma_T+(n-mu_T)^2*ihn(n);
    end
    ETA(k)=omega_0*omega_1*(mu_1-mu_0)^2/sigma_T;
end

figure
plot(ETA)
hold on
plot(ETA,'*r')
grid on

% Try the max value of ETA as threshold

figure
imshow(iwth>52);

% This weeks exercise

% Find ONLY THE TEXT in the folowing image

clear all
close all

i=imread('unevenillumination.tif');
figure
imshow(i,[])

% And the rest is up to you





