# Thirteen histogram-based global binarization algorithms for images are described in terms of their principles, implementation, codes and effects.

The purpose of image binarization is to maximize the preservation of the parts of the image that are of interest. In many cases, it is also a necessary image preprocessing process before image analysis, feature extraction and pattern recognition.This seemingly simple problem has attracted wide attention from scholars at home and abroad over the past 40 years, and has generated hundreds of threshold selection methods. However, like other image segmentation algorithms, none of the existing methods can achieve satisfactory results for a variety of images.

Among these huge classification methods, the histogram-based global binary algorithm has an absolute market share. Each of these algorithms has its own implementation scheme from different scientific levels, and they all have some common features:

1. Simple;

2. The algorithm is easy to implement;

3. Fast execution.

Several of these methods are outlined in this paper.

1: Gray level value method:

1. Description: Even if the average gray level of the whole image is used as the threshold for binarization, this method can be used as the initial guess value of other methods.

2. Principle:

3. Implementation code:

``````public static int GetMeanThreshold(int[] HistGram)
{
int Sum = 0, Amount = 0;
for (int Y = 0; Y < 256; Y++)
{
Amount += HistGram[Y];
Sum += Y * HistGram[Y];
}
return Sum / Amount;
}``````

2. Percentage threshold (P-Tile method)

1. Description

Doyle's P-Tile, or P-Quantile, introduced in 1962, is arguably the oldest threshold selection method.This method sets a threshold based on a priori probability so that the binary target or background pixel scale equals a priori probability. This method is simple and efficient, but it has no power for images with a priori probability that is difficult to estimate.

2. The principle is simple and implemented directly in code.

``````/// <summary>
/// Percentage threshold
/// </summary>
/// <param name="HistGram">Histogram of Gray Images</param>
/// <param name="Tile">Percentage of the area of the background in the image</param>
/// <returns></returns>
public static int GetPTileThreshold(int[] HistGram, int Tile = 50)
{
int Y, Amount = 0, Sum = 0;
for (Y = 0; Y < 256; Y++) Amount += HistGram[Y];        //  Total number of pixels
for (Y = 0; Y < 256; Y++)
{
Sum = Sum + HistGram[Y];
if (Sum >= Amount * Tile / 100) return Y;
}
return -1;
}``````

3. Threshold based on valley floor minimum

1. Description:

This method is useful for images with distinct bimodal histograms that look for valley bottoms with bimodal peaks as thresholds, but it does not necessarily obtain thresholds. For those with flat histograms or unimodal images, this method is not appropriate.

2. Implementation process:

The implementation of this function is an iterative process. Before each processing, the histogram data is judged to see if it is already a bimodal histogram. If not, the histogram data is smoothed with a radius of 1 (window size of 3). If a certain number of iterations, such as 1000 times before a bimodal histogram is obtained, the function fails to execute.The final threshold value is the valley floor value between the two peaks.

Note that during encoding, smoothing requires information before the current pixel, so a backup of the data before smoothing is required.In addition, the accuracy of the first data type is limited, and the reshaped histogram data should not be applied, but must be converted to floating-point data for processing, or the correct results will not be obtained.

The algorithm reference papers are as follows:

J. M. S. Prewitt and M. L. Mendelsohn, "The analysis of cell images," innnals of the New York Academy of Sciences, vol. 128, pp. 1035-1053, 1966.
C. A. Glasbey, "An analysis of histogram-based thresholding algorithms," CVGIP: Graphical Models and Image Processing, vol. 55, pp. 532-537, 1993.

3. Implementation code:

