大津法阈值分割+直方图

  1 #include <iostream>
  2 #include <stdio.h>
  3 #include <math.h>
  4 
  5 using namespace std;
  6 
  7 int Histogram(unsigned short* a, int* hist,int h,int w)
  8 {
  9     if (a == NULL || hist == NULL){
 10         return -1;
 11     }
 12     for (int i = 0; i < pow(2,12); i++){
 13         hist[i] = 0;
 14     }
 15     for (int i = 0; i < h; i++){
 16         for (int j = 0; j < w; j++){
 17             hist[a[j + i*w]]++;
 18         }
 19     }
 20     return 0;
 21 }
 22 
 23 int OTSU(int* hist, int &thresh,int max)
 24 {
 25     if (hist == NULL)
 26         return -1;
 27     double sum0=0;
 28     double sum1=0;
 29     double n0=0;
 30     double n1=0;
 31     double u0=0;
 32     double u1=0;
 33     double w0=0;
 34     double w1=0;
 35     double u=0;
 36     double n = 0;
 37     double g = -1;
 38     double tempg = -1;
 39     for (int i = 0; i < max; i++){
 40         n+=hist[i];
 41     }
 42     for (int i = 0; i < max; i++){
 43         sum0 = 0;
 44         sum1 = 0;
 45         n0 += hist[i];
 46         n1 = n - n0;
 47         if (n1 == 0)
 48             break;
 49         w0 = n0 / n;
 50         w1 = 1 - w0;
 51         for (int j = 0; j <= i; j++){
 52             sum0 += j*hist[j];
 53         }
 54         u0 = sum0 / n0;
 55         for (int k = i + 1; k < max; k++){
 56             sum1 += k*hist[k];
 57         }
 58         u1 = sum1 / n1;
 59         g = w0*w1*(u0 - u1)*(u0 - u1);
 60         if (tempg < g)
 61         {
 62             tempg = g;
 63             thresh = i;
 64         }
 65     }
 66     return thresh;
 67 }
 68 
 69 int main()
 70 {
 71     int height=3072;
 72     int width=3072;
 73     int count;
 74     unsigned short* img = (unsigned short*)malloc(sizeof(unsigned short)*height*width);
 75     FILE* fp;
 76     fp = fopen("test_image.raw", "rb");
 77     if (fp == NULL){
 78         cout << "打开文件失败" << endl;
 79         return -1;
 80     }
 81     else{
 82         count = fread(img, sizeof(unsigned short), height*width, fp);
 83     }
 84     fclose(fp);
 85 
 86     int max = pow(2, 12);
 87 
 88     int* hist = (int*)malloc(max* sizeof(int));
 89 
 90     Histogram(img, hist, height, width);
 91 
 92 
 93     fp = fopen("hist.raw", "wb");
 94     fwrite(hist, sizeof(int), max, fp);
 95     fclose(fp);
 96 
 97     int thresh = 0;
 98     OTSU(hist, thresh, max);
 99     cout << thresh << endl;
100 
101     return 0;
102 }