c语言数字图像处理（十）：阈值处理

 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 }

  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 }