``````public static int GetMinimumThreshold(int[] HistGram)
{
int Y, Iter = 0;
double[] HistGramC = new double[256];           // Floating-point numbers must be used to solve accuracy problems, otherwise correct results will not be obtained
double[] HistGramCC = new double[256];          // The averaging process destroys the previous data, so two pieces of data are required
for (Y = 0; Y < 256; Y++)
{
HistGramC[Y] = HistGram[Y];
HistGramCC[Y] = HistGram[Y];
}

// Smoothing the histogram by averaging three points
while (IsDimodal(HistGramCC) == false)                                        // Determine if it's already a bimodal image
{
HistGramCC[0] = (HistGramC[0] + HistGramC[0] + HistGramC[1]) / 3;                 // First point
for (Y = 1; Y < 255; Y++)
HistGramCC[Y] = (HistGramC[Y - 1] + HistGramC[Y] + HistGramC[Y + 1]) / 3;     // Middle Point
HistGramCC[255] = (HistGramC[254] + HistGramC[255] + HistGramC[255]) / 3;         // Last point
System.Buffer.BlockCopy(HistGramCC, 0, HistGramC, 0, 256 * sizeof(double));
Iter++;
if (Iter >= 1000) return -1;                                                   // Histogram cannot be smoothed to bimodal, error code returned
}
// The minimum value between the extreme two peaks of the threshold
bool Peakfound = false;
for (Y = 1; Y < 255; Y++)
{
if (HistGramCC[Y - 1] < HistGramCC[Y] && HistGramCC[Y + 1] < HistGramCC[Y]) Peakfound = true;
if (Peakfound == true && HistGramCC[Y - 1] >= HistGramCC[Y] && HistGramCC[Y + 1] >= HistGramCC[Y])
return Y - 1;
}
return -1;
}``````

The IsDimodal function is used to determine whether the histogram is bimodal, and the code is as follows:

``````    private static bool IsDimodal(double[] HistGram)       // Detect whether the histogram is bimodal
{
// The histogram's peaks are counted, and only the number of peaks, 2, is a bimodal
int Count = 0;
for (int Y = 1; Y < 255; Y++)
{
if (HistGram[Y - 1] < HistGram[Y] && HistGram[Y + 1] < HistGram[Y])
{
Count++;
if (Count > 2) return false;
}
}
if (Count == 2)
return true;
else
return false;
}``````

4. Effect:

Original Histogram * Smoothed Histogram

For this kind of image with obvious bimodal, the algorithm can still achieve good results.

4. Threshold based on bimodal mean

1. Description:

The algorithm is similar to the threshold method based on the minimum of the valley floor, except that the last step is not to get the valley floor value between the two peaks, but to take the average of the two peaks as the threshold value.

2. Reference code:

``````public static int GetIntermodesThreshold(int[] HistGram)
{
int Y, Iter = 0, Index;
double[] HistGramC = new double[256];           // Floating-point numbers must be used to solve accuracy problems, otherwise correct results will not be obtained
double[] HistGramCC = new double[256];          // The averaging process destroys the previous data, so two pieces of data are required
for (Y = 0; Y < 256; Y++)
{
HistGramC[Y] = HistGram[Y];
HistGramCC[Y] = HistGram[Y];
}
// Smoothing the histogram by averaging three points
while (IsDimodal(HistGramCC) == false)                                                  // Determine if it's already a bimodal image
{
HistGramCC[0] = (HistGramC[0] + HistGramC[0] + HistGramC[1]) / 3;                   // First point
for (Y = 1; Y < 255; Y++)
HistGramCC[Y] = (HistGramC[Y - 1] + HistGramC[Y] + HistGramC[Y + 1]) / 3;       // Middle Point
HistGramCC[255] = (HistGramC[254] + HistGramC[255] + HistGramC[255]) / 3;           // Last point
System.Buffer.BlockCopy(HistGramCC, 0, HistGramC, 0, 256 * sizeof(double));         // Back up your data in preparation for the next iteration
Iter++;
if (Iter >= 10000) return -1;                                                       // It seems that the histogram cannot be smoothed to bimodal, returning an error code
}
// Threshold is the average of two peaks
int[] Peak = new int[2];
for (Y = 1, Index = 0; Y < 255; Y++)
if (HistGramCC[Y - 1] < HistGramCC[Y] && HistGramCC[Y + 1] < HistGramCC[Y]) Peak[Index++] = Y - 1;
return ((Peak[0] + Peak[1]) / 2);
}``````

3. Effect:

Original Histogram * Smoothed Histogram

5. Optimal threshold for iteration

