Summarization and VC Implementation of Common Threshold Segmentation Methods for Gray Image

Classification: image processing OpenCV 2011-11-11 23:20 609 people read comment(11) Collection Report

In the field of image processing, binary images are widely used because of their low computational complexity and their ability to represent the key features of images. The common method to change gray image into binary image is to select a threshold value, and then to process each pixel of the image to be processed by a single point, that is, to compare its gray value with the threshold set, so as to obtain a binary black-and-white image. Because of its intuitiveness and easy implementation, this method has been in the central position in the field of image segmentation. This paper summarizes the thresholding image segmentation algorithms that the author has learned in recent years, describes the author's understanding of each algorithm, and implements these algorithms based on OpenCV and VC 6.0. Finally, the source code will be open, I hope you can make progress together. (The code in this article does not consider execution efficiency for the time being)

Firstly, the image to be segmented is given as follows:

1. Otsu Method (Maximum Interclass Variance Method)

This algorithm is a dynamic threshold segmentation algorithm proposed by Japanese Otsu. Its main idea is to divide the image into two parts: background and target according to the gray level characteristics, and to select the threshold value to maximize the variance between background and target. (The greater the variance between the background and the target, the greater the difference between the two parts. When part of the target is misclassified as background or part of the background is misclassified as target, the difference between the two parts will become smaller. Therefore, using the segmentation with the largest variance between classes means that the probability of misclassification is the smallest.) This is the main idea of this method. Its main realization principles are as follows:

1) Establishing gray histogram of image (there are L gray levels, each occurrence probability is p)

(2) Calculate the probability of background and target occurrence. The calculation method is as follows:

Suppose t is the selected threshold, A represents the background (gray level is 0~N). According to the elements in the histogram, Pa is the probability of background occurrence, B is the goal, and Pb is the probability of target occurrence.

3) Calculate the inter-class variance between regions A and B as follows:

 

The first expression calculates the average gray values of A and B regions respectively.

The second expression calculates the global average gray level of the gray image.

The third expression calculates the inter-class variance between regions A and B.

 

4) The inter-class variance of a single gray value is calculated in the above steps, so the optimal threshold value should be the gray value that can maximize the inter-class gray variance of A and B in the image. In the program, it is necessary to optimize the gray value of each occurrence.

My VC implementation code is as follows.

 

1. /***************************************************************************** 
2. * 
3. * \Function name: 
4. *   OneDimentionOtsu() 
5. * 
6. * \Input parameters: 
7. *   pGrayMat:      Binary image data 
8. *   width:         Graphic Dimension Width 
9. *   height:        Graphic Dimension Height 
10. *   nTlreshold:    Binary Segmentation Threshold after Algorithmic Processing 
11. * \Return value: 
12. *   nothing 
13. * \Function Description: Realize the Binary Segmentation of Gray Gray Gray Gray Image - Maximum Interclass Variance Method (Otsu algorithm, commonly known as Otsu algorithm) 
14. * 
15. ****************************************************************************/  
16.   
17. void CBinarizationDlg::OneDimentionOtsu(CvMat *pGrayMat, int width, int height, BYTE &nThreshold)  
18. {  
19.     double nHistogram[256];         //Gray histogram  
20.     double dVariance[256];          //Interclass variance  
21.     int N = height*width;           //Total Pixel Number  
22.     for(int i=0; i<256; i++)  
23.     {  
24.         nHistogram[i] = 0.0;  
25.         dVariance[i] = 0.0;  
26.     }  
27.     for(i=0; i<height; i++)  
28.     {  
29.         for(int j=0; j<width; j++)  
30.         {  
31.             unsigned char nData = (unsigned char)cvmGet(pGrayMat, i, j);  
32.             nHistogram[nData]++;     //Establishing Histogram  
33.         }  
34.     }  
35.     double Pa=0.0;      //Background occurrence probability  
36.     double Pb=0.0;      //Target occurrence probability  
37.     double Wa=0.0;      //Background average gray value  
38.     double Wb=0.0;      //Target average gray value  
39.     double W0=0.0;      //Global average gray value  
40.     double dData1=0.0,  dData2=0.0;  
41.     for(i=0; i<256; i++)     //Computing Global Average Gray Level  
42.     {  
43.         nHistogram[i] /= N;  
44.         W0 += i*nHistogram[i];  
45.     }  
46.     for(i=0; i<256; i++)     //Calculate the variance between classes for each gray value  
47.     {  
48.         Pa += nHistogram[i];  
49.         Pb = 1-Pa;  
50.         dData1 += i*nHistogram[i];  
51.         dData2 = W0-dData1;  
52.         Wa = dData1/Pa;  
53.         Wb = dData2/Pb;  
54.         dVariance[i] = (Pa*Pb* pow((Wb-Wa), 2));  
55.     }  
56.     //Traverse each variance to get the gray value corresponding to the maximum variance between classes  
57.     double temp=0.0;  
58.     for(i=0; i<256; i++)  
59.     {  
60.         if(dVariance[i]>temp)  
61.         {  
62.             temp = dVariance[i];  
63.             nThreshold = i;  
64.         }  
65.     }  
66. }  

