【OpenCV】找圆方法（阈值分割：大律算法otsu）

【OpenCV】找圆方法（阈值分割：大律算法otsu）

2018-04-13 Taily老段 OpenCV学习交流

原文地址：https://blog.csdn.net/taily_duan/article/details/51217777

作者：Taily老段

 

opencv查找轮廓---cvFindContours && cvDrawCountours 用法及例子

 

 

 

1、找到轮廓；

2、每个连通域的均值中心；

3、求连通域半径（平均）；

4、相似度——最小半径/最大半径；

5、根据相似度阈值、半径阈值来判断是否是圆

[cpp] view plain copy

1. #include <iostream>    

2. #include <opencv2/imgproc/imgproc.hpp>  

3. #include <opencv2/core/core.hpp>    

4. #include <opencv2/highgui/highgui.hpp>    

5. using namespace cv;      

6. using namespace std;    

7.   

8.   

9. int main()    

10. {    

11.     Mat q_MatImage;  

12. Mat q_MatImageGray;  

13. Mat q_MatImageShow;  

14. Mat q_MatImageShow2;  

15. q_MatImage=imread("1.png");//读入一张图片  

16. q_MatImage.copyTo(q_MatImageShow);  

17. q_MatImage.copyTo(q_MatImageShow2);  

18. cvtColor(q_MatImage,q_MatImageGray,CV_RGB2GRAY);  

19.     double q_dEpsilon = 10E-9;  

20.     unsigned int q_iReturn=0;  

21.   

22.     int q_iX,q_iY,q_iWidth,q_iHeight;  

23.     q_iX=20;  

24.     q_iY=40;  

25.     q_iWidth=600;  

26.     q_iHeight=420;  

27.   

28.     double q_dThresholdSimilarity=60;  

29.     double q_dThresholdMin=35;  

30.     double q_dThresholdMax=75;  

31.   

32.     //      Rect q_RectROI = Rect(q_iX,q_iY,q_iWidth,q_iHeight);  

33.     //      Mat q_MatROI = q_MatImageGray(q_RectROI);  

34.     //    

35.     //      threshold(q_MatROI, q_MatROI, 150, 255, CV_THRESH_BINARY);  

36.     threshold(q_MatImageGray, q_MatImageGray, 150, 255, CV_THRESH_BINARY);  

37.   

38.     namedWindow("Test1");       //创建一个名为Test窗口  

39.     imshow("Test1",q_MatImageGray);         //窗口中显示图像  

40.   

41.     vector<vector<Point>> q_vPointContours;  

42.   

43.     //findContours(q_MatROI, q_vPointContours, CV_RETR_EXTERNAL, CV_CHAIN_APPROX_NONE,Point(q_iX,q_iY));  

44.     findContours(q_MatImageGray, q_vPointContours, CV_RETR_EXTERNAL, CV_CHAIN_APPROX_NONE,Point(0,0));  

45.   

46.     size_t q_iAmountContours = q_vPointContours.size();  

47.   

48.     for (size_t i = 0;i < q_iAmountContours; i++)  

49.     {  

50.         size_t q_perNum = q_vPointContours[i].size();  

51.         for (size_t j = 0;j < q_iAmountContours; j++)  

52.         {  

53.             circle( q_MatImageGray, q_vPointContours[i][j] ,3 , CV_RGB(0,255,0),1, 8, 3 );  

54.         }  

55.     }  

56.   

57.     namedWindow("findContours");  

58.     imshow("findContours",q_MatImageGray);    

59.   

60.     std::vector<cv::Point2f> q_vPointCentersContours(q_iAmountContours);  

61.     std::vector<double> q_vdRadiusesContours(q_iAmountContours);  

62.     std::vector<double> q_vdSimilarityContours(q_iAmountContours);  

63.     std::vector<bool> q_vbFlagCircles(q_iAmountContours);  

64.   

65.     std::vector<double> q_vdRadiusesContour;  

66.     double q_dRadiusMax,q_dRadiusMin;  

67.     double q_dSumX,q_dSumY;  

68.     size_t q_iAmountPoints;  

69.   

70.     for(size_t q_iCycleContours=0;q_iCycleContours<q_iAmountContours;q_iCycleContours++)  

71.     {  

72.         q_dSumX=0.0;  

73.         q_dSumY=0.0;  

74.         q_iAmountPoints=q_vPointContours[q_iCycleContours].size();  

75.         if(0>=q_iAmountPoints)  

76.         {  

77.             continue;  

78.         }  

79.         for(size_t q_iCyclePoints=0;q_iCyclePoints<q_iAmountPoints;q_iCyclePoints++)  

80.         {  

81.             q_dSumX+=q_vPointContours[q_iCycleContours].at(q_iCyclePoints).x;  

82.             q_dSumY+=q_vPointContours[q_iCycleContours].at(q_iCyclePoints).y;  

83.         }  

84.   

85.         q_vPointCentersContours[q_iCycleContours].x=(float)(q_dSumX/q_iAmountPoints);//均值中心点X</span>  

86.         q_vPointCentersContours[q_iCycleContours].y=(float)(q_dSumY/q_iAmountPoints);//均值中心点Y</span>  

87.   

88.   

89.         q_vdRadiusesContour.resize(q_iAmountPoints);  

90.         double q_dDifferenceX,q_dDifferenceY;  

91.         double q_dSumRadius=0.0;  

92.         q_dRadiusMax=0.0;  

93.         q_dRadiusMin=DBL_MAX;;  

94.         for(size_t q_iCyclePoints=0;q_iCyclePoints<q_iAmountPoints;q_iCyclePoints++)  

95.         {  

96.             q_dDifferenceX=q_vPointCentersContours[q_iCycleContours].x-q_vPointContours[q_iCycleContours].at(q_iCyclePoints).x;  

97.             q_dDifferenceY=q_vPointCentersContours[q_iCycleContours].y-q_vPointContours[q_iCycleContours].at(q_iCyclePoints).y;  

98.             q_vdRadiusesContour[q_iCyclePoints]=sqrt(q_dDifferenceX*q_dDifferenceX+q_dDifferenceY*q_dDifferenceY);  

99.   

100.             if(q_vdRadiusesContour[q_iCyclePoints]>q_dRadiusMax)  

101.             {  

102.                 q_dRadiusMax=q_vdRadiusesContour[q_iCyclePoints];  

103.             }  

104.             if(q_vdRadiusesContour[q_iCyclePoints]<q_dRadiusMin)  

105.             {  

106.                 q_dRadiusMin=q_vdRadiusesContour[q_iCyclePoints];  

107.             }  

108.   

109.             q_dSumRadius+=q_vdRadiusesContour[q_iCyclePoints];  

110.         }  

111.         q_vdRadiusesContours[q_iCycleContours]=q_dSumRadius/q_iAmountPoints;   //均值半径  

112.   

113.         q_vdSimilarityContours[q_iCycleContours]=100.0*q_dRadiusMin/q_dRadiusMax;  //相似度  

114.         if((q_dThresholdSimilarity<q_vdSimilarityContours[q_iCycleContours])&&  

115.             (q_dThresholdMin<q_vdRadiusesContours[q_iCycleContours])&&  

116.             (q_dThresholdMax>q_vdRadiusesContours[q_iCycleContours]))    //判断是否是圆  

117.         {  

118.             q_vbFlagCircles[q_iCycleContours]=true;  

119.         }  

120.         else  

121.         {  

122.             q_vbFlagCircles[q_iCycleContours]=false;  

123.         }  

124.     }  

125.   

126.   

127.     if(q_dEpsilon < 10)  

128.     {  

129.         cv::Point q_PointCenterCurrent;  

130.         for(size_t q_iCycleContours=0;q_iCycleContours<q_iAmountContours;q_iCycleContours++)  

131.         {  

132.             if(q_vbFlagCircles[q_iCycleContours])  

133.             {  

134.                 q_PointCenterCurrent.x=(int)(q_vPointCentersContours[q_iCycleContours].x);  

135.                 q_PointCenterCurrent.y=(int)(q_vPointCentersContours[q_iCycleContours].y);  

136.                 circle(q_MatImageShow,q_PointCenterCurrent,3,Scalar(0.0,0.0,255.0),0);  

137.             }  

138.         }  

139.     }  

140.   

141.     int q_iIndexResultBegin=4;  

142.     int q_iAmountCircleResult=4;  

143.     int q_iIndexCiecleCurrent;  

144.   

145.     int q_iCountCircles=0;  

146.   

147.     for(size_t q_iCycleContours=0;q_iCycleContours<q_iAmountContours;q_iCycleContours++)  

148.     {  

149.         if(q_vbFlagCircles[q_iCycleContours])  

150.         {  

151.             q_iIndexCiecleCurrent=q_iIndexResultBegin+q_iAmountCircleResult*q_iCountCircles;  

152.             //          match_result[q_iIndexCiecleCurrent]=(float)(q_vdSimilarityContours[q_iCycleContours]);  

153.             //          match_result[q_iIndexCiecleCurrent+1]=(float)(q_vdRadiusesContours[q_iCycleContours]);  

154.             //          match_result[q_iIndexCiecleCurrent+2]=(float)(q_vPointCentersContours[q_iCycleContours].x);  

155.             //          match_result[q_iIndexCiecleCurrent+3]=(float)(q_vPointCentersContours[q_iCycleContours].y);  

156.             q_iCountCircles++;  

157.         }  

158.     }  

159.     cout << "总共找到 " <<  q_iCountCircles << "个圆！" << endl;  

160.   

161.   

162.     namedWindow("Test");        //创建一个名为Test窗口  

163.     imshow("Test",q_MatImageShow);//窗口中显示图像  

164.     waitKey();              //等待5000ms后窗口自动关闭  

165. }  

 

 

 

