rosinthreshold.m

image thresholding otsu method

function T = RosinThreshold(imhist, picknonempty)

% should I ensure the chosen(threshold) bin is non empty?
if (nargin < 2)
    picknonempty = 0;
end

% find best threshold

[mmax2, mpos] = max(imhist);
p1 = [mpos, mmax2];

% find last non-empty bin
L = length(imhist);
lastbin = mpos;
for i = mpos:L
    if (imhist(i) > 0)
        lastbin=i;
    end
end
    
p2 = [lastbin, imhist(lastbin)];
DD = sqrt((p2(1)-p1(1))^2 + (p2(2)-p1(2))^2);

if (DD ~= 0)
	best = -1;
	found = -1;
	for i = mpos:lastbin
        p0 = [i,  imhist(i)];
        d = abs((p2(1)-p1(1))*(p1(2)-p0(2)) - (p1(1)-p0(1))*(p2(2)-p1(2)));
        d = d / DD;
        
        if ((d > best) && ((imhist(i)>0) || (picknonempty==0)))
            best=d;
            found = i;
        end
	end
	
	if (found == -1)
        found = lastbin+1;
	end
else
    found = lastbin+1;
end

% plot(imhist);
% hold on;
% line([p1(1) p2(1)],[p1(2) p2(2)]);
% p3 = [found imhist(found)];
% m1 = (p2(2) - p1(2)) / (p2(1) - p1(1))
% m2 = -0.7*m1
% c0 = p1(2) - m1*p1(1)
% c = p3(2) - m2*p3(1)
% x4 = (c0-c) / (m2-m1);
% y4 = m2*x4+c;
% p4 = [x4 y4];
% line([p3(1) p4(1)],[p3(2) p4(2)]);
% plot(lastbin,0,'ro')
% plot(mpos,0,'go')
% plot(found,0,'ks')
% hold off
% pause

%T = found;
T = min(found, L);