1. Description:

The algorithm first assumes a threshold value, then calculates the center value of foreground and background under that threshold. When the mean value of current and background center values is the same as the assumed threshold value, the iteration aborts and the value is used to binarize the threshold value.

2. Implementation process:

(1) The maximum and minimum gray values of the image are calculated and recorded as gl and gu, respectively, so that the initial threshold is:

(2) Segmenting the image into foreground and background according to threshold T0, and calculating average gray values Ab and Af:

(3) Order

If Tk=Tk+1, take Tk as the threshold value, otherwise, go to 2 to continue iteration.

3. Reference code:

``````public static int GetIterativeBestThreshold(int[] HistGram)
{
int X, Iter = 0;
int MeanValueOne, MeanValueTwo, SumOne, SumTwo, SumIntegralOne, SumIntegralTwo;
int MinValue, MaxValue;
int Threshold, NewThreshold;

for (MinValue = 0; MinValue < 256 && HistGram[MinValue] == 0; MinValue++) ;
for (MaxValue = 255; MaxValue > MinValue && HistGram[MinValue] == 0; MaxValue--) ;

if (MaxValue == MinValue) return MaxValue;          // There is only one color in the image
if (MinValue + 1 == MaxValue) return MinValue;      // There are only two colors in the image

Threshold = MinValue;
NewThreshold = (MaxValue + MinValue) >> 1;
while (Threshold != NewThreshold)    // End iteration when the current last two iterations have the same threshold
{
SumOne = 0; SumIntegralOne = 0;
SumTwo = 0; SumIntegralTwo = 0;
Threshold = NewThreshold;
for (X = MinValue; X <= Threshold; X++)         //The image is divided into target and background parts according to the threshold value, and the average gray value of the two parts is calculated.
{
SumIntegralOne += HistGram[X] * X;
SumOne += HistGram[X];
}
MeanValueOne = SumIntegralOne / SumOne;
for (X = Threshold + 1; X <= MaxValue; X++)
{
SumIntegralTwo += HistGram[X] * X;
SumTwo += HistGram[X];
}
MeanValueTwo = SumIntegralTwo / SumTwo;
NewThreshold = (MeanValueOne + MeanValueTwo) >> 1;       //Find a new threshold
Iter++;
if (Iter >= 1000) return -1;
}
return Threshold;
}``````

4. Effect:

Original Map * Binary Map * Histogram * Histogram

6. OSTU Law

1. Description:

This algorithm was proposed by Otsu, Japan, in 1979. The main idea is to take a threshold value to maximize the variance between foreground and background classes. graythresh in matlab is based on this algorithm.

2. Principle:

There are a lot of things on the network about how this algorithm works, and for the sake of space limitation, let's not go into details here.

3. Reference code:

``````    public static int GetOSTUThreshold(int[] HistGram)
{
int X, Y, Amount = 0;
int PixelBack = 0, PixelFore = 0, PixelIntegralBack = 0, PixelIntegralFore = 0, PixelIntegral = 0;
double OmegaBack, OmegaFore, MicroBack, MicroFore, SigmaB, Sigma;              // Interclass Variance;
int MinValue, MaxValue;
int Threshold = 0;

for (MinValue = 0; MinValue < 256 && HistGram[MinValue] == 0; MinValue++) ;
for (MaxValue = 255; MaxValue > MinValue && HistGram[MinValue] == 0; MaxValue--) ;
if (MaxValue == MinValue) return MaxValue;          // There is only one color in the image
if (MinValue + 1 == MaxValue) return MinValue;      // There are only two colors in the image

for (Y = MinValue; Y <= MaxValue; Y++) Amount += HistGram[Y];        //  Total number of pixels

PixelIntegral = 0;
for (Y = MinValue; Y <= MaxValue; Y++) PixelIntegral += HistGram[Y] * Y;
SigmaB = -1;
for (Y = MinValue; Y < MaxValue; Y++)
{
PixelBack = PixelBack + HistGram[Y];
PixelFore = Amount - PixelBack;
OmegaBack = (double)PixelBack / Amount;
OmegaFore = (double)PixelFore / Amount;
PixelIntegralBack += HistGram[Y] * Y;
PixelIntegralFore = PixelIntegral - PixelIntegralBack;
MicroBack = (double)PixelIntegralBack / PixelBack;
MicroFore = (double)PixelIntegralFore / PixelFore;
Sigma = OmegaBack * OmegaFore * (MicroBack - MicroFore) * (MicroBack - MicroFore);
if (Sigma > SigmaB)
{
SigmaB = Sigma;
Threshold = Y;
}
}
return Threshold;
}``````

