im_color = imread('fox.jpg');
im_gray = rgb2gray(im_color);

[X, Y] = size(im_gray);
[intensity,x] = imhist(im_gray);
% Normalizálhatnám a hisztogramot, de nem változtat az eredményen:
% intensity = intensity / (X*Y); 

% A bal oldalon szereplő változókba bemásolja a zeros(256,1) tömböt
[w1, w2, u1, u2, sigma] = deal(zeros(256,1));

% Kihasználom, hogy az értékek kumulálódnak, vagyis hogy nem kell
% kiszámolnom az egész szummát a képletekben minden egyes k-ra, hanem az
% előző értékből számolható
for k = 2:256
    w1(k) = w1(k-1) + intensity(k);
    u1(k) = (u1(k-1)*w1(k-1) + intensity(k)*k) / w1(k);
    w2(256-k+1) = w2(256-k+2) + intensity(256-k+1);
    u2(256-k+1) = (u2(256-k+2)*w2(256-k+2) + intensity(256-k+1)*(256-k+1)) / w2(256-k+1);
end
for k= 2:256
    sigma(k) = w1(k) * w2(k) * (u1(k)-u2(k))^2;
end
bar(u1)

[val, place] = max(sigma)

thresholded = im_gray;
thresholded(im_gray < place) = 0;
thresholded(im_gray >= place) = 255;

% figure(1)
% subplot(2,2,1)
% imshow(im_color)
% title('Original')
% 
% subplot(2,2,2)
% imshow(thresholded)
% title('my Otsu')
% 
% subplot(2,2,3)
% thresh = multithresh(im_gray,4);
% seg_I = imquantize(im_gray,thresh);
% imshow(seg_I, [])
% title('built in')