## 基本全局閾值處理方法

1. 為全局閾值T選擇一個初始的估計值

2. 用T分割圖像，產生兩組像素：G1由大於T的像素組成，G2由小於T的像素組成

3. 對G1和G2的像素分別計算平均灰度值m1和m2

4. 計算新的閾值T = 1/2 * (m1 + m2)

5. 重復步驟2-4，直到連續迭代中的T值差小於一個預定義的參數ΔT

``` 1 void threshold(short** in_array, short** out_array, long height, long width, int delt_t)
2 {
3     double T = 0xff / 2.0;
4     double m1 = 0.0, m2 = 0.0;
5     int m1_num = 0, m2_num = 0;
6
7     while(dabs(T, 0.5*(m1 + m2)) > delt_t){
8         T = 0.5 * (m1 + m2);
9         m1 = 0.0;
10         m2 = 0.0;
11         m1_num = 0;
12         m2_num = 0;
13
14         for (int i = 0; i < height; i++){
15             for (int j = 0; j < width; j++){
16                 if (in_array[i][j] <= T){
17                     m1 += in_array[i][j];
18                     m1_num++;
19                 }
20                 else{
21                     m2 += in_array[i][j];
22                     m2_num++;
23                 }
24             }
25         }
26         if (m1_num != 0)
27             m1 /= m1_num;
28         if (m2_num != 0)
29             m2 /= m2_num;
30         printf("%lf\n", T);
31     }
32     for (int i = 0; i < height; i++){
33         for (int j = 0; j < width; j++){
34             if (in_array[i][j] <= T)
35                 out_array[i][j] = 0xff;
36             else
37                 out_array[i][j] = 0x00;
38         }
39     }
40 }```

從實驗結果看出，第二組閾值處理的效果並不好，因此考慮更優的算法實現

## Otsu方法進行最佳全局閾值處理

Otsu方法在類間方差最大的情況下是最佳的

### 算法執行流程

```  1 double dabs(double a, double b)
2 {
3     if (a < b)
4         return b-a;
5     else
6         return a-b;
7 }
8
9 void calculate_histogram(long height, long width, short **image, unsigned long histogram[])
10 {
11     short k;
12     for(int i=0; i < height; i++){
13         for(int j=0; j < width; j++){
14             k = image[i][j];
15             histogram[k] = histogram[k] + 1;
16         }
17     }
18 }
19
20 void calculate_pi(long height, long width, unsigned long histogram[], double pi[])
21 {
22     for (int i = 0; i < GRAY_LEVELS; ++i){
23         pi[i] = (double)histogram[i] / (double)(height * width);
24     }
25 }
26
27 double p1(int k, double pi[])
28 {
29     double sum = 0.0;
30
31     for (int i = 0; i <= k; i++){
32         sum += pi[i];
33     }
34
35     return sum;
36 }
37
38 double m(int k, double pi[])
39 {
40     double sum = 0.0;
41
42     for (int i = 0; i <= k; i++){
43         sum += i * pi[i];
44     }
45
46     return sum;
47 }
48
49 double calculate_sigma(int k, double mg, double pi[])
50 {
51     double p1k = p1(k, pi);
52     double mk = m(k, pi);
53
54     if (p1k < 1e-10 || (1 - p1k) < 1e-10)
55         return 0.0;
56     else
57         return pow(mg * p1k - mk, 2) / (p1k * (1 - p1k));
58 }
59
60 void otsu(short** in_array, short** out_array, long height, long width)
61 {
62     unsigned long histogram[GRAY_LEVELS] = {};
63     double pi[GRAY_LEVELS] = {};
64     double sigma[GRAY_LEVELS] = {};
65     double mg;
66     short max_count = 0;
67     short k = 0;
68     double max_value = 0.0;
69     double k_star;
70
71     calculate_histogram(height, width, in_array, histogram);
72     calculate_pi(height, width, histogram, pi);
73     mg = m(GRAY_LEVELS-1, pi);
74
75     for (int i = 0; i < GRAY_LEVELS; i++)
76         sigma[i] = calculate_sigma(i, mg, pi);
77
78     max_value = sigma[0];
79     max_count = 1;
80     k = 0;
81     for (int i = 1; i < GRAY_LEVELS; i++){
82         if (dabs(sigma[i], max_value) < 1e-10){
83             k += i;
84             max_count++;
85         }
86         else if (sigma[i] > max_value){
87             max_value = sigma[i];
88             max_count = 1;
89             k = i;
90         }
91     }
92     k_star = k / max_count;
93
94     printf("%lf\n", k_star);
95     for (int i = 0; i < height; i++){
96         for (int j = 0; j < width; j++){
97             if (in_array[i][j] <= k_star)
98                 out_array[i][j] = 0x00;
99             else
100                 out_array[i][j] = 0xff;
101         }
102     }
103 }```