4. Effect:

This algorithm has some adaptability to images with flat histograms.

7. 1-D Maximum Entropy

1. Description:

The algorithm introduces the concept of entropy in information theory into the image, and calculates the sum of the two parts of the entropy after the threshold segmentation to determine whether the threshold is the optimal threshold.

2. Algorithmic principles

There are also many articles in this area, leaving the reader to look for relevant information on his own.

3. Reference code:

``````public static int Get1DMaxEntropyThreshold(int[] HistGram)
{
int  X, Y,Amount=0;
double[] HistGramD = new double[256];
double SumIntegral, EntropyBack, EntropyFore, MaxEntropy;
int MinValue = 255, MaxValue = 0;
int Threshold = 0;

for (MinValue = 0; MinValue < 256 && HistGram[MinValue] == 0; MinValue++) ;
for (MaxValue = 255; MaxValue > MinValue && HistGram[MinValue] == 0; MaxValue--) ;
if (MaxValue == MinValue) return MaxValue;          // There is only one color in the image
if (MinValue + 1 == MaxValue) return MinValue;      // There are only two colors in the image

for (Y = MinValue; Y <= MaxValue; Y++) Amount += HistGram[Y];        //  Total number of pixels

for (Y = MinValue; Y <= MaxValue; Y++)   HistGramD[Y] = (double)HistGram[Y] / Amount+1e-17;

MaxEntropy = double.MinValue; ;
for (Y = MinValue + 1; Y < MaxValue; Y++)
{
SumIntegral = 0;
for (X = MinValue; X <= Y; X++) SumIntegral += HistGramD[X];
EntropyBack = 0;
for (X = MinValue; X <= Y; X++) EntropyBack += (-HistGramD[X] / SumIntegral * Math.Log(HistGramD[X] / SumIntegral));
EntropyFore = 0;
for (X = Y + 1; X <= MaxValue; X++) EntropyFore += (-HistGramD[X] / (1 - SumIntegral) * Math.Log(HistGramD[X] / (1 - SumIntegral)));
if (MaxEntropy < EntropyBack + EntropyFore)
{
Threshold = Y;
MaxEntropy = EntropyBack + EntropyFore;
}
}
return Threshold;
}``````

8. Torque Maintenance Method
1. Description:

The algorithm chooses an appropriate threshold value so that the binary image and the original gray image have three identical initial moment values.

2. Principle:

Reference paper: W. Tsai, "Moment-preserving thresholding: a new approach,"Comput.VisionGraphics Image Process., Vol. 29, PP. 377-393, 1985.

Since the paper cannot be downloaded (for a fee), just share the formulas found from other sources.

The function A\BC in this is visible in the code section.

3. Reference code:

``````public static byte GetMomentPreservingThreshold(int[] HistGram)
{
int X, Y, Index = 0, Amount=0;
double[] Avec = new double[256];
double X2, X1, X0, Min;

for (Y = 0; Y <= 255; Y++) Amount += HistGram[Y];        //  Total number of pixels
for (Y = 0; Y < 256; Y++) Avec[Y] = (double)A(HistGram, Y) / Amount;       // The threshold is chosen such that A(y,t)/A(y,n) is closest to x0.

// The following finds x0.

X2 = (double)(B(HistGram, 255) * C(HistGram, 255) - A(HistGram, 255) * D(HistGram, 255)) / (double)(A(HistGram, 255) * C(HistGram, 255) - B(HistGram, 255) * B(HistGram, 255));
X1 = (double)(B(HistGram, 255) * D(HistGram, 255) - C(HistGram, 255) * C(HistGram, 255)) / (double)(A(HistGram, 255) * C(HistGram, 255) - B(HistGram, 255) * B(HistGram, 255));
X0 = 0.5 - (B(HistGram, 255) / A(HistGram, 255) + X2 / 2) / Math.Sqrt(X2 * X2 - 4 * X1);

for (Y = 0, Min = double.MaxValue; Y < 256; Y++)
{
if (Math.Abs(Avec[Y] - X0) < Min)
{
Min = Math.Abs(Avec[Y] - X0);
Index = Y;
}
}
return (byte)Index;
}

