% Ridler Calvard

clear all
close all

% Go to right folder (replace with your catalogue)

cd 'C:\Documents and Settings\aurdal\My Courses\UIO\INF3300\'

% Noisy square

i=150*ones(256);
%i(:,129:256)=200;
i(:,220:256)=200;
figure
imshow(i,[0 255])

% Add noise

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

% Start by calculating the initial mean and max value

t_old=round(mean(mean(i)));
G=max(max(i));

% Generate histogram and display

ih=histc(i(1:prod(size(i))),0:G);
figure
plot(ih)

% Initialise

stop=0;
T=[t_old];
index=2;

% Now iterate

while(stop==0)
    mu_1=sum((1:t_old).*(ih(1:t_old)))/(sum(ih(1:t_old)));
    mu_2=sum(((t_old+1):G).*(ih((t_old+1):G)))/(sum(ih((t_old+1):G)));
    t_new=round((mu_1+mu_2)/2);
    if(t_new==t_old)
        stop=1;
    end
    t_old=t_new;
    T(index)=t_old;
    index=index+1;
end

% Display image thresholded at the selected value

figure
imshow(i>t_new,[]);

figure(3)
a=axis;
l=line([t_new t_new],[a(3) a(4)]);
set(l,'color','r')



% Otsu

clear all
close all

% Go to right folder

cd 'C:\Documents and Settings\aurdal\My Courses\UIO\INF3300\'

% Histogram of letter

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

% Generate histogram and display

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

% Initialise L

L=256;

% Generate a normalized histogram and display

ihn=ih/sum(ih);
figure
plot(ihn)
sum(ihn)

% Now look at what happens at different values for k

k=130; % this is where you can set different values

% Show thresholded image

figure
imshow(i<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+ihn(n);
end
omega_0
omega_1=1-omega_0

% Calculate and display mus

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

% Calculate sigmas

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

% 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.....

