hough椭圆变换——opencv

参考：http://blog.csdn.net/u012507022/article/details/50979005
直接下载了参考文中的opencv代码，存在编译错误，对代码进行了一些小的修改得到如下：

主程序为：

#include "opencv2/core/core.hpp"  
#include "opencv2/highgui/highgui.hpp"  
#include "opencv2/imgproc/imgproc.hpp"  
#include "opencv2/features2d/features2d.hpp"    //需要添加该头文件  
#include <iostream>  
#include <vector>
#include "Math.h"
#include "ellipse.h"


using namespace cv;
using namespace std;

int main()
{

    Mat src = imread("C:\\Users\\aoe\\Desktop\\hough_ellipse\\hough_ellipse\\1.jpg");
    if (!src.data)
    {
        cout << "Read image error" << endl;
        return -1;
    }

    Mat dst;
    cvtColor(src, dst, CV_RGB2GRAY);
    GaussianBlur(dst, dst, Size(3, 3), 0, 0);       //高斯模糊（Gaussian Blur）
    Canny(dst, dst, 100, 200, 3);                  //Canny边缘检测
    namedWindow("Canny", CV_WINDOW_AUTOSIZE);
    imshow("Canny", dst);

    //提取轮廓************************************************
    vector<vector<Point>>  contours;
    vector<Vec4i> hierarchy;
    findContours(dst, contours, hierarchy, CV_RETR_EXTERNAL, CV_CHAIN_APPROX_NONE);
    //Hough变换************************************************
    int accumulate;   // Hough空间最大累积量
    Ellipse myellipse;
    Mat result(src.size(), CV_8UC3, Scalar(0));

    for (int i = 0; i< contours.size(); i++)
    {
        myellipse.Computer_axy(contours[i], dst.size());
        accumulate = myellipse.hough_ellipse(contours[i]);
        if (accumulate >= contours[i].size()*0.25)    // 判断是否超过给定阈值，判断是否为椭圆
            result = myellipse.draw_Eliipse(src);
        else
        {
            cout << "This profile is not an ellipse" << endl;
            break;
        }
    }

    namedWindow("Hough_result", CV_WINDOW_AUTOSIZE);
    imshow("Hough_result", result);

    waitKey();
    return 0;
}

ellipse.h为：

#include "opencv2/core/core.hpp"  
#include "opencv2/highgui/highgui.hpp"  
#include "opencv2/imgproc/imgproc.hpp"  
#include "opencv2/features2d/features2d.hpp"    //需要添加该头文件  
#include <iostream>  
#include <vector>


using namespace cv;
using namespace std;

class Ellipse
{
public:
    void Computer_axy(vector<Point> contour, Size imgsize);
    int hough_ellipse(vector<Point> contour);
    Mat draw_Eliipse(Mat);
private:
    double a;    //半长轴
    double b;    //短轴
    double theta;   //椭圆的旋转角度
    Point point_center;   //椭圆中心的坐标    
};

inline void Ellipse::Computer_axy(vector<Point> contour, Size imgsize)
{
    float Ly, Lx, LL;
    double maxVal;
    Mat distance(1, contour.size(), CV_32FC1);      //每一点到轮廓的所有距离
    Mat max_distance(imgsize, CV_32FC1, Scalar(0));  //每一点到轮廓的最大距离

    for (int i = 0; i<max_distance.rows; i++)
    {
        for (int j = 0; j<max_distance.cols; j++)
        {
            for (int n = 0; n<contour.size(); n++)
            {
                Ly = (i - contour.at(n).y)*(i - contour.at(n).y);
                Lx = (j - contour.at(n).x)*(j - contour.at(n).x);
                LL = sqrt(Ly + Lx);
                distance.at<float>(n) = LL;
            }
            minMaxLoc(distance, NULL, &maxVal, NULL, NULL);
            max_distance.at<float>(i, j) = maxVal;
        }
    }
    double minVal = 0; //最大值一定要赋初值，否则运行时会报错
    Point minLoc;
    minMaxLoc(max_distance, &minVal, NULL, &minLoc, NULL);
    a = minVal;
    point_center = minLoc;
}

inline int Ellipse::hough_ellipse(vector<Point> contour)
{
    double G, XX, YY;
    int B;
    Mat hough_space(floor(a + 1), 180, CV_8UC1, Scalar(0));  //高度:a，宽度180

    for (int k = 0; k<contour.size(); k++)
    {
        for (int w = 0; w<180; w++)
        {
            G = w*CV_PI / 180; //角度转换为弧度
            XX = pow(((contour.at(k).y - point_center.y)*cos(G) + (contour.at(k).x - point_center.x)*sin(G)), 2) / (a*a);
            YY = pow((-(contour.at(k).y - point_center.y)*sin(G) + (contour.at(k).x - point_center.x)*cos(G)), 2);

            B = floor(sqrt(abs(YY / (1 - XX))) + 1);
            if (B>0 && B <= a)
            {
                hough_space.at<uchar>(B, w) += 1;
            }
        }
    }
    double Circumference;
    double maxVal = 0; //最大值一定要赋初值，否则运行时会报错
    Point maxLoc;
    minMaxLoc(hough_space, NULL, &maxVal, NULL, &maxLoc);
    b = maxLoc.y;
    theta = maxLoc.x;
    Circumference = 2 * CV_PI*b + 4 * (a - b);

    return maxVal;
}

inline Mat Ellipse::draw_Eliipse(Mat src)
{
    cout << "长轴：" << a << endl;
    cout << "短轴：" << b << endl;
    cout << "椭圆中心：" << point_center << endl;
    cout << "theta：" << theta << endl;
    ellipse(src, point_center, Size(b, a), theta, 0, 360, Scalar(0, 255, 0), 3);
    return  src;
}