轉自：http://download.csdn.net/detail/hupeng810/1511870

function imagBW = kittlerMet(imag)

% KITTLERMET binarizes a gray scale image 'imag' into a binary image

% Input:

%?? imag: the gray scale image, with black foreground(0), and white

%?? background(255).

% Output:

%?? imagBW: the binary image of the gray scale image 'imag', with kittler's

%?? minimum error thresholding algorithm.

?

% Reference:

%?? J. Kittler and J. Illingworth. Minimum Error Thresholding. Pattern

%?? Recognition. 1986. 19(1):41-47

?

MAXD = 100000;

imag = imag(:,:,1);

[counts, x] = imhist(imag);? % counts are the histogram. x is the intensity level.

GradeI = length(x);?? % the resolusion of the intensity. i.e. 256 for uint8.

J_t = zeros(GradeI, 1);? % criterion function

prob = counts ./ sum(counts);? % Probability distribution

meanT = x' * prob;? % Total mean level of the picture

% Initialization

w0 = prob(1);?? % Probability of the first class

miuK = 0;?? % First-order cumulative moments of the histogram up to the kth level.

J_t(1) = MAXD;

n = GradeI-1;

for i = 1 : n

??? w0 = w0 + prob(i+1);

??? miuK = miuK + i * prob(i+1);? % first-order cumulative moment

??? if (w0 < eps) || (w0 > 1-eps)

??????? J_t(i+1) = MAXD;??? % T = i

??? else

??????? miu1 = miuK / w0;

??????? miu2 = (meanT-miuK) / (1-w0);

??????? var1 = (((0 : i)'-miu1).^2)' * prob(1 : i+1);

??????? var1 = var1 / w0;? % variance

??????? var2 = (((i+1 : n)'-miu2).^2)' * prob(i+2 : n+1);

??????? var2 = var2 / (1-w0);

?? ?????if var1 > eps && var2 > eps?? % in case of var1=0 or var2 =0

??????????? J_t(i+1) = 1+w0 * log(var1)+(1-w0) * log(var2)-2*w0*log(w0)-2*(1-w0)*log(1-w0);

??????? else

??????????? J_t(i+1) = MAXD;

??????? end

??? end

end

minJ = min(J_t);

index = find(J_t == minJ);

th = mean(index);

th = (th-1)/n

imagBW = im2bw(imag, th);

?

% figure, imshow(imagBW), title('kittler binary');