大律算法otsu：

 

[cpp] view plain copy

1. int thresh = Otsu(q_MatImageGray);      

2. threshold(q_MatImageGray, q_MatImageGray, thresh, 255, CV_THRESH_BINARY);    

3.   

4. for(int i=0; i < q_MatImageGray.rows; i++)        

5. {     

6.     for(int j = 0; j < q_MatImageGray.cols; j++)         

7.     {     

8.         q_MatImageGray.at<uchar>(i,j) = 255 -q_MatImageGray.at<uchar>(i,j);  

9.     }          

10. }  

 

[cpp] view plain copy

1. int Otsu(Mat src)          

2. {  

3.     int height=src.rows;        

4.     int width =src.cols;              

5.   

6.   

7.     //histogram          

8.     float histogram[256] = {0};          

9.     for(int i=0; i < height; i++)        

10.     {       

11.         unsigned char* p=(unsigned char*)src.ptr<uchar>(i);   

12.         for(int j = 0; j < width; j++)         

13.         {     

14.             histogram[*p++]++;          

15.         }          

16.     }          

17.     //normalize histogram          

18.     int size = height * width;          

19.     for(int i = 0; i < 256; i++)        

20.     {          

21.         histogram[i] = histogram[i] / size;          

22.     }          

23.   

24.   

25.     //average pixel value          

26.     float avgValue=0;          

27.     for(int i=0; i < 256; i++)        

28.     {          

29.         avgValue += i * histogram[i];  //整幅图像的平均灰度        

30.     }           

31.   

32.   

33.     int threshold;            

34.     float maxVariance=0;          

35.     float w = 0, u = 0;          

36.     for(int i = 0; i < 256; i++)         

37.     {          

38.         w += histogram[i];  //假设当前灰度i为阈值, 0~i 灰度的像素(假设像素值在此范围的像素叫做前景像素) 所占整幅图像的比例        

39.         u += i * histogram[i];  // 灰度i 之前的像素(0~i)的平均灰度值： 前景像素的平均灰度值        

40.   

41.   

42.         float t = avgValue * w - u;          

43.         float variance = t * t / (w * (1 - w) );          

44.         if(variance > maxVariance)         

45.         {          

46.             maxVariance = variance;          

47.             threshold = i;          

48.         }          

49.     }          

50.   

51.   

52.     return threshold;          

53. }  

 

 