The result of threshold segmentation is shown in the following figure, and the threshold value is 116.

 

2. One-dimensional cross-entropy method

 

Similar to the maximum variance between classes, this method was developed by Li and Lee applying the theory of entropy in information theory. Firstly, the concept of cross-entropy is briefly introduced.

For two distributions P and Q, the information cross-entropy D is defined as follows:

 

The physical meaning of this representation is the theoretical distance of information between two distributions. Another understanding is the change of information brought about by changing the distribution P to Q. For image segmentation, if the original image is to be replaced by the segmented image, the optimal segmentation basis should be to minimize the cross-entropy between the two images. The following is a brief summary of the process of minimum cross-entropy method.

It can be assumed that P is the gray distribution of the source image and Q is the gray distribution of the segmented image.

      

H is the statistical histogram in the formula above.

N is the total number of pixels in the image.

L is the total gray level of the source image.

P represents the source image, and each element represents the gray distribution (average gray value) at each gray level.

Q is the segmented binary image, and two u represent the average gray value of the two segmented regions respectively, where t is the threshold used for the segmented image.

According to the above definition, it is easy to deduce the quantitative expression of cross-entropy between P and Q based on the formula of cross-entropy with the sum of gray levels as the calculation amount.

According to the idea mentioned above, the minimum t of D is the optimal threshold in the sense of minimum cross-entropy.

The author's VC implementation code is as follows.

 

