c語言數字影象處理（十）：閾值處理
摘要：
定義
全域性閾值處理
假設某一副灰度圖有如下的直方圖，該影象由暗色背景下的較亮物體組成，從背景中提取這一物體時，將閾值T作為分割點，分割後的影象g(x, y)由下述公式給出，稱為全域性閾值處理
多閾值處理
本文僅完成全域性...
定義
全域性閾值處理
假設某一副灰度圖有如下的直方圖，該影象由暗色背景下的較亮物體組成，從背景中提取這一物體時，將閾值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 {
3double T = 0xff / 2.0;
4double m1 = 0.0, m2 = 0.0;
5int m1_num = 0, m2_num = 0;
6
7while(dabs(T, 0.5*(m1 + m2)) > delt_t){
8T = 0.5 * (m1 + m2);
9m1 = 0.0;
10m2 = 0.0;
11m1_num = 0;
12m2_num = 0;
13
14for (int i = 0; i < height; i++){
15for (int j = 0; j < width; j++){
16if (in_array[i][j] <= T){
17m1 += in_array[i][j];
18m1_num++;
19}
20else{
21m2 += in_array[i][j];
22m2_num++;
23}
24}
25}
26if (m1_num != 0)
27m1 /= m1_num;
28if (m2_num != 0)
29m2 /= m2_num;
30printf("%lf\n", T);
31}
32for (int i = 0; i < height; i++){
33for (int j = 0; j < width; j++){
34if (in_array[i][j] <= T)
35out_array[i][j] = 0xff;
36else
37out_array[i][j] = 0x00;
38}
39}
40 }
測試結果
從實驗結果看出，第二組閾值處理的效果並不好，因此考慮更優的演算法實現
Otsu方法進行最佳全域性閾值處理
閾值處理可視為一種統計決策理論問題，其目的是在把畫素分配給兩個或多個組的過程中引入的平均誤差最小。這一問題有個閉合形式的解，稱為貝葉斯決策規則。
Otsu方法在類間方差最大的情況下是最佳的
演算法執行流程
程式碼實現
1 double dabs(double a, double b)
2 {
3if (a < b)
4return b-a;
5else
6return a-b;
7 }
8
9 void calculate_histogram(long height, long width, short **image, unsigned long histogram[])
10 {
11short k;
12for(int i=0; i < height; i++){
13for(int j=0; j < width; j++){
14k = image[i][j];
15histogram[k] = histogram[k] + 1;
16}
17}
18 }
19
20 void calculate_pi(long height, long width, unsigned long histogram[], double pi[])
21 {
22for (int i = 0; i < GRAY_LEVELS; ++i){
23pi[i] = (double)histogram[i] / (double)(height * width);
24}
25 }
26
27 double p1(int k, double pi[])
28 {
29double sum = 0.0;
30
31for (int i = 0; i <= k; i++){
32sum += pi[i];
33}
34
35return sum;
36 }
37
38 double m(int k, double pi[])
39 {
40double sum = 0.0;
41
42for (int i = 0; i <= k; i++){
43sum += i * pi[i];
44}
45
46return sum;
47 }
48
49 double calculate_sigma(int k, double mg, double pi[])
50 {
51double p1k = p1(k, pi);
52double mk = m(k, pi);
53
54if (p1k < 1e-10 || (1 - p1k) < 1e-10)
55return 0.0;
56else
57return pow(mg * p1k - mk, 2) / (p1k * (1 - p1k));
58 }
59
60 void otsu(short** in_array, short** out_array, long height, long width)
61 {
62unsigned long histogram[GRAY_LEVELS] = {};
63double pi[GRAY_LEVELS] = {};
64double sigma[GRAY_LEVELS] = {};
65double mg;
66short max_count = 0;
67short k = 0;
68double max_value = 0.0;
69double k_star;
70
71calculate_histogram(height, width, in_array, histogram);
72calculate_pi(height, width, histogram, pi);
73mg = m(GRAY_LEVELS-1, pi);
74
75for (int i = 0; i < GRAY_LEVELS; i++)
76sigma[i] = calculate_sigma(i, mg, pi);
77
78max_value = sigma[0];
79max_count = 1;
80k = 0;
81for (int i = 1; i < GRAY_LEVELS; i++){
82if (dabs(sigma[i], max_value) < 1e-10){
83k += i;
84max_count++;
85}
86else if (sigma[i] > max_value){
87max_value = sigma[i];
88max_count = 1;
89k = i;
90}
91}
92k_star = k / max_count;
93
94printf("%lf\n", k_star);
95for (int i = 0; i < height; i++){
96for (int j = 0; j < width; j++){
97if (in_array[i][j] <= k_star)
98out_array[i][j] = 0x00;
99else
100out_array[i][j] = 0xff;
101}
102}
103 }
結果