[cpp] view plain copy

1. int Otsu2(Mat src)  

2. {  

3.     int height=src.rows;        

4.     int width =src.cols;   

5.   

6.     int x=0,y=0;  

7.     int pixelCount[256];  

8.     float pixelPro[256];  

9.     int i, j, pixelSum = width * height, threshold = 0;  

10.   

11.     //初始化  

12.     for(i = 0; i < 256; i++)  

13.     {  

14.         pixelCount[i] = 0;  

15.         pixelPro[i] = 0;  

16.     }  

17.   

18.     //统计灰度级中每个像素在整幅图像中的个数  

19.     for(i = y; i < height; i++)  

20.     {  

21.         for(j = x;j <width;j++)  

22.         {  

23.             pixelCount[src.at<uchar>(i,j)]++;  

24.         }  

25.     }  

26.   

27.   

28.     //计算每个像素在整幅图像中的比例  

29.     for(i = 0; i < 256; i++)  

30.     {  

31.         pixelPro[i] = (float)(pixelCount[i]) / (float)(pixelSum);  

32.     }  

33.   

34.     //经典ostu算法,得到前景和背景的分割  

35.     //遍历灰度级[0,255],计算出方差最大的灰度值,为最佳阈值  

36.     float w0, w1, u0tmp, u1tmp, u0, u1, u,deltaTmp, deltaMax = 0;  

37.     for(i = 0; i < 256; i++)  

38.     {  

39.         w0 = w1 = u0tmp = u1tmp = u0 = u1 = u = deltaTmp = 0;  

40.   

41.         for(j = 0; j < 256; j++)  

42.         {  

43.             if(j <= i) //背景部分  

44.             {  

45.                 //以i为阈值分类，第一类总的概率  

46.                 w0 += pixelPro[j];        

47.                 u0tmp += j * pixelPro[j];  

48.             }  

49.             else       //前景部分  

50.             {  

51.                 //以i为阈值分类，第二类总的概率  

52.                 w1 += pixelPro[j];        

53.                 u1tmp += j * pixelPro[j];  

54.             }  

55.         }  

56.   

57.         u0 = u0tmp / w0;        //第一类的平均灰度  

58.         u1 = u1tmp / w1;        //第二类的平均灰度  

59.         u = u0tmp + u1tmp;      //整幅图像的平均灰度  

60.         //计算类间方差  

61.         deltaTmp = w0 * (u0 - u)*(u0 - u) + w1 * (u1 - u)*(u1 - u);  

62.         //找出最大类间方差以及对应的阈值  

63.         if(deltaTmp > deltaMax)  

64.         {     

65.             deltaMax = deltaTmp;  

66.             threshold = i;  

67.         }  

68.     }  

69.     //返回最佳阈值;  

70.     return threshold;  

71. }  

 

 

 