private static double A(int[] HistGram, int Index)
{
double Sum = 0;
for (int Y = 0; Y <= Index; Y++)
Sum += HistGram[Y];
return Sum;
}

private static double B(int[] HistGram, int Index)
{
double Sum = 0;
for (int Y = 0; Y <= Index; Y++)
Sum += (double)Y * HistGram[Y];
return Sum;
}

private static double C(int[] HistGram, int Index)
{
double Sum = 0;
for (int Y = 0; Y <= Index; Y++)
Sum += (double)Y * Y * HistGram[Y];
return Sum;
}

private static double D(int[] HistGram, int Index)
{
double Sum = 0;
for (int Y = 0; Y <= Index; Y++)
Sum += (double)Y * Y * Y * HistGram[Y];
return Sum;
}``````

For many images, the algorithm can also achieve satisfactory results.

9. Threshold based on fuzzy set theory

The specific analysis of the algorithm is as follows: Principle, effect and code of an image binarization algorithm based on fuzzy set theory

This method also borrows the concept of Shannon's entropy, which generally achieves better segmentation results, whether for bimodal or unimodal images.

10. Kittler Minimum Error Classification

Due to limited effort, the following algorithms only give the algorithm's paper and related code.

The algorithm can be analyzed in detail as follows:

• Kittler, J & Illingworth, J (1986), "Minimum error thresholding", Pattern Recognition 19: 41-47

Reference code:

``````public static int GetKittlerMinError(int[] HistGram)
{
int X, Y;
int MinValue, MaxValue;
int Threshold ;
int PixelBack, PixelFore;
double OmegaBack, OmegaFore, MinSigma, Sigma, SigmaBack, SigmaFore;
for (MinValue = 0; MinValue < 256 && HistGram[MinValue] == 0; MinValue++) ;
for (MaxValue = 255; MaxValue > MinValue && HistGram[MinValue] == 0; MaxValue--) ;
if (MaxValue == MinValue) return MaxValue;          // There is only one color in the image
if (MinValue + 1 == MaxValue) return MinValue;      // There are only two colors in the image
Threshold = -1;
MinSigma = 1E+20;
for (Y = MinValue; Y < MaxValue; Y++)
{
PixelBack = 0; PixelFore = 0;
OmegaBack = 0; OmegaFore = 0;
for (X = MinValue; X <= Y; X++)
{
PixelBack += HistGram[X];
OmegaBack = OmegaBack + X * HistGram[X];
}
for (X = Y + 1; X <= MaxValue; X++)
{
PixelFore += HistGram[X];
OmegaFore = OmegaFore + X * HistGram[X];
}
OmegaBack = OmegaBack / PixelBack;
OmegaFore = OmegaFore / PixelFore;
SigmaBack = 0; SigmaFore = 0;
for (X = MinValue; X <= Y; X++) SigmaBack = SigmaBack + (X - OmegaBack) * (X - OmegaBack) * HistGram[X];
for (X = Y + 1; X <= MaxValue; X++) SigmaFore = SigmaFore + (X - OmegaFore) * (X - OmegaFore) * HistGram[X];
if (SigmaBack == 0 || SigmaFore == 0)
{
if (Threshold == -1)
Threshold = Y;
}
else
{
SigmaBack = Math.Sqrt(SigmaBack / PixelBack);
SigmaFore = Math.Sqrt(SigmaFore / PixelFore);
Sigma = 1 + 2 * (PixelBack * Math.Log(SigmaBack / PixelBack) + PixelFore * Math.Log(SigmaFore / PixelFore));
if (Sigma < MinSigma)
{
MinSigma = Sigma;
Threshold = Y;
}
}
}
return Threshold;
}``````

