% Clear everything

clear all
close all

% Go to your image catalog and read your image

cd 'C:\Documents and Settings\aurdal\My Courses\UIO\INF3300' %Exchange with your own path
i=imread('monalisa.jpg');

% Take a look at the worlds most famous woman

imshow(i)

% Threshold at all values greater than a limit

it=i>10;

% Open a new figure window and display

figure
imshow(it)

% Band thresholding is just as simple

itb=(i<50)|(i>100);

% Open a new figure window and display

figure
imshow(itb)

% Semithresholding

its=double(i).*double((i>100));

% Open a new figure window and display

figure
imshow(its,[min(min(its)) max(max(its))])

% Smart colormap manipulation

map=gray(256);
map(1,:)=[1 0 0]
colormap(map)

% Tests on disks, load image

clear all
close all

i=imread('disk.png'); % Try with disk1 also, it'a a smaller disk
i=double(i(:,:,1));
figure;
imshow(i,[0 255])

% Add low noise

in=i+4*randn(size(i));
figure;
imshow(in,[0 255])

% Generate histogram

inh=histc(in(1:prod(size(in))),0:255);

figure
bar(inh)

% Otsu's method

% Read image, add some noise and show histogram

clear all
close all

% i=imread('disk.png');
% i=double(i(:,:,1));
% figure;
% imshow(i,[0 255])
% in=i+5*randn(size(i));
% figure;
% imshow(in,[0 255])

i=imread('francis.jpg');
i=rgb2gray(i);
figure;
imshow(i,[0 255])
in=i;

% Initialise L

L=256;

% Generate a normalized histogram and display

inh=histc(in(1:prod(size(in))),0:255);
inhn=inh/sum(inh);
figure
bar(inhn)
sum(inhn)

% Now look at what happens at different values for k

k=130; % Set this to whatever you like in the range 1....L

% Show thresholded image

figure
imshow(in>=k)

% Initialise

omega_0=0;
omega_1=0;
mu_0=0;
mu_1=0;
mu_T=0;
sigma_0=0;
sigma_1=0;
sigma_T=0;

% Calculate and display omegas

for n=1:k
  omega_0=omega_0+inhn(n);
end
omega_0
omega_1=1-omega_0

% Calculate and display mus

for n=1:k
  mu_0=mu_0+n*inhn(n)/omega_0;
end
for n=k+1:L
  mu_1=mu_1+n*inhn(n)/omega_1;
end
for n=1:L
  mu_T=mu_T+n*inhn(n);
end
mu_0
mu_1
mu_T

% Calculate sigmas

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

% Finally calculate eta

eta=omega_0*omega_1*(mu_1-mu_0)^2/sigma_T

% The interesting thing is obviously to do this for all k

OMEGA_0=[];
OMEGA_1=[];
MU_0=[];
MU_1=[];
SIGMA_0=[];
SIGMA_1=[];
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+inhn(n);
    end
    OMEGA_0(k)=omega_0;
    omega_1=1-omega_0;
    OMEGA_1(k)=omega_1;
    if (omega_0>0)
        for n=1:k
            mu_0=mu_0+n*inhn(n)/omega_0;
        end
    end
    MU_0(k)=mu_0;
    if (omega_1>0)
        for n=k+1:L
            mu_1=mu_1+n*inhn(n)/omega_1;
        end
    end
    MU_1(k)=mu_1;
    for n=1:L
        mu_T=mu_T+n*inhn(n);
    end
    if (omega_0>0)
        for n=1:k
            sigma_0=sigma_0+(n-mu_0)^2*inhn(n)/omega_0;
        end
    end
    SIGMA_0(k)=sigma_0;
    if (omega_1>0)
        for n=k+1:L
            sigma_1=sigma_1+(n-mu_1)^2*inhn(n)/omega_1;
        end
    end
    SIGMA_1(k)=sigma_1;
    for n=1:L
        sigma_T=sigma_T+(n-mu_T)^2*inhn(n);
    end
    ETA(k)=omega_0*omega_1*(mu_1-mu_0)^2/sigma_T;
end

% Plot eta

figure
plot(ETA)
hold on
plot(ETA,'r.')
plot(OMEGA_0,'y')
plot(OMEGA_1,'c')
plot(MU_0,'k')
plot(MU_1,'b')
plot(SIGMA_0,'g')
plot(SIGMA_1,'m')
grid on 


% Niblack's method

% Read image, add some noise and show histogram

clear all
close all

% i=imread('disk.png');
% i=double(i(:,:,1));
% figure;
% imshow(i,[0 255])
% in=i+5*randn(size(i));
% figure;
% imshow(in,[0 255])

i=imread('francis.jpg');
i=double(rgb2gray(i));
figure;
imshow(i,[0 255])

figure
ih=histc(i(1:prod(size(i))),0:255);
bar(ih)

% And the rest is up to you.....