c語言數字圖像處理（十）：閾值處理

定義

全局閾值處理

假設某一副灰度圖有以下的直方圖，該圖像由暗色背景下的較亮物體組成，從背景中提取這一物體時，將閾值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 }

結果