The algorithm does not work very well in practice.

Eleven: ISODATA (also called intermeans method)

Reference papers:

Reference code (uncluttered):

``````// Also called intermeans
// Iterative procedure based on the isodata algorithm [T.W. Ridler, S. Calvard, Picture
// thresholding using an iterative selection method, IEEE Trans. System, Man and
// Cybernetics, SMC-8 (1978) 630-632.]
// The procedure divides the image into objects and background by taking an initial threshold,
// then the averages of the pixels at or below the threshold and pixels above are computed.
// The averages of those two values are computed, the threshold is incremented and the
// process is repeated until the threshold is larger than the composite average. That is,
//  threshold = (average background + average objects)/2
// The code in ImageJ that implements this function is the getAutoThreshold() method in the ImageProcessor class.
//
// From: Tim Morris (dtm@ap.co.umist.ac.uk)
// Subject: Re: Thresholding method?
// posted to sci.image.processing on 1996/06/24
// The algorithm implemented in NIH Image sets the threshold as that grey
// value, G, for which the average of the averages of the grey values
// below and above G is equal to G. It does this by initialising G to the
// lowest sensible value and iterating:

// L = the average grey value of pixels with intensities < G
// H = the average grey value of pixels with intensities > G
// is G = (L + H)/2?
// yes => exit
// no => increment G and repeat
//
// There is a discrepancy with IJ because they are slightly different methods

public static int GetIsoDataThreshold(int[] HistGram)
{
int i, l, toth, totl, h, g = 0;
for (i = 1; i < HistGram.Length; i++)
{
if (HistGram[i] > 0)
{
g = i + 1;
break;
}
}
while (true)
{
l = 0;
totl = 0;
for (i = 0; i < g; i++)
{
totl = totl + HistGram[i];
l = l + (HistGram[i] * i);
}
h = 0;
toth = 0;
for (i = g + 1; i < HistGram.Length; i++)
{
toth += HistGram[i];
h += (HistGram[i] * i);
}
if (totl > 0 && toth > 0)
{
l /= totl;
h /= toth;
if (g == (int)Math.Round((l + h) / 2.0))
break;
}
g++;
if (g > HistGram.Length - 2)
{
return 0;
}
}
return g;
}``````

12. Shanbhag Method

Reference paper:

Shanbhag, Abhijit G. (1994), "Utilization of information measure as a means of image thresholding", Graph. Models Image Process. (Academic Press, Inc.) 56 (5): 414--419, ISSN 1049-9652, DOI 10.1006/cgip.1994.1037

Reference code (uncluttered):

