图像处理——（源）大律法阈值分割（threshold）函数编程实现




 1 #include <iostream>
 2 #include <opencv2/core.hpp>
 3 #include <opencv2/highgui.hpp>
 4 #include <opencv2/imgproc.hpp>
 5  
 6 int Otsu(cv::Mat& src, cv::Mat& dst, int thresh){
 7     const int Grayscale = 256;
 8     int graynum[Grayscale] = { 0 };
 9     int r = src.rows;
10     int c = src.cols;
11     for (int i = 0; i < r; ++i){
12         const uchar* ptr = src.ptr<uchar>(i);
13         for (int j = 0; j < c; ++j){        //直方图统计
14             graynum[ptr[j]]++;
15         }
16     }
17  
18     double P[Grayscale] = { 0 };   
19     double PK[Grayscale] = { 0 };
20     double MK[Grayscale] = { 0 };
21     double srcpixnum = r*c, sumtmpPK = 0, sumtmpMK = 0;
22     for (int i = 0; i < Grayscale; ++i){
23         P[i] = graynum[i] / srcpixnum;   //每个灰度级出现的概率
24         PK[i] = sumtmpPK + P[i];         //概率累计和 
25         sumtmpPK = PK[i];
26         MK[i] = sumtmpMK + i*P[i];       //灰度级的累加均值                                                                                                                                                                                                                                                                                                                                                                                                        
27         sumtmpMK = MK[i];
28     }
29     
30     //计算类间方差
31     double Var=0;
32     for (int k = 0; k < Grayscale; ++k){
33         if ((MK[Grayscale-1] * PK[k] - MK[k])*(MK[Grayscale-1] * PK[k] - MK[k]) / (PK[k] * (1 - PK[k])) > Var){
34             Var = (MK[Grayscale-1] * PK[k] - MK[k])*(MK[Grayscale-1] * PK[k] - MK[k]) / (PK[k] * (1 - PK[k]));
35             thresh = k;
36         }
37     }
38  
39     //阈值处理
40     src.copyTo(dst);
41     for (int i = 0; i < r; ++i){
42         uchar* ptr = dst.ptr<uchar>(i);
43         for (int j = 0; j < c; ++j){
44             if (ptr[j]> thresh)
45                 ptr[j] = 255;
46             else
47                 ptr[j] = 0;
48         }
49     }
50     return thresh;
51 }
52  
53  
54 int main(){
55     cv::Mat src = cv::imread("E://lena.jpg");
56     if (src.empty()){
57         return -1;
58     }
59     if (src.channels() > 1)
60         cv::cvtColor(src, src, CV_RGB2GRAY);
61  
62     cv::Mat dst,dst2;
63     int thresh=0;
64     double t2 = (double)cv::getTickCount();
65     thresh=Otsu(src , dst, thresh); //Otsu
66     std::cout << "Mythresh=" << thresh << std::endl;
67     t2 = (double)cv::getTickCount() - t2;
68     double time2 = (t2 *1000.) / ((double)cv::getTickFrequency());
69     std::cout << "my_process=" << time2 << " ms. " << std::endl << std::endl;
70     double  Otsu = 0;
71  
72     Otsu=cv::threshold(src, dst2, Otsu, 255, CV_THRESH_OTSU + CV_THRESH_BINARY);
73     std::cout << "OpenCVthresh=" << Otsu << std::endl;
74  
75     cv::namedWindow("src", CV_WINDOW_NORMAL);
76     cv::imshow("src", src);
77     cv::namedWindow("dst", CV_WINDOW_NORMAL);
78     cv::imshow("dst", dst);
79     cv::namedWindow("dst2", CV_WINDOW_NORMAL);
80     cv::imshow("dst2", dst2);
81     
82     cv::waitKey(0);
83 }
萍水相逢逢萍水，浮萍之水水浮萍!