2

I'm trying to implement the following Minimum Error Thresholding (By J. Kittler and J. Illingworth) method in MATLAB.

You may have a look at the PDF:

My code is:

function [ Level ] = MET( IMG )
%Maximum Error Thresholding By Kittler
%   Finding the Min of a cost function J in any possible thresholding. The
%   function output is the Optimal Thresholding.

for t = 0:255 % Assuming 8 bit image
    I1 = IMG;
    I1 = I1(I1 <= t);
    q1 = sum(hist(I1, 256));

    I2 = IMG;
    I2 = I2(I2 > t);
    q2 = sum(hist(I2, 256));

    % J is proportional to the Overlapping Area of the 2 assumed Gaussians
    J(t + 1) = 1 + 2 * (q1 * log(std(I1, 1)) + q2 * log(std(I2, 1)))...
        -2 * (q1 * log(q1) + q2 * log(q2));
end

[~, Level] = min(J);

%Level = (IMG <= Level);

end

I've tried it on the following image: Letters

Original size image.

The target is to extract a binary image of the letters (Hebrew Letters). I applied the code on sub blocks of the image (40 x 40). Yet I got results which are inferior to K-Means Clustering method.

Did I miss something? Anyone has a better idea?

Thanks.

P.S. Would anyone add "Adaptive-Thresholding" to the subject tags (I can't as I'm new).

Royi
  • 4,640
  • 6
  • 46
  • 64

2 Answers2

4

Thresholding is a rather tricky business. In the many years I've been thresholding images I have not found one single technique that always performs well, and I have come to distrust the claims of universally excellent performance in CS journals.

The maximum error thresholding method only works on nicely bimodal histogram (but it works well on those). Your image looks like signal and background may not be clearly separated enough for this thresholding method to work.

If you want to make sure that the code works fine, you could create a test program like this and check both whether you get good initial segmentation, as well as at what level of 'bimodality' the code breaks down.

Jonas
  • 74,690
  • 10
  • 137
  • 177
  • Which method would you suggest for this kind of image? Thanks. – Royi Jan 14 '10 at 17:15
  • 1
    Honestly? I don't know. I have never worked with this kind of images, and thus I have very little intuition for this problem. I have looked at the histogram for this, and I noticed that even though the background is horribly inhomogeneous, the foreground is mostly <0.12. In other words, your letters are fairly uniform in blackness, which means that the peak in the histogram corresponding to them is very pronounced, so you can simply cut the histogram at the first minimum after the highest peak and get a reasonable first approximation. – Jonas Jan 14 '10 at 22:22
  • This first approximation gives you a pretty good idea of which pixels belong to the background. You can then segment the background into the brighter and darker portions. Within each portion, you can then do hysteresis thresholding to improve your initial guess, i.e. determine a local threshold, so that pixels that are below that threshold and that are connected to any pixel of the initial guess (you may have to clean up the initial guess a bit) are kept, whereas pixels that make the threshold but are unconnected are discarded. – Jonas Jan 14 '10 at 22:28
  • Obviously, there are many, many other ways to segment this image. Just find one that robustly gives you a good starting point, and come up with some pre-processing and post-processing routines that allow you to apply your prior knowledge of the images you're analyzing. If you want to solve problems, as opposed to write a CS paper, you should not get caught up too much hunting for the perfect core algorithm. You want to write a program that does the job. – Jonas Jan 14 '10 at 22:32
2

I think your code is not fully correct. You use the absolute histogram of the image instead of the relative histogram which is used in the paper. In addition, your code is rather inefficient as it computes two histograms per possible threshold. I implemented the algorithm myself. Maybe, someone can make use of it:

function [ optimalThreshold, J ] = kittlerMinimimErrorThresholding( img )
%KITTLERMINIMIMERRORTHRESHOLDING Compute an optimal image threshold.
%   Computes the Minimum Error Threshold as described in
%   
%   'J. Kittler and J. Illingworth, "Minimum Error Thresholding," Pattern
%   Recognition 19, 41-47 (1986)'.
%   
%   The image 'img' is expected to have integer values from 0 to 255.
%   'optimalThreshold' holds the found threshold. 'J' holds the values of
%   the criterion function.

%Initialize the criterion function
J = Inf * ones(255, 1);

%Compute the relative histogram
histogram = double(histc(img(:), 0:255)) / size(img(:), 1);

%Walk through every possible threshold. However, T is interpreted
%differently than in the paper. It is interpreted as the lower boundary of
%the second class of pixels rather than the upper boundary of the first
%class. That is, an intensity of value T is treated as being in the same
%class as higher intensities rather than lower intensities.
for T = 1:255

    %Split the hostogram at the threshold T.
    histogram1 = histogram(1:T);
    histogram2 = histogram((T+1):end);

    %Compute the number of pixels in the two classes.
    P1 = sum(histogram1);
    P2 = sum(histogram2);

    %Only continue if both classes contain at least one pixel.
    if (P1 > 0) && (P2 > 0)

        %Compute the standard deviations of the classes.
        mean1 = sum(histogram1 .* (1:T)') / P1;
        mean2 = sum(histogram2 .* (1:(256-T))') / P2;
        sigma1 = sqrt(sum(histogram1 .* (((1:T)' - mean1) .^2) ) / P1);
        sigma2 = sqrt(sum(histogram2 .* (((1:(256-T))' - mean2) .^2) ) / P2);

        %Only compute the criterion function if both classes contain at
        %least two intensity values.
        if (sigma1 > 0) && (sigma2 > 0)

            %Compute the criterion function.
            J(T) = 1 + 2 * (P1 * log(sigma1) + P2 * log(sigma2)) ...
                     - 2 * (P1 * log(P1) + P2 * log(P2));

        end
    end

end

%Find the minimum of J.
[~, optimalThreshold] = min(J);
optimalThreshold = optimalThreshold - 0.5;
Kocki
  • 123
  • 5
  • I want to use your code (kittler Minimim Error Thresholding) for CT images. and typical Hounsfield Unit values are between -1000 : >1000. So, How Can I set the range of the threshold (`T`)? – Ellie Jan 01 '19 at 14:24