``````public static int GetShanbhagThreshold(int[] HistGram)
{
int threshold;
int ih, it;
int first_bin;
int last_bin;
double term;
double tot_ent;  /* total entropy */
double min_ent;  /* max entropy */
double ent_back; /* entropy of the background pixels at a given threshold */
double ent_obj;  /* entropy of the object pixels at a given threshold */
double[] norm_histo = new double[HistGram.Length]; /* normalized histogram */
double[] P1 = new double[HistGram.Length]; /* cumulative normalized histogram */
double[] P2 = new double[HistGram.Length];

int total = 0;
for (ih = 0; ih < HistGram.Length; ih++)
total += HistGram[ih];

for (ih = 0; ih < HistGram.Length; ih++)
norm_histo[ih] = (double)HistGram[ih] / total;

P1[0] = norm_histo[0];
P2[0] = 1.0 - P1[0];
for (ih = 1; ih < HistGram.Length; ih++)
{
P1[ih] = P1[ih - 1] + norm_histo[ih];
P2[ih] = 1.0 - P1[ih];
}

/* Determine the first non-zero bin */
first_bin = 0;
for (ih = 0; ih < HistGram.Length; ih++)
{
if (!(Math.Abs(P1[ih]) < 2.220446049250313E-16))
{
first_bin = ih;
break;
}
}

/* Determine the last non-zero bin */
last_bin = HistGram.Length - 1;
for (ih = HistGram.Length - 1; ih >= first_bin; ih--)
{
if (!(Math.Abs(P2[ih]) < 2.220446049250313E-16))
{
last_bin = ih;
break;
}
}

// Calculate the total entropy each gray-level
// and find the threshold that maximizes it
threshold = -1;
min_ent = Double.MaxValue;

for (it = first_bin; it <= last_bin; it++)
{
/* Entropy of the background pixels */
ent_back = 0.0;
term = 0.5 / P1[it];
for (ih = 1; ih <= it; ih++)
{ //0+1?
ent_back -= norm_histo[ih] * Math.Log(1.0 - term * P1[ih - 1]);
}
ent_back *= term;

/* Entropy of the object pixels */
ent_obj = 0.0;
term = 0.5 / P2[it];
for (ih = it + 1; ih < HistGram.Length; ih++)
{
ent_obj -= norm_histo[ih] * Math.Log(1.0 - term * P2[ih]);
}
ent_obj *= term;

/* Total entropy */
tot_ent = Math.Abs(ent_back - ent_obj);

if (tot_ent < min_ent)
{
min_ent = tot_ent;
threshold = it;
}
}
return threshold;
}``````

Thirteen, Yen Law

Reference papers:

1) Yen J.C., Chang F.J., and Chang S. (1995) "A New Criterion  for Automatic Multilevel Thresholding" IEEE Trans. on Image  Processing, 4(3): 370-378
2) Sezgin M. and Sankur B. (2004) "Survey over Image Thresholding Techniques and Quantitative Performance Evaluation" Journal of  Electronic Imaging, 13(1): 146-165

Reference code (uncluttered):

``````// M. Emre Celebi
// 06.15.2007
// Ported to ImageJ plugin by G.Landini from E Celebi's fourier_0.8 routines
public static int GetYenThreshold(int[] HistGram)
{
int threshold;
int ih, it;
double crit;
double max_crit;
double[] norm_histo = new double[HistGram.Length]; /* normalized histogram */
double[] P1 = new double[HistGram.Length]; /* cumulative normalized histogram */
double[] P1_sq = new double[HistGram.Length];
double[] P2_sq = new double[HistGram.Length];

int total = 0;
for (ih = 0; ih < HistGram.Length; ih++)
total += HistGram[ih];

for (ih = 0; ih < HistGram.Length; ih++)
norm_histo[ih] = (double)HistGram[ih] / total;

P1[0] = norm_histo[0];
for (ih = 1; ih < HistGram.Length; ih++)
P1[ih] = P1[ih - 1] + norm_histo[ih];

P1_sq[0] = norm_histo[0] * norm_histo[0];
for (ih = 1; ih < HistGram.Length; ih++)
P1_sq[ih] = P1_sq[ih - 1] + norm_histo[ih] * norm_histo[ih];

P2_sq[HistGram.Length - 1] = 0.0;
for (ih = HistGram.Length - 2; ih >= 0; ih--)
P2_sq[ih] = P2_sq[ih + 1] + norm_histo[ih + 1] * norm_histo[ih + 1];

/* Find the threshold that maximizes the criterion */
threshold = -1;
max_crit = Double.MinValue;
for (it = 0; it < HistGram.Length; it++)
{
crit = -1.0 * ((P1_sq[it] * P2_sq[it]) > 0.0 ? Math.Log(P1_sq[it] * P2_sq[it]) : 0.0) + 2 * ((P1[it] * (1.0 - P1[it])) > 0.0 ? Math.Log(P1[it] * (1.0 - P1[it])) : 0.0);
if (crit > max_crit)
{
max_crit = crit;
threshold = it;
}
}
return threshold;
}``````

Many of the code in this line is taken from ImageJ, an open source software, for your reference: http://fiji.sc/wiki/index.php/Auto_Threshold Get more information here.

Finally, I have made a simple UI interface for these algorithms for interested readers.