1. /***************************************************************************** 
2. * 
3. * \Function name: 
4. *   MiniCross() 
5. * 
6. * \Input parameters: 
7. *   pGrayMat:      Binary image data 
8. *   width:         Graphic Dimension Width 
9. *   height:        Graphic Dimension Height 
10. *   nTlreshold:    Binary Segmentation Threshold after Algorithmic Processing 
11. * \Return value: 
12. *   nothing 
13. * \Function Description: Realizing Binary Segmentation of Gray Gray Image-Minimum Cross Entropy Method 
14. * 
15. ****************************************************************************/  
16.   
17. void CBinarizationDlg::MiniCross(CvMat *pGrayMat, int width, int height, BYTE &nThreshold)  
18. {  
19.     double dHistogram[256];         //Gray histogram  
20.     double dEntropy[256];           //Cross Entropy of Each Pixel  
21.     int N = height*width;           //Total Pixel Number  
22.     for(int i=0; i<256; i++)  
23.     {  
24.         dHistogram[i] = 0.0;  
25.         dEntropy[i] = 0.0;  
26.     }  
27.     for(i=0; i<height; i++)  
28.     {  
29.         for(int j=0; j<width; j++)  
30.         {  
31.             unsigned char nData = (unsigned char)cvmGet(pGrayMat, i, j);  
32.             dHistogram[nData]++;     //Establishing Histogram  
33.         }  
34.     }  
35.     double Pa=0.0;      //Area 1 Average Gray Value  
36.     double Pb=0.0;      //Area 2 Average Gray Value  
37.     double P0=0.0;      //Global average gray value  
38.     double Wa=0.0;      //Part I Entropy  
39.     double Wb=0.0;      //Entropy of Part Two  
40.     double dData1=0.0, dData2=0.0;  //Intermediate value  
41.     double dData3=0.0, dData4=0.0;  //Intermediate value  
42.   
43.     for(i=0; i<256; i++)     //Computing Global Average Gray Level  
44.     {  
45.         dHistogram[i] /= N;  
46.         P0 += i*dHistogram[i];  
47.     }  
48.   
49.     for(i=0; i<256; i++)  
50.     {  
51.         Wa=Wb=dData1=dData2=dData3=dData4=Pa=Pb=0.0;  
52.         for(int j=0; j<256; j++)  
53.         {  
54.             if(j<=i)  
55.             {  
56.                 dData1 += dHistogram[j];  
57.                 dData2 += j*dHistogram[j];  
58.             }  
59.             else  
60.             {  
61.                 dData3 += dHistogram[j];  
62.                 dData4 += j*dHistogram[j];  
63.             }  
64.         }  
65.         Pa = dData2/dData1;  
66.         Pb = dData4/dData3;  
67.         for(j=0; j<256; j++)  
68.         {  
69.             if(j<=i)  
70.             {  
71.                 if((Pa!=0)&&(dHistogram[j]!=0))  
72.                 {  
73.                     double d1 = log(dHistogram[j]/Pa);  
74.                     Wa += j*dHistogram[j]*d1/log(2);  
75.                 }  
76.             }  
77.             else  
78.             {  
79.                 if((Pb!=0)&&(dHistogram[j]!=0))  
80.                 {  
81.                     double d2 = log(dHistogram[j]/Pb);  
82.                     Wb += j*dHistogram[j]*d2/log(2);  
83.                 }  
84.             }  
85.         }  
86.         dEntropy[i] = Wa+Wb;  
87.     }  
88.     //Ergodic Entropy Value and Gray Value of Minimum Cross Entropy  
89.     double temp=dEntropy[0];  
90.     for(i=1; i<256; i++)  
91.     {  
92.         if(dEntropy[i]<temp)  
93.         {  
94.             temp = dEntropy[i];  
95.             nThreshold = i;  
96.         }  
97.     }  
98. }  

The result of threshold segmentation is as follows: the threshold value obtained by solving the problem is 106.

 

 

3. Two-dimensional OTSU method

 

This method is an extension of the method of maximum variance between classes, which extends the maximum variance between two one-dimensional distributions to the maximum value of the trace of the matrix of interspecific dispersion, and adds the average pixel value of the neighborhood of the pixel point on the basis of considering the gray level of the pixel point.

Following is an analysis of the method's thinking and the pushing-down process according to my understanding:

1) First of all, we need to establish two-dimensional gray statistical histogram P(f, g);

If the gray level of the image is L-level, then the gray level of the average gray level of 8 neighborhoods of each pixel is L-level, and the histogram P is constructed accordingly. The horizontal axis of the two-dimensional statistical histogram is the gray value f(i, j) of each pixel, and the vertical coordinate is the average neighborhood value g(i, j) corresponding to the same point, where (0 < I < height, 0 < J < width, 0 < f(i, j) < L), while the gray value f of the corresponding P(f,g) image is the whole gray value of f, and the statistical value of the neighborhood gray mean value G is the proportion of the total pixels (i.e., the ratio of the neighborhood gray value of g). It is the joint probability density of gray value. Each element of P satisfies the following formula:

 

n is the statistic value of the gray value f in the whole image and the gray mean value g in the neighborhood.

N is the total number of pixels in the image.

2) For the two-dimensional statistical histogram shown in the following figure, t represents abscissa (gray value) and s represents ordinate (gray mean value of neighborhood of pixels)

 