#include <iostream>    
#include <opencv2/imgproc/imgproc.hpp>  
#include <opencv2/core/core.hpp>    
#include <opencv2/highgui/highgui.hpp>    
using namespace cv;      
using namespace std;    
  
  
int main()    
{    
    Mat q_MatImage;  
    Mat q_MatImageGray;  
    Mat q_MatImageShow;  
    //Mat q_MatImageShow2;  
    //读入一张图片
    q_MatImage=imread("1.png");  
    q_MatImage.copyTo(q_MatImageShow);  
    //q_MatImage.copyTo(q_MatImageShow2);  
    cvtColor(q_MatImage,q_MatImageGray,CV_RGB2GRAY);  
    double q_dEpsilon = 10E-9;  
    unsigned int q_iReturn=0;  
  
    int q_iX,q_iY,q_iWidth,q_iHeight;  
    q_iX=20;  
    q_iY=40;  
    q_iWidth=600;  
    q_iHeight=420;  
  
    double q_dThresholdSimilarity=60;  
    double q_dThresholdMin=35;  
    double q_dThresholdMax=75;  

    namedWindow("Test1");       //创建一个名为Test窗口  
    imshow("Test1",q_MatImageGray);         //窗口中显示图像  
    waitKey();
  
    //      Rect q_RectROI = Rect(q_iX,q_iY,q_iWidth,q_iHeight);  
    //      Mat q_MatROI = q_MatImageGray(q_RectROI);  
    //    
    //      threshold(q_MatROI, q_MatROI, 150, 255, CV_THRESH_BINARY);  
    threshold(q_MatImageGray, q_MatImageGray, 150, 255, CV_THRESH_BINARY);  
    imshow("Test1",q_MatImageGray);         //窗口中显示图像  
    waitKey();
    vector<vector<Point>> q_vPointContours;  
  
    //findContours(q_MatROI, q_vPointContours, CV_RETR_EXTERNAL, CV_CHAIN_APPROX_NONE,Point(q_iX,q_iY));  
    findContours(q_MatImageGray, q_vPointContours, CV_RETR_EXTERNAL, CV_CHAIN_APPROX_NONE,Point(0,0));  
  
    size_t q_iAmountContours = q_vPointContours.size();  
  
    for (size_t i = 0;i < q_iAmountContours; i++)  
    {  
        size_t q_perNum = q_vPointContours[i].size();  
        for (size_t j = 0;j < q_iAmountContours; j++)  
        {  
            circle( q_MatImageGray, q_vPointContours[i][j] ,3 , CV_RGB(0,255,0),1, 8, 3 );  
        }  
    }  
  
    namedWindow("findContours");  
    imshow("findContours",q_MatImageGray);    
    waitKey();

    std::vector<cv::Point2f> q_vPointCentersContours(q_iAmountContours);  
    std::vector<double> q_vdRadiusesContours(q_iAmountContours);  
    std::vector<double> q_vdSimilarityContours(q_iAmountContours);  
    std::vector<bool> q_vbFlagCircles(q_iAmountContours);  
  
    std::vector<double> q_vdRadiusesContour;  
    double q_dRadiusMax,q_dRadiusMin;  
    double q_dSumX,q_dSumY;  
    size_t q_iAmountPoints;  
  
    for(size_t q_iCycleContours=0;q_iCycleContours<q_iAmountContours;q_iCycleContours++)  
    {  
        q_dSumX=0.0;  
        q_dSumY=0.0;  
        q_iAmountPoints=q_vPointContours[q_iCycleContours].size();  
        if(0>=q_iAmountPoints)  
        {  
            continue;  
        }  
        for(size_t q_iCyclePoints=0;q_iCyclePoints<q_iAmountPoints;q_iCyclePoints++)  
        {  
            q_dSumX+=q_vPointContours[q_iCycleContours].at(q_iCyclePoints).x;  
            q_dSumY+=q_vPointContours[q_iCycleContours].at(q_iCyclePoints).y;  
        }  
  
        q_vPointCentersContours[q_iCycleContours].x=(float)(q_dSumX/q_iAmountPoints);//均值中心点X</span>  
        q_vPointCentersContours[q_iCycleContours].y=(float)(q_dSumY/q_iAmountPoints);//均值中心点Y</span>  
  
  
        q_vdRadiusesContour.resize(q_iAmountPoints);  
        double q_dDifferenceX,q_dDifferenceY;  
        double q_dSumRadius=0.0;  
        q_dRadiusMax=0.0;  
        q_dRadiusMin=DBL_MAX;;  
        for(size_t q_iCyclePoints=0;q_iCyclePoints<q_iAmountPoints;q_iCyclePoints++)  
        {  
            q_dDifferenceX=q_vPointCentersContours[q_iCycleContours].x-q_vPointContours[q_iCycleContours].at(q_iCyclePoints).x;  
            q_dDifferenceY=q_vPointCentersContours[q_iCycleContours].y-q_vPointContours[q_iCycleContours].at(q_iCyclePoints).y;  
            q_vdRadiusesContour[q_iCyclePoints]=sqrt(q_dDifferenceX*q_dDifferenceX+q_dDifferenceY*q_dDifferenceY);  
  
            if(q_vdRadiusesContour[q_iCyclePoints]>q_dRadiusMax)  
            {  
                q_dRadiusMax=q_vdRadiusesContour[q_iCyclePoints];  
            }  
            if(q_vdRadiusesContour[q_iCyclePoints]<q_dRadiusMin)  
            {  
                q_dRadiusMin=q_vdRadiusesContour[q_iCyclePoints];  
            }  
  
            q_dSumRadius+=q_vdRadiusesContour[q_iCyclePoints];  
        }  
        q_vdRadiusesContours[q_iCycleContours]=q_dSumRadius/q_iAmountPoints;   //均值半径  
  
        q_vdSimilarityContours[q_iCycleContours]=100.0*q_dRadiusMin/q_dRadiusMax;  //相似度  
        if((q_dThresholdSimilarity<q_vdSimilarityContours[q_iCycleContours])&&  
            (q_dThresholdMin<q_vdRadiusesContours[q_iCycleContours])&&  
            (q_dThresholdMax>q_vdRadiusesContours[q_iCycleContours]))    //判断是否是圆  
        {  
            q_vbFlagCircles[q_iCycleContours]=true;  
        }  
        else  
        {  
            q_vbFlagCircles[q_iCycleContours]=false;  
        }  
    }  
  
    namedWindow("Test");        //创建一个名为Test窗口  
    int q_iCountCircles=0;
    if(q_dEpsilon < 10)  
    {  
        cv::Point q_PointCenterCurrent;  
        for(size_t q_iCycleContours=0;q_iCycleContours<q_iAmountContours;q_iCycleContours++)  
        {  
            if(q_vbFlagCircles[q_iCycleContours])  
            {  
                q_PointCenterCurrent.x=(int)(q_vPointCentersContours[q_iCycleContours].x);  
                q_PointCenterCurrent.y=(int)(q_vPointCentersContours[q_iCycleContours].y);  
                waitKey();
                circle(q_MatImageShow,q_PointCenterCurrent,3,Scalar(0.0,0.0,255.0),0); 
                imshow("Test",q_MatImageShow);//窗口中显示图像
                q_iCountCircles++; 
                cout << "共找到第 " <<  q_iCountCircles << "个圆！" << endl; 
            }  
        }  
    }  
  
    cout << "总共找到 " <<  q_iCountCircles << "个圆！" << endl;        
    waitKey();              //等待5000ms后窗口自动关闭 
    /*
    int q_iIndexResultBegin=4;  
    int q_iAmountCircleResult=4;  
    int q_iIndexCiecleCurrent;  
 
    for(size_t q_iCycleContours=0;q_iCycleContours<q_iAmountContours;q_iCycleContours++)  
    {  
        if(q_vbFlagCircles[q_iCycleContours])  
        {  
            q_iIndexCiecleCurrent=q_iIndexResultBegin+q_iAmountCircleResult*q_iCountCircles;  
            //          match_result[q_iIndexCiecleCurrent]=(float)(q_vdSimilarityContours[q_iCycleContours]);  
            //          match_result[q_iIndexCiecleCurrent+1]=(float)(q_vdRadiusesContours[q_iCycleContours]);  
            //          match_result[q_iIndexCiecleCurrent+2]=(float)(q_vPointCentersContours[q_iCycleContours].x);  
            //          match_result[q_iIndexCiecleCurrent+3]=(float)(q_vPointCentersContours[q_iCycleContours].y);  
            //q_iCountCircles++;  
        }  
    } 
    */
 
}