定義

全域性閾值處理

假設某一副灰度圖有如下的直方圖，該影象由暗色背景下的較亮物體組成，從背景中提取這一物體時，將閾值T作為分割點，分割後的影象g(x, y)由下述公式給出，稱為全域性閾值處理

多閾值處理

本文僅完成全域性閾值處理的演算法實現

基本全域性閾值處理方法

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 }

結果