For the threshold point (t,s) in an image, the difference between the gray value T and the gray mean s in its neighborhood should not be too large. If t is much larger than s (the point is located in the area II of the image above), it means that the gray value of the pixel is much larger than the gray mean of its neighborhood, so the point is likely to be a noise point, otherwise, if t is larger than S. It is much smaller (the point is located in the IV area on the way), that is, the pixel value of the point is much smaller than its neighborhood mean value, which means that it is an edge point. Therefore, we neglect these interference factors in background foreground segmentation, and consider that Pi,j=0 in these two regions. The remaining I and III regions represent the foreground and background respectively. Based on this, the optimal derivation of the discreteness criterion for the selected threshold (t,s) is derived.

3) Derivation of Discreteness Matrix Criteria at Threshold (t, s) Points

According to the analysis above, the probability of the occurrence of foreground and background segmented by threshold (t, s) is as follows:

Define two intermediate variables to facilitate the following derivation:

Accordingly, the gray mean vectors of these two parts can be deduced as follows (the two components are calculated according to the gray value and the gray mean of each point, respectively):

The gray mean vector of the whole image is:

In the same way as the one-dimensional Otsu method, we derive the variance between classes, which is two-dimensional, so we need to use the matrix form. Referring to the one-dimensional method, the "variance" matrix between classes can also be defined as follows:

In order to easily determine the "maximum" of such a matrix when it is implemented, the trace of the matrix (the sum of diagonals) is used to measure the "size" of the matrix in mathematics. Therefore, the trace of the matrix is used as the discreteness measure, and the derivation is as follows:

 

In this way, the optimal threshold combination can be obtained when the parameter is maximized (t,s).

The following is the implementation process of the algorithm:

1) Establishment of two-dimensional histogram

2) The histogram is traversed to calculate the matrix dispersion of each (t,s) combination, which is the so-called maximum inter-class variance in the one-dimensional Otsu method.

3) To get the maximum variance between classes (t,s). Since T represents the gray value and s represents the gray mean of the change point in its neighborhood, I think that s is the best choice when choosing the threshold, of course, t is also the best choice, because it can be seen from the solution results that this value is often very close.

The specific implementation code is as follows:

 

1. /***************************************************************************** 
2. * 
3. * \Function name: 
4. *   TwoDimentionOtsu() 
5. * 
6. * \Input parameters: 
7. *   pGrayMat:      Binary image data 
8. *   width:         Graphic Dimension Width 
9. *   height:        Graphic Dimension Height 
10. *   nTlreshold:    Binary Segmentation Threshold after Algorithmic Processing 
11. * \Return value: 
12. *   nothing 
13. * \Function Description: Realizing Binary Segmentation of Gray Gray Gray Image - Maximum Interclass Variance Method (Two-Dimensional Otsu Algorithms) 
14. * \Note: When constructing two-dimensional histogram, the 3*3 neighborhood mean of gray points is used. 
15. ******************************************************************************/  
16. void CBinarizationDlg::TwoDimentionOtsu(CvMat *pGrayMat, int width, int height, BYTE &nThreshold)  
17. {  
18.     double dHistogram[256][256];        //Establishment of two-dimensional gray histogram  
19.     double dTrMatrix = 0.0;             //Traces of Discrete Matrix  
20.     int N = height*width;               //Total Pixel Number  
21.     for(int i=0; i<256; i++)  
22.     {  
23.         for(int j=0; j<256; j++)  
24.             dHistogram[i][j] = 0.0;      //initialize variable  
25.     }  
26.     for(i=0; i<height; i++)  
27.     {  
28.         for(int j=0; j<width; j++)  
29.         {  
30.             unsigned char nData1 = (unsigned char)cvmGet(pGrayMat, i, j);  //Current gray value  
31.             unsigned char nData2 = 0;  
32.             int nData3 = 0;         //Note that the sum of nine values may exceed one byte  
33.             for(int m=i-1; m<=i+1; m++)  
34.             {  
35.                 for(int n=j-1; n<=j+1; n++)  
36.                 {  
37.                     if((m>=0)&&(m<height)&&(n>=0)&&(n<width))  
38.                         nData3 += (unsigned char)cvmGet(pGrayMat, m, n); //Current gray value  
39.                 }  
40.             }  
41.             nData2 = (unsigned char)(nData3/9);    //Zero-filling and Neighborhood Means for Overbound Index Values  
42.             dHistogram[nData1][nData2]++;  
43.         }  
44.     }  
45.     for(i=0; i<256; i++)  
46.         for(int j=0; j<256; j++)  
47.             dHistogram[i][j] /= N;  //Get the normalized probability distribution  
48.   
49.     double Pai = 0.0;      //Mean Vector i Component in Target Area  
50.     double Paj = 0.0;      //j component of mean vector in target area  
51.     double Pbi = 0.0;      //Mean Vector i Component in Background Area  
52.     double Pbj = 0.0;      //Background mean vector j component  
53.     double Pti = 0.0;      //Global mean vector i component  
54.     double Ptj = 0.0;      //Global mean vector j component  
55.     double W0 = 0.0;       //Joint probability density of target area  
56.     double W1 = 0.0;       //Joint probability density of background region  
57.     double dData1 = 0.0;  
58.     double dData2 = 0.0;  
59.     double dData3 = 0.0;  
60.     double dData4 = 0.0;   //Intermediate variable  
61.     int nThreshold_s = 0;  
62.     int nThreshold_t = 0;  
63.     double temp = 0.0;     //Seeking the Maximum  
64.     for(i=0; i<256; i++)  
65.     {  
66.         for(int j=0; j<256; j++)  
67.         {  
68.             Pti += i*dHistogram[i][j];  
69.             Ptj += j*dHistogram[i][j];  
70.         }  
71.     }  
72.     for(i=0; i<256; i++)  
73.     {  
74.         for(int j=0; j<256; j++)  
75.         {  
76.             W0 += dHistogram[i][j];  
77.             dData1 += i*dHistogram[i][j];  
78.             dData2 += j*dHistogram[i][j];  
79.   
80.             W1 = 1-W0;  
81.             dData3 = Pti-dData1;  
82.             dData4 = Ptj-dData2;  
83. /*          W1=dData3=dData4=0.0;   //Initialization of internal loop data 
84.             for(int s=i+1; s<256; s++) 
85.             { 
86.                 for(int t=j+1; t<256; t++) 
87.                 { 
88.                     W1 += dHistogram[s][t]; 
89.                     dData3 += s*dHistogram[s][t];  //Method 2 
90.                     dData4 += t*dHistogram[s][t];  //You can also add loops to your calculations. 
91.                 } 
92.             }*/  
93.   
94.             Pai = dData1/W0;  
95.             Paj = dData2/W0;  
96.             Pbi = dData3/W1;  
97.             Pbj = dData4/W1;   //Two mean vectors are obtained and expressed by four components.  
98.             dTrMatrix = ((W0*Pti-dData1)*(W0*Pti-dData1)+(W0*Ptj-dData1)*(W0*Ptj-dData2))/(W0*W1);  
99.             if(dTrMatrix > temp)  
100.             {  
101.                 temp = dTrMatrix;  
102.                 nThreshold_s = i;  
103.                 nThreshold_t = j;  
104.             }  
105.         }  
106.     }  
107.     nThreshold = nThreshold_t;   //Returns the gray value in the result  
108.     //nThreshold = 100;  
109. }  

 

The result of threshold segmentation is as follows: the threshold obtained by solving is 114. (s=114, t=117)

References

 

      [1] Nobuyuki Otsu. A Threshold SelectionMethod from Gray-Level Histograms.

      [2] C. H. Li and C. K. Lee. Minimum CrossEntropy Thresholding

      [3] Jian Gong, LiYuan Liand WeiNan Chen. Fast Recursive Algorithms For Two-Dimensional Thresholding.

 

Because the speed of the network is not enough, this article's formula and picture upload is slow. It's not easy to use the author for a long time. Therefore, if you want to reprint, please indicate the address of reprinting, in order to calculate the hard work of the author. At the same time, I hope to leave more messages and give the author points.