斑点SimpleBlobDetector检测



明：本文为博主原创文章，未经博主允许不得转载。



Opencv中提供了SimpleBlobDetector的特征点检测方法，正如它的名称，该算法使用最简单的方式来检测斑点类的特征点。下面我们就来分析一下该算法。

首先通过一系列连续的阈值把输入的灰度图像转换为一个二值图像的集合，阈值范围为[T₁, T₂]，步长为t，则所有阈值为：

T₁，T₁+t，T₁+2t，T₁+3t，……，T₂               （1）

第二步是利用Suzuki提出的算法通过检测每一幅二值图像的边界的方式提取出每一幅二值图像的连通区域，我们可以认为由边界所围成的不同的连通区域就是该二值图像的斑点；

第三步是根据所有二值图像斑点的中心坐标对二值图像斑点进行分类，从而形成灰度图像的斑点，属于一类的那些二值图像斑点最终形成灰度图像的斑点，具体来说就是，灰度图像的斑点是由中心坐标间的距离小于阈值T_b的那些二值图像斑点所组成的，即这些二值图像斑点属于该灰度图像斑点；

最后就是确定灰度图像斑点的信息——位置和尺寸。位置是属于该灰度图像斑点的所有二值图像斑点中心坐标的加权和，即公式2，权值q等于该二值图像斑点的惯性率的平方，它的含义是二值图像的斑点的形状越接近圆形，越是我们所希望的斑点，因此对灰度图像斑点位置的贡献就越大。尺寸则是属于该灰度图像斑点的所有二值图像斑点中面积大小居中的半径长度。

           （2）

在第二步中，并不是所有的二值图像的连通区域都可以认为是二值图像的斑点，我们往往通过一些限定条件来得到更准确的斑点。这些限定条件包括颜色，面积和形状，斑点的形状又可以用圆度，偏心率，或凸度来表示。

对于二值图像来说，只有两种斑点颜色——白色斑点和黑色斑点，我们只需要一种颜色的斑点，通过确定斑点的灰度值就可以区分出斑点的颜色。

连通区域的面积太大和太小都不是斑点，所以我们需要计算连通区域的面积，只有当该面积在我们所设定的最大面积和最小面积之间时，该连通区域才作为斑点被保留下来。

圆形的斑点是最理想的，任意形状的圆度C定义为：

C = 4pi(S/p^2)  （3）

其中，S和p分别表示该形状的面积和周长，当C为1时，表示该形状是一个完美的圆形，而当C为0时，表示该形状是一个逐渐拉长的多边形。

偏心率是指某一椭圆轨道与理想圆形的偏离程度，长椭圆轨道的偏心率高，而近于圆形的轨道的偏心率低。圆形的偏心率等于0，椭圆的偏心率介于0和1之间，而偏心率等于1表示的是抛物线。直接计算斑点的偏心率较为复杂，但利用图像矩的概念计算图形的惯性率，再由惯性率计算偏心率较为方便。偏心率E和惯性率I之间的关系为：

       E^2+I^2=1  （4）

因此圆形的惯性率等于1，惯性率越接近1，圆形的程度越高。

最后一个表示斑点形状的量是凸度。在平面中，凸形图指的是图形的所有部分都在由该图形切线所围成的区域的内部。我们可以用凸度来表示斑点凹凸的程度，凸度V的定义为：

V=S/H （5）

其中，H表示该斑点的凸壳面积

在计算斑点的面积，中心处的坐标，尤其是惯性率时，都可以应用图像矩的方法。由于公式较多，这里就不在详细讲解。

下面给出SimpleBlobDetector的源码分析。我们先来看看SimpleBlobDetector类的默认参数的设置：

[cpp] view plain copy print ?

1. SimpleBlobDetector::Params::Params()  
2. {  
3.     thresholdStep = 10;    //二值化的阈值步长，即公式1的t  
4.     minThreshold = 50;   //二值化的起始阈值，即公式1的T1  
5.     maxThreshold = 220;    //二值化的终止阈值，即公式1的T2  
6.     //重复的最小次数，只有属于灰度图像斑点的那些二值图像斑点数量大于该值时，该灰度图像斑点才被认为是特征点  
7.     minRepeatability = 2;     
8.     //最小的斑点距离，不同二值图像的斑点间距离小于该值时，被认为是同一个位置的斑点，否则是不同位置上的斑点  
9.     minDistBetweenBlobs = 10;  
10.   
11.     filterByColor = true;    //斑点颜色的限制变量  
12.     blobColor = 0;    //表示只提取黑色斑点；如果该变量为255，表示只提取白色斑点  
13.   
14.     filterByArea = true;    //斑点面积的限制变量  
15.     minArea = 25;    //斑点的最小面积  
16.     maxArea = 5000;    //斑点的最大面积  
17.   
18.     filterByCircularity = false;    //斑点圆度的限制变量，默认是不限制  
19.     minCircularity = 0.8f;    //斑点的最小圆度  
20.     //斑点的最大圆度，所能表示的float类型的最大值  
21.     maxCircularity = std::numeric_limits<float>::max();  
22.   
23.     filterByInertia = true;    //斑点惯性率的限制变量  
24.     //minInertiaRatio = 0.6;  
25.     minInertiaRatio = 0.1f;    //斑点的最小惯性率  
26.     maxInertiaRatio = std::numeric_limits<float>::max();    //斑点的最大惯性率  
27.   
28.     filterByConvexity = true;    //斑点凸度的限制变量  
29.     //minConvexity = 0.8;  
30.     minConvexity = 0.95f;    //斑点的最小凸度  
31.     maxConvexity = std::numeric_limits<float>::max();    //斑点的最大凸度  
32. }  

SimpleBlobDetector::Params::Params()
{
    thresholdStep = 10;    //二值化的阈值步长，即公式1的t
    minThreshold = 50;   //二值化的起始阈值，即公式1的T1
    maxThreshold = 220;    //二值化的终止阈值，即公式1的T2
    //重复的最小次数，只有属于灰度图像斑点的那些二值图像斑点数量大于该值时，该灰度图像斑点才被认为是特征点
    minRepeatability = 2;   
    //最小的斑点距离，不同二值图像的斑点间距离小于该值时，被认为是同一个位置的斑点，否则是不同位置上的斑点
    minDistBetweenBlobs = 10;

    filterByColor = true;    //斑点颜色的限制变量
    blobColor = 0;    //表示只提取黑色斑点；如果该变量为255，表示只提取白色斑点

    filterByArea = true;    //斑点面积的限制变量
    minArea = 25;    //斑点的最小面积
    maxArea = 5000;    //斑点的最大面积

    filterByCircularity = false;    //斑点圆度的限制变量，默认是不限制
    minCircularity = 0.8f;    //斑点的最小圆度
    //斑点的最大圆度，所能表示的float类型的最大值
    maxCircularity = std::numeric_limits<float>::max();

    filterByInertia = true;    //斑点惯性率的限制变量
    //minInertiaRatio = 0.6;
    minInertiaRatio = 0.1f;    //斑点的最小惯性率
    maxInertiaRatio = std::numeric_limits<float>::max();    //斑点的最大惯性率

    filterByConvexity = true;    //斑点凸度的限制变量
    //minConvexity = 0.8;
    minConvexity = 0.95f;    //斑点的最小凸度
    maxConvexity = std::numeric_limits<float>::max();    //斑点的最大凸度
}

我们再来介绍检测二值图像斑点的函数findBlobs。

[cpp] view plain copy print ?

1. //image为输入的灰度图像  
2. //binaryImage为二值图像  
3. //centers表示该二值图像的斑点  
4. void SimpleBlobDetector::findBlobs(const cv::Mat &image, const cv::Mat &binaryImage, vector<Center> ¢ers) const  
5. {  
6.     (void)image;  
7.     centers.clear();    //斑点变量清零  
8.   
9.     vector < vector<Point> > contours;    //定义二值图像的斑点的边界像素变量  
10.     Mat tmpBinaryImage = binaryImage.clone();    //复制二值图像  
11.     //调用findContours函数，找到当前二值图像的所有斑点的边界  
12.     findContours(tmpBinaryImage, contours, CV_RETR_LIST, CV_CHAIN_APPROX_NONE);  
13.   
14. #ifdef DEBUG_BLOB_DETECTOR  
15.     //  Mat keypointsImage;  
16.     //  cvtColor( binaryImage, keypointsImage, CV_GRAY2RGB );  
17.     //  
18.     //  Mat contoursImage;  
19.     //  cvtColor( binaryImage, contoursImage, CV_GRAY2RGB );  
20.     //  drawContours( contoursImage, contours, -1, Scalar(0,255,0) );  
21.     //  imshow("contours", contoursImage );  
22. #endif  
23.     //遍历当前二值图像的所有斑点  
24.     for (size_t contourIdx = 0; contourIdx < contours.size(); contourIdx++)  
25.     {  
26.         //结构类型Center代表着斑点，它包括斑点的中心位置，半径和权值  
27.         Center center;    //斑点变量  
28.         //初始化斑点中心的置信度，也就是该斑点的权值  
29.         center.confidence = 1;  
30.         //调用moments函数，得到当前斑点的矩  
31.         Moments moms = moments(Mat(contours[contourIdx]));  
32.         if (params.filterByArea)    //斑点面积的限制  
33.         {  
34.             double area = moms.m00;    //零阶矩即为二值图像的面积  
35.             //如果面积超出了设定的范围，则不再考虑该斑点  
36.             if (area < params.minArea || area >= params.maxArea)  
37.                 continue;  
38.         }  
39.   
40.         if (params.filterByCircularity)    //斑点圆度的限制  
41.         {  
42.             double area = moms.m00;    //得到斑点的面积  
43.             //计算斑点的周长  
44.             double perimeter = arcLength(Mat(contours[contourIdx]), true);  
45.             //由公式3得到斑点的圆度  
46.             double ratio = 4 * CV_PI * area / (perimeter * perimeter);  
47.             //如果圆度超出了设定的范围，则不再考虑该斑点  
48.             if (ratio < params.minCircularity || ratio >= params.maxCircularity)  
49.                 continue;  
50.         }  
51.   
52.         if (params.filterByInertia)    //斑点惯性率的限制  
53.         {  
54.             //计算公式13中最右侧等式中的开根号的值  
55.             double denominator = sqrt(pow(2 * moms.mu11, 2) + pow(moms.mu20 - moms.mu02, 2));  
56.             const double eps = 1e-2;    //定义一个极小值  
57.             double ratio;  
58.             if (denominator > eps)  
59.             {  
60.                 //cosmin和sinmin用于计算图像协方差矩阵中较小的那个特征值λ2  
61.                 double cosmin = (moms.mu20 - moms.mu02) / denominator;  
62.                 double sinmin = 2 * moms.mu11 / denominator;  
63.                 //cosmin和sinmin用于计算图像协方差矩阵中较大的那个特征值λ1  
64.                 double cosmax = -cosmin;  
65.                 double sinmax = -sinmin;  
66.                 //imin为λ2乘以零阶中心矩μ00  
67.                 double imin = 0.5 * (moms.mu20 + moms.mu02) - 0.5 * (moms.mu20 - moms.mu02) * cosmin - moms.mu11 * sinmin;  
68.                 //imax为λ1乘以零阶中心矩μ00  
69.                 double imax = 0.5 * (moms.mu20 + moms.mu02) - 0.5 * (moms.mu20 - moms.mu02) * cosmax - moms.mu11 * sinmax;  
70.                 ratio = imin / imax;    //得到斑点的惯性率  
71.             }  
72.             else  
73.             {  
74.                 ratio = 1;    //直接设置为1，即为圆  
75.             }  
76.             //如果惯性率超出了设定的范围，则不再考虑该斑点  
77.             if (ratio < params.minInertiaRatio || ratio >= params.maxInertiaRatio)  
78.                 continue;  
79.             //斑点中心的权值定义为惯性率的平方  
80.             center.confidence = ratio * ratio;  
81.         }  
82.   
83.         if (params.filterByConvexity)    //斑点凸度的限制  
84.         {  
85.             vector < Point > hull;    //定义凸壳变量  
86.             //调用convexHull函数，得到该斑点的凸壳  
87.             convexHull(Mat(contours[contourIdx]), hull);  
88.             //分别得到斑点和凸壳的面积，contourArea函数本质上也是求图像的零阶矩  
89.             double area = contourArea(Mat(contours[contourIdx]));  
90.             double hullArea = contourArea(Mat(hull));  
91.             double ratio = area / hullArea;    //公式5，计算斑点的凸度  
92.             //如果凸度超出了设定的范围，则不再考虑该斑点  
93.             if (ratio < params.minConvexity || ratio >= params.maxConvexity)  
94.                 continue;  
95.         }  
96.   
97.         //根据公式7，计算斑点的质心  
98.         center.location = Point2d(moms.m10 / moms.m00, moms.m01 / moms.m00);  
99.   
100.         if (params.filterByColor)    //斑点颜色的限制  
101.         {  
102.             //判断一下斑点的颜色是否正确  
103.             if (binaryImage.at<uchar> (cvRound(center.location.y), cvRound(center.location.x)) != params.blobColor)  
104.                 continue;  
105.         }  
106.   
107.         //compute blob radius  
108.         {  
109.             vector<double> dists;    //定义距离队列  
110.             //遍历该斑点边界上的所有像素  
111.             for (size_t pointIdx = 0; pointIdx < contours[contourIdx].size(); pointIdx++)  
112.             {  
113.                 Point2d pt = contours[contourIdx][pointIdx];    //得到边界像素坐标  
114.                 //计算该点坐标与斑点中心的距离，并放入距离队列中  
115.                 dists.push_back(norm(center.location - pt));  
116.             }  
117.             std::sort(dists.begin(), dists.end());    //距离队列排序  
118.             //计算斑点的半径，它等于距离队列中中间两个距离的平均值  
119.             center.radius = (dists[(dists.size() - 1) / 2] + dists[dists.size() / 2]) / 2.;  
120.         }  
121.   
122.         centers.push_back(center);    //把center变量压入centers队列中  
123.   
124. #ifdef DEBUG_BLOB_DETECTOR  
125.         //    circle( keypointsImage, center.location, 1, Scalar(0,0,255), 1 );  
126. #endif  
127.     }  
128. #ifdef DEBUG_BLOB_DETECTOR  
129.     //  imshow("bk", keypointsImage );  
130.     //  waitKey();  
131. #endif  
132. }  

//image为输入的灰度图像
//binaryImage为二值图像
//centers表示该二值图像的斑点
void SimpleBlobDetector::findBlobs(const cv::Mat &image, const cv::Mat &binaryImage, vector<Center> ¢ers) const
{
    (void)image;
    centers.clear();    //斑点变量清零

    vector < vector<Point> > contours;    //定义二值图像的斑点的边界像素变量
    Mat tmpBinaryImage = binaryImage.clone();    //复制二值图像
    //调用findContours函数，找到当前二值图像的所有斑点的边界
    findContours(tmpBinaryImage, contours, CV_RETR_LIST, CV_CHAIN_APPROX_NONE);

#ifdef DEBUG_BLOB_DETECTOR
    //  Mat keypointsImage;
    //  cvtColor( binaryImage, keypointsImage, CV_GRAY2RGB );
    //
    //  Mat contoursImage;
    //  cvtColor( binaryImage, contoursImage, CV_GRAY2RGB );
    //  drawContours( contoursImage, contours, -1, Scalar(0,255,0) );
    //  imshow("contours", contoursImage );
#endif
    //遍历当前二值图像的所有斑点
    for (size_t contourIdx = 0; contourIdx < contours.size(); contourIdx++)
    {
        //结构类型Center代表着斑点，它包括斑点的中心位置，半径和权值
        Center center;    //斑点变量
        //初始化斑点中心的置信度，也就是该斑点的权值
        center.confidence = 1;
        //调用moments函数，得到当前斑点的矩
        Moments moms = moments(Mat(contours[contourIdx]));
        if (params.filterByArea)    //斑点面积的限制
        {
            double area = moms.m00;    //零阶矩即为二值图像的面积
            //如果面积超出了设定的范围，则不再考虑该斑点
            if (area < params.minArea || area >= params.maxArea)
                continue;
        }

        if (params.filterByCircularity)    //斑点圆度的限制
        {
            double area = moms.m00;    //得到斑点的面积
            //计算斑点的周长
            double perimeter = arcLength(Mat(contours[contourIdx]), true);
            //由公式3得到斑点的圆度
            double ratio = 4 * CV_PI * area / (perimeter * perimeter);
            //如果圆度超出了设定的范围，则不再考虑该斑点
            if (ratio < params.minCircularity || ratio >= params.maxCircularity)
                continue;
        }

        if (params.filterByInertia)    //斑点惯性率的限制
        {
            //计算公式13中最右侧等式中的开根号的值
            double denominator = sqrt(pow(2 * moms.mu11, 2) + pow(moms.mu20 - moms.mu02, 2));
            const double eps = 1e-2;    //定义一个极小值
            double ratio;
            if (denominator > eps)
            {
                //cosmin和sinmin用于计算图像协方差矩阵中较小的那个特征值λ2
                double cosmin = (moms.mu20 - moms.mu02) / denominator;
                double sinmin = 2 * moms.mu11 / denominator;
                //cosmin和sinmin用于计算图像协方差矩阵中较大的那个特征值λ1
                double cosmax = -cosmin;
                double sinmax = -sinmin;
                //imin为λ2乘以零阶中心矩μ00
                double imin = 0.5 * (moms.mu20 + moms.mu02) - 0.5 * (moms.mu20 - moms.mu02) * cosmin - moms.mu11 * sinmin;
                //imax为λ1乘以零阶中心矩μ00
                double imax = 0.5 * (moms.mu20 + moms.mu02) - 0.5 * (moms.mu20 - moms.mu02) * cosmax - moms.mu11 * sinmax;
                ratio = imin / imax;    //得到斑点的惯性率
            }
            else
            {
                ratio = 1;    //直接设置为1，即为圆
            }
            //如果惯性率超出了设定的范围，则不再考虑该斑点
            if (ratio < params.minInertiaRatio || ratio >= params.maxInertiaRatio)
                continue;
            //斑点中心的权值定义为惯性率的平方
            center.confidence = ratio * ratio;
        }

        if (params.filterByConvexity)    //斑点凸度的限制
        {
            vector < Point > hull;    //定义凸壳变量
            //调用convexHull函数，得到该斑点的凸壳
            convexHull(Mat(contours[contourIdx]), hull);
            //分别得到斑点和凸壳的面积，contourArea函数本质上也是求图像的零阶矩
            double area = contourArea(Mat(contours[contourIdx]));
            double hullArea = contourArea(Mat(hull));
            double ratio = area / hullArea;    //公式5，计算斑点的凸度
            //如果凸度超出了设定的范围，则不再考虑该斑点
            if (ratio < params.minConvexity || ratio >= params.maxConvexity)
                continue;
        }

        //根据公式7，计算斑点的质心
        center.location = Point2d(moms.m10 / moms.m00, moms.m01 / moms.m00);

        if (params.filterByColor)    //斑点颜色的限制
        {
            //判断一下斑点的颜色是否正确
            if (binaryImage.at<uchar> (cvRound(center.location.y), cvRound(center.location.x)) != params.blobColor)
                continue;
        }

        //compute blob radius
        {
            vector<double> dists;    //定义距离队列
            //遍历该斑点边界上的所有像素
            for (size_t pointIdx = 0; pointIdx < contours[contourIdx].size(); pointIdx++)
            {
                Point2d pt = contours[contourIdx][pointIdx];    //得到边界像素坐标
                //计算该点坐标与斑点中心的距离，并放入距离队列中
                dists.push_back(norm(center.location - pt));
            }
            std::sort(dists.begin(), dists.end());    //距离队列排序
            //计算斑点的半径，它等于距离队列中中间两个距离的平均值
            center.radius = (dists[(dists.size() - 1) / 2] + dists[dists.size() / 2]) / 2.;
        }

        centers.push_back(center);    //把center变量压入centers队列中

#ifdef DEBUG_BLOB_DETECTOR
        //    circle( keypointsImage, center.location, 1, Scalar(0,0,255), 1 );
#endif
    }
#ifdef DEBUG_BLOB_DETECTOR
    //  imshow("bk", keypointsImage );
    //  waitKey();
#endif
}

最后介绍检测特征点的函数detectImpl。

[cpp] view plain copy print ?

1. void SimpleBlobDetector::detectImpl(const cv::Mat& image, std::vector<cv::KeyPoint>& keypoints, const cv::Mat&) const  
2. {  
3.     //TODO: support mask  
4.     keypoints.clear();    //特征点变量清零  
5.     Mat grayscaleImage;  
6.     //把彩色图像转换为二值图像  
7.     if (image.channels() == 3)  
8.         cvtColor(image, grayscaleImage, CV_BGR2GRAY);  
9.     else  
10.         grayscaleImage = image;  
11.     //二维数组centers表示所有得到的斑点，第一维数据表示的是灰度图像斑点，第二维数据表示的是属于该灰度图像斑点的所有二值图像斑点   
12.     //结构类型Center代表着斑点，它包括斑点的中心位置，半径和权值  
13.     vector < vector<Center> > centers;  
14.     //遍历所有阈值，进行二值化处理  
15.     for (double thresh = params.minThreshold; thresh < params.maxThreshold; thresh += params.thresholdStep)  
16.     {  
17.         Mat binarizedImage;  
18.         //调用threshold函数，把灰度图像grayscaleImage转换为二值图像binarizedImage  
19.         threshold(grayscaleImage, binarizedImage, thresh, 255, THRESH_BINARY);  
20.   
21. #ifdef DEBUG_BLOB_DETECTOR  
22.         //    Mat keypointsImage;  
23.         //    cvtColor( binarizedImage, keypointsImage, CV_GRAY2RGB );  
24. #endif  
25.         //变量curCenters表示该二值图像内的所有斑点  
26.         vector < Center > curCenters;  
27.         //调用findBlobs函数，对二值图像binarizedImage检测斑点，得到curCenters  
28.         findBlobs(grayscaleImage, binarizedImage, curCenters);  
29.         //newCenters表示在当前二值图像内检测到的不属于已有灰度图像斑点的那些二值图像斑点  
30.         vector < vector<Center> > newCenters;  
31.         //遍历该二值图像内的所有斑点  
32.         for (size_t i = 0; i < curCenters.size(); i++)  
33.         {  
34. #ifdef DEBUG_BLOB_DETECTOR  
35.             //      circle(keypointsImage, curCenters[i].location, curCenters[i].radius, Scalar(0,0,255),-1);  
36. #endif  
37.             // isNew表示的是当前二值图像斑点是否为新出现的斑点  
38.             bool isNew = true;  
39.             //遍历已有的所有灰度图像斑点，判断该二值图像斑点是否为新的灰度图像斑点  
40.             for (size_t j = 0; j < centers.size(); j++)  
41.             {  
42.                 //属于该灰度图像斑点的中间位置的那个二值图像斑点代表该灰度图像斑点，并把它的中心坐标与当前二值图像斑点的中心坐标比较，计算它们的距离dist  
43.                 double dist = norm(centers[j][ centers[j].size() / 2 ].location - curCenters[i].location);  
44.                 //如果距离大于所设的阈值，并且距离都大于上述那两个二值图像斑点的半径，则表示该二值图像的斑点是新的斑点  
45.                 isNew = dist >= params.minDistBetweenBlobs && dist >= centers[j][ centers[j].size() / 2 ].radius && dist >= curCenters[i].radius;  
46.                 //如果不是新的斑点，则需要把它添加到属于它的当前（即第j个）灰度图像的斑点中，因为通过上面的距离比较可知，当前二值图像斑点属于当前灰度图像斑点  
47.                 if (!isNew)  
48.                 {  
49.                     //把当前二值图像斑点存入当前（即第j个）灰度图像斑点数组的最后位置  
50.                     centers[j].push_back(curCenters[i]);  
51.                     //得到构成该灰度图像斑点的所有二值图像斑点的数量  
52.                     size_t k = centers[j].size() - 1;  
53.                     //按照半径由小至大的顺序，把新得到的当前二值图像斑点放入当前灰度图像斑点数组的适当位置处，由于二值化阈值是按照从小到大的顺序设置，所以二值图像斑点本身就是按照面积的大小顺序被检测到的，因此此处的排序处理要相对简单一些  
54.                     while( k > 0 && centers[j][k].radius < centers[j][k-1].radius )  
55.                     {  
56.                         centers[j][k] = centers[j][k-1];  
57.                         k--;  
58.                     }  
59.                     centers[j][k] = curCenters[i];  
60.                     //由于当前二值图像斑点已经找到了属于它的灰度图像斑点，因此退出for循环，无需再遍历已有的灰度图像斑点  
61.                     break;  
62.                 }  
63.             }  
64.             if (isNew)    //当前二值图像斑点是新的斑点  
65.             {  
66.                 //把当前斑点存入newCenters数组内  
67.                 newCenters.push_back(vector<Center> (1, curCenters[i]));  
68.                 //centers.push_back(vector<Center> (1, curCenters[i]));  
69.             }  
70.         }  
71.         //把当前二值图像内的所有newCenters复制到centers内  
72.         std::copy(newCenters.begin(), newCenters.end(), std::back_inserter(centers));  
73.   
74. #ifdef DEBUG_BLOB_DETECTOR  
75.         //    imshow("binarized", keypointsImage );  
76.         //waitKey();  
77. #endif  
78.     }    //所有二值图像斑点检测完毕  
79.     //遍历所有灰度图像斑点，得到特征点信息  
80.     for (size_t i = 0; i < centers.size(); i++)  
81.     {  
82.         //如果属于当前灰度图像斑点的二值图像斑点的数量小于所设阈值，则该灰度图像斑点不是特征点  
83.         if (centers[i].size() < params.minRepeatability)  
84.             continue;  
85.         Point2d sumPoint(0, 0);  
86.         double normalizer = 0;  
87.         //遍历属于当前灰度图像斑点的所有二值图像斑点  
88.         for (size_t j = 0; j < centers[i].size(); j++)  
89.         {  
90.             sumPoint += centers[i][j].confidence * centers[i][j].location;    //公式2的分子  
91.             normalizer += centers[i][j].confidence;    //公式2的分母  
92.         }  
93.         sumPoint *= (1. / normalizer);    //公式2，得到特征点的坐标位置  
94.         //保存该特征点的坐标位置和尺寸大小  
95.         KeyPoint kpt(sumPoint, (float)(centers[i][centers[i].size() / 2].radius));  
96.         keypoints.push_back(kpt);    //保存该特征点  
97.     }  
98.   
99. #ifdef DEBUG_BLOB_DETECTOR  
100.     namedWindow("keypoints", CV_WINDOW_NORMAL);  
101.     Mat outImg = image.clone();  
102.     for(size_t i=0; i<keypoints.size(); i++)  
103.     {  
104.         circle(outImg, keypoints[i].pt, keypoints[i].size, Scalar(255, 0, 255), -1);  
105.     }  
106.     //drawKeypoints(image, keypoints, outImg);  
107.     imshow("keypoints", outImg);  
108.     waitKey();  
109. #endif  
110. }  

void SimpleBlobDetector::detectImpl(const cv::Mat& image, std::vector<cv::KeyPoint>& keypoints, const cv::Mat&) const
{
    //TODO: support mask
    keypoints.clear();    //特征点变量清零
    Mat grayscaleImage;
    //把彩色图像转换为二值图像
    if (image.channels() == 3)
        cvtColor(image, grayscaleImage, CV_BGR2GRAY);
    else
        grayscaleImage = image;
    //二维数组centers表示所有得到的斑点，第一维数据表示的是灰度图像斑点，第二维数据表示的是属于该灰度图像斑点的所有二值图像斑点 
    //结构类型Center代表着斑点，它包括斑点的中心位置，半径和权值
    vector < vector<Center> > centers;
    //遍历所有阈值，进行二值化处理
    for (double thresh = params.minThreshold; thresh < params.maxThreshold; thresh += params.thresholdStep)
    {
        Mat binarizedImage;
        //调用threshold函数，把灰度图像grayscaleImage转换为二值图像binarizedImage
        threshold(grayscaleImage, binarizedImage, thresh, 255, THRESH_BINARY);

#ifdef DEBUG_BLOB_DETECTOR
        //    Mat keypointsImage;
        //    cvtColor( binarizedImage, keypointsImage, CV_GRAY2RGB );
#endif
        //变量curCenters表示该二值图像内的所有斑点
        vector < Center > curCenters;
        //调用findBlobs函数，对二值图像binarizedImage检测斑点，得到curCenters
        findBlobs(grayscaleImage, binarizedImage, curCenters);
        //newCenters表示在当前二值图像内检测到的不属于已有灰度图像斑点的那些二值图像斑点
        vector < vector<Center> > newCenters;
        //遍历该二值图像内的所有斑点
        for (size_t i = 0; i < curCenters.size(); i++)
        {
#ifdef DEBUG_BLOB_DETECTOR
            //      circle(keypointsImage, curCenters[i].location, curCenters[i].radius, Scalar(0,0,255),-1);
#endif
            // isNew表示的是当前二值图像斑点是否为新出现的斑点
            bool isNew = true;
            //遍历已有的所有灰度图像斑点，判断该二值图像斑点是否为新的灰度图像斑点
            for (size_t j = 0; j < centers.size(); j++)
            {
                //属于该灰度图像斑点的中间位置的那个二值图像斑点代表该灰度图像斑点，并把它的中心坐标与当前二值图像斑点的中心坐标比较，计算它们的距离dist
                double dist = norm(centers[j][ centers[j].size() / 2 ].location - curCenters[i].location);
                //如果距离大于所设的阈值，并且距离都大于上述那两个二值图像斑点的半径，则表示该二值图像的斑点是新的斑点
                isNew = dist >= params.minDistBetweenBlobs && dist >= centers[j][ centers[j].size() / 2 ].radius && dist >= curCenters[i].radius;
                //如果不是新的斑点，则需要把它添加到属于它的当前（即第j个）灰度图像的斑点中，因为通过上面的距离比较可知，当前二值图像斑点属于当前灰度图像斑点
                if (!isNew)
                {
                    //把当前二值图像斑点存入当前（即第j个）灰度图像斑点数组的最后位置
                    centers[j].push_back(curCenters[i]);
                    //得到构成该灰度图像斑点的所有二值图像斑点的数量
                    size_t k = centers[j].size() - 1;
                    //按照半径由小至大的顺序，把新得到的当前二值图像斑点放入当前灰度图像斑点数组的适当位置处，由于二值化阈值是按照从小到大的顺序设置，所以二值图像斑点本身就是按照面积的大小顺序被检测到的，因此此处的排序处理要相对简单一些
                    while( k > 0 && centers[j][k].radius < centers[j][k-1].radius )
                    {
                        centers[j][k] = centers[j][k-1];
                        k--;
                    }
                    centers[j][k] = curCenters[i];
                    //由于当前二值图像斑点已经找到了属于它的灰度图像斑点，因此退出for循环，无需再遍历已有的灰度图像斑点
                    break;
                }
            }
            if (isNew)    //当前二值图像斑点是新的斑点
            {
                //把当前斑点存入newCenters数组内
                newCenters.push_back(vector<Center> (1, curCenters[i]));
                //centers.push_back(vector<Center> (1, curCenters[i]));
            }
        }
        //把当前二值图像内的所有newCenters复制到centers内
        std::copy(newCenters.begin(), newCenters.end(), std::back_inserter(centers));

#ifdef DEBUG_BLOB_DETECTOR
        //    imshow("binarized", keypointsImage );
        //waitKey();
#endif
    }    //所有二值图像斑点检测完毕
    //遍历所有灰度图像斑点，得到特征点信息
    for (size_t i = 0; i < centers.size(); i++)
    {
        //如果属于当前灰度图像斑点的二值图像斑点的数量小于所设阈值，则该灰度图像斑点不是特征点
        if (centers[i].size() < params.minRepeatability)
            continue;
        Point2d sumPoint(0, 0);
        double normalizer = 0;
        //遍历属于当前灰度图像斑点的所有二值图像斑点
        for (size_t j = 0; j < centers[i].size(); j++)
        {
            sumPoint += centers[i][j].confidence * centers[i][j].location;    //公式2的分子
            normalizer += centers[i][j].confidence;    //公式2的分母
        }
        sumPoint *= (1. / normalizer);    //公式2，得到特征点的坐标位置
        //保存该特征点的坐标位置和尺寸大小
        KeyPoint kpt(sumPoint, (float)(centers[i][centers[i].size() / 2].radius));
        keypoints.push_back(kpt);    //保存该特征点
    }

#ifdef DEBUG_BLOB_DETECTOR
    namedWindow("keypoints", CV_WINDOW_NORMAL);
    Mat outImg = image.clone();
    for(size_t i=0; i<keypoints.size(); i++)
    {
        circle(outImg, keypoints[i].pt, keypoints[i].size, Scalar(255, 0, 255), -1);
    }
    //drawKeypoints(image, keypoints, outImg);
    imshow("keypoints", outImg);
    waitKey();
#endif
}

下面我们给出一个具体的应用实例。

[cpp] view plain copy print ?

1. #include "opencv2/core/core.hpp"  
2. #include "highgui.h"  
3. #include "opencv2/imgproc/imgproc.hpp"  
4. #include "opencv2/features2d/features2d.hpp"  
5. #include "opencv2/nonfree/nonfree.hpp"  
6.   
7. using namespace cv;  
8. //using namespace std;  
9.   
10. int main(int argc, char** argv)  
11. {  
12.    Mat img = imread("box_in_scene.png");  
13.   
14.    SimpleBlobDetector::Params params;    
15.    params.minThreshold = 40;    
16.    params.maxThreshold = 160;    
17.    params.thresholdStep = 5;    
18.    params.minArea = 100;    
19.    params.minConvexity = .05f;    
20.    params.minInertiaRatio = .05f;    
21.    params.maxArea = 8000;    
22.   
23.    SimpleBlobDetector detector(params);  
24.   
25.    vector<KeyPoint> key_points;  
26.   
27.    detector.detect(img,key_points);  
28.      
29.    Mat output_img;  
30.   
31.    drawKeypoints( img, key_points, output_img, Scalar(0,0,255), DrawMatchesFlags::DRAW_RICH_KEYPOINTS );  
32.   
33.    namedWindow("SimpleBlobDetector");  
34.    imshow("SimpleBlobDetector", output_img);  
35.    waitKey(0);  
36.   
37.    return 0;  
38. }  

#include "opencv2/core/core.hpp"
#include "highgui.h"
#include "opencv2/imgproc/imgproc.hpp"
#include "opencv2/features2d/features2d.hpp"
#include "opencv2/nonfree/nonfree.hpp"

using namespace cv;
//using namespace std;

int main(int argc, char** argv)
{
   Mat img = imread("box_in_scene.png");

   SimpleBlobDetector::Params params;  
   params.minThreshold = 40;  
   params.maxThreshold = 160;  
   params.thresholdStep = 5;  
   params.minArea = 100;  
   params.minConvexity = .05f;  
   params.minInertiaRatio = .05f;  
   params.maxArea = 8000;  

   SimpleBlobDetector detector(params);

   vector<KeyPoint> key_points;

   detector.detect(img,key_points);
   
   Mat output_img;

   drawKeypoints( img, key_points, output_img, Scalar(0,0,255), DrawMatchesFlags::DRAW_RICH_KEYPOINTS );

   namedWindow("SimpleBlobDetector");
   imshow("SimpleBlobDetector", output_img);
   waitKey(0);

   return 0;
}

在上述程序中，为了更好的检测到图像的特征点，我们修改了一些默认参数，如果仅仅使用系统的默认参数，则在实例化SimpleBlobDetector时，只需下列代码即可：

SimpleBlobDetector detector;

Opencv2.4.9源码分析DenseFeatureDetector

http://blog.youkuaiyun.com/zhaocj/article/details/45198965

DenseFeatureDetector可以生成具有一定密度和规律分布的图像特征点，它主要用于3D VIZ。

DenseFeatureDetector的原理是，把输入图像分割成大小相等的网格，每一个网格提取一个像素作为特征点。类似于图像尺度金字塔，该方法也可以生成不同层图像的特征点，每一层图像所分割的网格大小是不同的，即表示各层的尺度不同。

 

下面我们就来分析它的源码。

DenseFeatureDetector类的构造函数：

[cpp] view plain copy print ?

1. DenseFeatureDetector::DenseFeatureDetector( float _initFeatureScale, int _featureScaleLevels,  
2.                                       float _featureScaleMul, int _initXyStep,  
3.                                       int _initImgBound, bool _varyXyStepWithScale,  
4.                                       bool _varyImgBoundWithScale ) :  
5.     initFeatureScale(_initFeatureScale), featureScaleLevels(_featureScaleLevels),  
6.     featureScaleMul(_featureScaleMul), initXyStep(_initXyStep), initImgBound(_initImgBound),  
7.     varyXyStepWithScale(_varyXyStepWithScale), varyImgBoundWithScale(_varyImgBoundWithScale)  
8. {}  

DenseFeatureDetector::DenseFeatureDetector( float _initFeatureScale, int _featureScaleLevels,
                                      float _featureScaleMul, int _initXyStep,
                                      int _initImgBound, bool _varyXyStepWithScale,
                                      bool _varyImgBoundWithScale ) :
    initFeatureScale(_initFeatureScale), featureScaleLevels(_featureScaleLevels),
    featureScaleMul(_featureScaleMul), initXyStep(_initXyStep), initImgBound(_initImgBound),
    varyXyStepWithScale(_varyXyStepWithScale), varyImgBoundWithScale(_varyImgBoundWithScale)
{}

initFeatureScale表示初始图像层特征点的尺度，默认为1

featureScaleLevels表示需要构建多少层图像，默认为1

featureScaleMul表示各层图像之间参数的比例系数，该系数等于相邻两层图像之间的网格宽度之比，尺度之比，以及预留边界宽度之比，默认为0.1

initXyStep表示初始图像层的网格宽度，默认为6

initImgBound表示初始图像层的预留边界宽度，默认为0

varyXyStepWithScale表示各层图像是否进行网格宽度的调整，如果为false，则表示各层图像网格宽度都是initXyStep，如果为true，则表示各层图像网格宽度不等，它们之间的比例系数为featureScaleMul，默认为true

varyImgBoundWithScale表示各层图像是否进行预留边界宽度的调整，如果为false，则表示各层图像预留边界宽度都是initImgBound，如果为true，则表示各层图像预留边界宽度不等，它们之间的比例系数为featureScaleMul，默认为false

下面是检测特征点函数detectImpl：

[cpp] view plain copy print ?

1. void DenseFeatureDetector::detectImpl( const Mat& image, vector<KeyPoint>& keypoints, const Mat& mask ) const  
2. {  
3.     // curScale表示当前层图像特征点的尺度  
4.     float curScale = static_cast<float>(initFeatureScale);  
5.     //curStep表示当前层图像的网格宽度  
6.     int curStep = initXyStep;  
7.     //curBound表示当前层图像的预留边界宽度  
8.     int curBound = initImgBound;  
9.     //遍历各层图像  
10.     for( int curLevel = 0; curLevel < featureScaleLevels; curLevel++ )  
11.     {  
12.         //遍历当前层图像的所有网格，图像四周的预留边界处是没有网格的，横、纵坐标的步长就是网格的宽度  
13.         for( int x = curBound; x < image.cols - curBound; x += curStep )  
14.         {  
15.             for( int y = curBound; y < image.rows - curBound; y += curStep )  
16.             {  
17.                 //把网格的左上角坐标处的像素作为该网格的特征点，并保存  
18.                 keypoints.push_back( KeyPoint(static_cast<float>(x), static_cast<float>(y), curScale) );  
19.             }  
20.         }  
21.         //调整下一层图像特征点的尺度  
22.         curScale = static_cast<float>(curScale * featureScaleMul);  
23.         //如果varyXyStepWithScale为true，则调整下一层图像的网格宽度  
24.         if( varyXyStepWithScale ) curStep = static_cast<int>( curStep * featureScaleMul + 0.5f );  
25.         //如果varyImgBoundWithScale为true，则调整下一层图像的预留边界宽度  
26.         if( varyImgBoundWithScale ) curBound = static_cast<int>( curBound * featureScaleMul + 0.5f );  
27.     }  
28.     //掩码矩阵的特征点处理  
29.     KeyPointsFilter::runByPixelsMask( keypoints, mask );  
30. }  

void DenseFeatureDetector::detectImpl( const Mat& image, vector<KeyPoint>& keypoints, const Mat& mask ) const
{
    // curScale表示当前层图像特征点的尺度
    float curScale = static_cast<float>(initFeatureScale);
    //curStep表示当前层图像的网格宽度
    int curStep = initXyStep;
    //curBound表示当前层图像的预留边界宽度
    int curBound = initImgBound;
    //遍历各层图像
    for( int curLevel = 0; curLevel < featureScaleLevels; curLevel++ )
    {
        //遍历当前层图像的所有网格，图像四周的预留边界处是没有网格的，横、纵坐标的步长就是网格的宽度
        for( int x = curBound; x < image.cols - curBound; x += curStep )
        {
            for( int y = curBound; y < image.rows - curBound; y += curStep )
            {
                //把网格的左上角坐标处的像素作为该网格的特征点，并保存
                keypoints.push_back( KeyPoint(static_cast<float>(x), static_cast<float>(y), curScale) );
            }
        }
        //调整下一层图像特征点的尺度
        curScale = static_cast<float>(curScale * featureScaleMul);
        //如果varyXyStepWithScale为true，则调整下一层图像的网格宽度
        if( varyXyStepWithScale ) curStep = static_cast<int>( curStep * featureScaleMul + 0.5f );
        //如果varyImgBoundWithScale为true，则调整下一层图像的预留边界宽度
        if( varyImgBoundWithScale ) curBound = static_cast<int>( curBound * featureScaleMul + 0.5f );
    }
    //掩码矩阵的特征点处理
    KeyPointsFilter::runByPixelsMask( keypoints, mask );
}

下面给出一个具体的应用实例：

[cpp] view plain copy print ?

1. #include "opencv2/core/core.hpp"  
2. #include "highgui.h"  
3. #include "opencv2/imgproc/imgproc.hpp"  
4. #include "opencv2/features2d/features2d.hpp"  
5. #include "opencv2/nonfree/nonfree.hpp"  
6.   
7. using namespace cv;  
8. //using namespace std;  
9.   
10. int main(int argc, char** argv)  
11. {  
12.    Mat img = imread("box_in_scene.png"), img1;;  
13.      
14.    cvtColor( img, img1, CV_BGR2GRAY );  
15.   
16.    DenseFeatureDetector dense;  
17.   
18.    vector<KeyPoint> key_points;  
19.    Mat output_img;  
20.   
21.    dense.detect(img1,key_points,Mat());  
22.    drawKeypoints(img, key_points, output_img, Scalar::all(-1), DrawMatchesFlags::DRAW_RICH_KEYPOINTS);  
23.   
24.    namedWindow("DENSE");  
25.    imshow("DENSE", output_img);  
26.    waitKey(0);  
27.   
28.    return 0;  
29. }  

#include "opencv2/core/core.hpp"
#include "highgui.h"
#include "opencv2/imgproc/imgproc.hpp"
#include "opencv2/features2d/features2d.hpp"
#include "opencv2/nonfree/nonfree.hpp"

using namespace cv;
//using namespace std;

int main(int argc, char** argv)
{
   Mat img = imread("box_in_scene.png"), img1;;
   
   cvtColor( img, img1, CV_BGR2GRAY );

   DenseFeatureDetector dense;

   vector<KeyPoint> key_points;
   Mat output_img;

   dense.detect(img1,key_points,Mat());
   drawKeypoints(img, key_points, output_img, Scalar::all(-1), DrawMatchesFlags::DRAW_RICH_KEYPOINTS);

   namedWindow("DENSE");
   imshow("DENSE", output_img);
   waitKey(0);

   return 0;
}

phaseCorrelate
本文为博主原创文章，未经博主允许不得转载。

相位相关法（phase correlate）可以用于检测两幅内容相同的图像之间的相对位移量。它是基于傅立叶变换的位移定理：一个平移过的函数的傅立叶变换仅仅是未平移函数的傅立叶变换与一个具有线性相位的指数因子的乘积，即空间域中的平移会造成频域中频谱的相移。它的公式定义为：设二维函数（图像）f(x,y)的傅立叶变换为F(u,v)，即DFT[f(x,y)]=F(u,v)，如果f(x,y)平移(a,b)，则平移后的傅立叶变换为：

        （1）

因此，当两幅函数f₁(x,y)和f₂(x,y)仅仅有位移的差异，即f₂(x,y)= f₁(x-a,y-b)，则它们的傅立叶变换F₁(u,v)和F₂(u,v)有如下关系：

         （2）

由上式很容易得到f₁(x,y)和f₂(x,y)的互功率谱为（这里还用到了f₁(x,y)和f₂(x,y)的频谱的模相等）：

         （3）

式中F^*表示F的共轭，上式表示平移定理保证了互功率谱的相位等于两幅图像之间的相移。

Opencv的文档给出了详细的用相位相关法求解位移量的过程，

1、对待处理的两幅图像src1和src2应用窗函数去除图像的边界效应，文档中推荐使用汉宁窗，它可用createHanningWindow函数生成；

2、求傅立叶变换：G_a=DFT[scr1]和G_a=DFT[scr1]；

3、计算互功率谱： ；

4、对互功率谱求傅立叶逆变换：r=DFT^-1[R]；

5、对r计算最大值的位置，并在以该位置为中心的5×5的窗体内应用下列公式获得亚像素级的精度位置：

                （4）

最终(a,b)为两个图像之间的位移量。

 

phaseCorrelate函数的原型为：

Point2d phaseCorrelate(InputArray src1,InputArray src2, InputArray window=noArray())

src1和src2为两个要比较的图像，window为步骤1中所要用到的窗函数，该参数可选，默认情况下不使用任何窗函数，函数返回是具有亚像素级精度的位移量。

下面就给出具体的源码分析，该函数在source/modules/imgproc/src/phasecorr.cpp文件内，

[cpp] view plain copy print ?

1. cv::Point2d cv::phaseCorrelate(InputArray _src1, InputArray _src2, InputArray _window)  
2. {  
3.     return phaseCorrelateRes(_src1, _src2, _window, 0);  
4. }  

cv::Point2d cv::phaseCorrelate(InputArray _src1, InputArray _src2, InputArray _window)
{
    return phaseCorrelateRes(_src1, _src2, _window, 0);
}

phaseCorrelateRes函数的第4个参数表示的是最大响应值，也就是公式4中分母部分的归一化后的值，在这里我们没有用到它。

[cpp] view plain copy print ?

1. cv::Point2d cv::phaseCorrelateRes(InputArray _src1, InputArray _src2, InputArray _window, double* response)  
2. {  
3.     //分别得到两幅输入图像和窗函数的矩阵形式  
4.     Mat src1 = _src1.getMat();  
5.     Mat src2 = _src2.getMat();  
6.     Mat window = _window.getMat();  
7.     //输入图像的类型和大小的判断，必须一致，而且类型必须是32位或64位浮点灰度图像  
8.     CV_Assert( src1.type() == src2.type());  
9.     CV_Assert( src1.type() == CV_32FC1 || src1.type() == CV_64FC1 );  
10.     CV_Assert( src1.size == src2.size);  
11.     //如果使用窗函数，则窗函数的大小和类型必须与输入图像的一致  
12.     if(!window.empty())  
13.     {  
14.         CV_Assert( src1.type() == window.type());  
15.         CV_Assert( src1.size == window.size);  
16.     }  
17.     //因为要进行离散傅立叶变换，所以为了提高效率，就要得到最佳的图像尺寸  
18.     int M = getOptimalDFTSize(src1.rows);  
19.     int N = getOptimalDFTSize(src1.cols);  
20.   
21.     Mat padded1, padded2, paddedWin;  
22.     //生成尺寸修改以后的矩阵  
23.     if(M != src1.rows || N != src1.cols)   //最佳尺寸不是原图像的尺寸  
24.     {  
25.         //通过补零的方式填充多出来的像素  
26.         copyMakeBorder(src1, padded1, 0, M - src1.rows, 0, N - src1.cols, BORDER_CONSTANT, Scalar::all(0));  
27.         copyMakeBorder(src2, padded2, 0, M - src2.rows, 0, N - src2.cols, BORDER_CONSTANT, Scalar::all(0));  
28.   
29.         if(!window.empty())  
30.         {  
31.             copyMakeBorder(window, paddedWin, 0, M - window.rows, 0, N - window.cols, BORDER_CONSTANT, Scalar::all(0));  
32.         }  
33.     }  
34.     else    //最佳尺寸与原图像的尺寸一致  
35.     {  
36.         padded1 = src1;  
37.         padded2 = src2;  
38.         paddedWin = window;  
39.     }  
40.   
41.     Mat FFT1, FFT2, P, Pm, C;  
42.   
43.     // perform window multiplication if available  
44.     //执行步骤1，两幅输入图像分别与窗函数逐点相乘  
45.     if(!paddedWin.empty())  
46.     {  
47.         // apply window to both images before proceeding...  
48.         multiply(paddedWin, padded1, padded1);  
49.         multiply(paddedWin, padded2, padded2);  
50.     }  
51.   
52.     // execute phase correlation equation  
53.     // Reference: http://en.wikipedia.org/wiki/Phase_correlation  
54.     //执行步骤2，分别对两幅图像取傅立叶变换  
55.     dft(padded1, FFT1, DFT_REAL_OUTPUT);  
56.     dft(padded2, FFT2, DFT_REAL_OUTPUT);  
57.     //执行步骤3  
58.     //计算互功率谱的分子部分，即公式3中的分子，其中P为输出结果，true表示的是对FF2取共轭，所以得到的结果为：P=FFT1×FFT2*，mulSpectrums函数为通用函数  
59.     mulSpectrums(FFT1, FFT2, P, 0, true);  
60.     //计算互功率谱的分母部分，即公式3中的分母，结果为：Pm=|P|，magSpectrums函数就是在phasecorr.cpp文件内给出的，它的作用是对复数取模。  
61.     magSpectrums(P, Pm);  
62.     //计算互功率谱，即公式3，结果为：C=P / Pm，divSpectrums函数也是在phasecorr.cpp文件内给出的，它仿照mulSpectrums函数的写法，其中参数false表示不取共轭  
63.     divSpectrums(P, Pm, C, 0, false); // FF* / |FF*| (phase correlation equation completed here...)  
64.     //执行步骤4，傅立叶逆变换  
65.     idft(C, C); // gives us the nice peak shift location...  
66.     /*平移处理，fftShift函数也是在phasecorr.cpp文件内给出的，它的作用是把图像平均分割成——左上、左下、右上、右下，把左上和右下对调，把右上和左下对调。它的目的是把能量调整到图像的中心，也就是图像的中心对应于两幅图像相频差为零的地方，即没有发生位移的地方。*/  
67.     fftShift(C); // shift the energy to the center of the frame.  
68.   
69.     //执行步骤5  
70.     // locate the highest peak  
71.     //找到最大点处的像素位置，minMaxLoc为通用函数  
72.     Point peakLoc;  
73.     minMaxLoc(C, NULL, NULL, NULL, &peakLoc);  
74.   
75.     // get the phase shift with sub-pixel accuracy, 5x5 window seems about right here...  
76.     //在5×5的窗体内确定亚像素精度的坐标位置  
77.     Point2d t;  
78.     // weightedCentroid也是在phasecorr.cpp文件内给出的，它是利用公式4来计算更精确的坐标位置  
79.     t = weightedCentroid(C, peakLoc, Size(5, 5), response);  
80.   
81.     // max response is M*N (not exactly, might be slightly larger due to rounding errors)  
82.     //求最大响应值  
83.     if(response)  
84.         *response /= M*N;  
85.   
86.     // adjust shift relative to image center...  
87.     //最终确定位移量  
88.     //先找到图像中点，然后用中点减去由步骤5得到的坐标位置  
89.     Point2d center((double)padded1.cols / 2.0, (double)padded1.rows / 2.0);  
90.   
91.     return (center - t);  
92. }  

cv::Point2d cv::phaseCorrelateRes(InputArray _src1, InputArray _src2, InputArray _window, double* response)
{
    //分别得到两幅输入图像和窗函数的矩阵形式
    Mat src1 = _src1.getMat();
    Mat src2 = _src2.getMat();
    Mat window = _window.getMat();
    //输入图像的类型和大小的判断，必须一致，而且类型必须是32位或64位浮点灰度图像
    CV_Assert( src1.type() == src2.type());
    CV_Assert( src1.type() == CV_32FC1 || src1.type() == CV_64FC1 );
    CV_Assert( src1.size == src2.size);
    //如果使用窗函数，则窗函数的大小和类型必须与输入图像的一致
    if(!window.empty())
    {
        CV_Assert( src1.type() == window.type());
        CV_Assert( src1.size == window.size);
    }
    //因为要进行离散傅立叶变换，所以为了提高效率，就要得到最佳的图像尺寸
    int M = getOptimalDFTSize(src1.rows);
    int N = getOptimalDFTSize(src1.cols);

    Mat padded1, padded2, paddedWin;
    //生成尺寸修改以后的矩阵
    if(M != src1.rows || N != src1.cols)   //最佳尺寸不是原图像的尺寸
    {
        //通过补零的方式填充多出来的像素
        copyMakeBorder(src1, padded1, 0, M - src1.rows, 0, N - src1.cols, BORDER_CONSTANT, Scalar::all(0));
        copyMakeBorder(src2, padded2, 0, M - src2.rows, 0, N - src2.cols, BORDER_CONSTANT, Scalar::all(0));

        if(!window.empty())
        {
            copyMakeBorder(window, paddedWin, 0, M - window.rows, 0, N - window.cols, BORDER_CONSTANT, Scalar::all(0));
        }
    }
    else    //最佳尺寸与原图像的尺寸一致
    {
        padded1 = src1;
        padded2 = src2;
        paddedWin = window;
    }

    Mat FFT1, FFT2, P, Pm, C;

    // perform window multiplication if available
    //执行步骤1，两幅输入图像分别与窗函数逐点相乘
    if(!paddedWin.empty())
    {
        // apply window to both images before proceeding...
        multiply(paddedWin, padded1, padded1);
        multiply(paddedWin, padded2, padded2);
    }

    // execute phase correlation equation
    // Reference: http://en.wikipedia.org/wiki/Phase_correlation
    //执行步骤2，分别对两幅图像取傅立叶变换
    dft(padded1, FFT1, DFT_REAL_OUTPUT);
    dft(padded2, FFT2, DFT_REAL_OUTPUT);
    //执行步骤3
    //计算互功率谱的分子部分，即公式3中的分子，其中P为输出结果，true表示的是对FF2取共轭，所以得到的结果为：P=FFT1×FFT2*，mulSpectrums函数为通用函数
    mulSpectrums(FFT1, FFT2, P, 0, true);
    //计算互功率谱的分母部分，即公式3中的分母，结果为：Pm=|P|，magSpectrums函数就是在phasecorr.cpp文件内给出的，它的作用是对复数取模。
    magSpectrums(P, Pm);
    //计算互功率谱，即公式3，结果为：C=P / Pm，divSpectrums函数也是在phasecorr.cpp文件内给出的，它仿照mulSpectrums函数的写法，其中参数false表示不取共轭
    divSpectrums(P, Pm, C, 0, false); // FF* / |FF*| (phase correlation equation completed here...)
    //执行步骤4，傅立叶逆变换
    idft(C, C); // gives us the nice peak shift location...
    /*平移处理，fftShift函数也是在phasecorr.cpp文件内给出的，它的作用是把图像平均分割成——左上、左下、右上、右下，把左上和右下对调，把右上和左下对调。它的目的是把能量调整到图像的中心，也就是图像的中心对应于两幅图像相频差为零的地方，即没有发生位移的地方。*/
    fftShift(C); // shift the energy to the center of the frame.

    //执行步骤5
    // locate the highest peak
    //找到最大点处的像素位置，minMaxLoc为通用函数
    Point peakLoc;
    minMaxLoc(C, NULL, NULL, NULL, &peakLoc);

    // get the phase shift with sub-pixel accuracy, 5x5 window seems about right here...
    //在5×5的窗体内确定亚像素精度的坐标位置
    Point2d t;
    // weightedCentroid也是在phasecorr.cpp文件内给出的，它是利用公式4来计算更精确的坐标位置
    t = weightedCentroid(C, peakLoc, Size(5, 5), response);

    // max response is M*N (not exactly, might be slightly larger due to rounding errors)
    //求最大响应值
    if(response)
        *response /= M*N;

    // adjust shift relative to image center...
    //最终确定位移量
    //先找到图像中点，然后用中点减去由步骤5得到的坐标位置
    Point2d center((double)padded1.cols / 2.0, (double)padded1.rows / 2.0);

    return (center - t);
}

 

phaseCorrelate函数的源码可读性较强，而且它的原理在文档中说明得也十分清晰，因此对相位相关的理解难度应该不大。下面就给出具体的实现程序。我是对彩色图像进行测试，因为phaseCorrelate函数只能处理32位或64位浮点型的单通道图像，因此如果不是处理这类图像，是需要进行转换的。

[cpp] view plain copy print ?

1. #include "opencv2/core/core.hpp"  
2. #include "opencv2/highgui/highgui.hpp"  
3. #include "opencv2/imgproc/imgproc.hpp"  
4. #include <iostream>  
5. using namespace cv;  
6. using namespace std;  
7. int main( int argc, char** argv )  
8. {  
9.     Mat src1, src2;  
10.     Mat dst1, dst2;  
11.     if( argc != 3 || !(src1=imread(argv[1], 1)).data || !(src2=imread(argv[2], 1)).data)  
12.         return -1;  
13.   
14.     cvtColor( src1, src1, CV_BGR2GRAY );     //转换为灰度图像  
15.         src1.convertTo(dst1,CV_32FC1);       //转换为32位浮点型  
16.     cvtColor( src2, src2, CV_BGR2GRAY );  
17.         src2.convertTo(dst2,CV_32FC1);  
18.   
19.     Point2d phase_shift;  
20.         phase_shift = phaseCorrelate(dst1,dst2);  
21.         cout<<endl<<"warp :"<<endl<<"\tX shift : "<<phase_shift.x<<"\tY shift : "<<phase_shift.y<<endl;  
22.   
23.         waitKey(0);  
24.   
25.         return 0;  
26. }  

#include "opencv2/core/core.hpp"
#include "opencv2/highgui/highgui.hpp"
#include "opencv2/imgproc/imgproc.hpp"
#include <iostream>
using namespace cv;
using namespace std;
int main( int argc, char** argv )
{
	Mat src1, src2;
	Mat dst1, dst2;
	if( argc != 3 || !(src1=imread(argv[1], 1)).data || !(src2=imread(argv[2], 1)).data)
		return -1;

	cvtColor( src1, src1, CV_BGR2GRAY );     //转换为灰度图像
        src1.convertTo(dst1,CV_32FC1);       //转换为32位浮点型
	cvtColor( src2, src2, CV_BGR2GRAY );
        src2.convertTo(dst2,CV_32FC1);

	Point2d phase_shift;
        phase_shift = phaseCorrelate(dst1,dst2);
        cout<<endl<<"warp :"<<endl<<"\tX shift : "<<phase_shift.x<<"\tY shift : "<<phase_shift.y<<endl;

        waitKey(0);

        return 0;
}

 

最终得到的结果是：

warp:

      X  shift :  3.803         Y  shift :  21.8661

我们来验证一下上述结果是否正确。我用的测试图像（如下所示）的原始大小为1761×1761，两幅图像的相对位移量为（16，96），然后我把两幅图像的大小调整为400×400，这样经过计算，相对位移量应该是（3.6343，21.8058）。可以看出两者之间还是存在一定的误差的，但我认为误差还是能够在可接受的范围内的。

Opencv2.4.9源码分析——HoughLines

霍夫变换是一种特征提取技术。经典的霍夫变换能够识别出图像中的直线，后来又发展到能够识别出任意形状，但更常见的识别形状是圆和椭圆。函数HoughLines的作用就是能够检测出边缘图像中的直线。

在直角坐标系下，直线被定义为：

y = mx + b         （1）

其中，m为斜率，b为与y轴的截距，只要确定了m和b，一条直线就可以被唯一地确定下来。如果用ρ₀表示原点到该直线的代数距离，θ₀表示该直线的正交线与x轴的夹角，则：

  （2）

      （3）

则该直线又可表示为：

     （4）

写成更一般的形式：

       （5）

很容易想到，（ρ，θ）是极坐标的表示形式。但如果把（ρ，θ）也用直角坐标的形式表示，即把ρ和θ做正交处理，则（ρ，θ）就被称为霍夫空间。

在直角坐标系中的一点，对应于霍夫空间的一条正弦曲线。直线是由无数个点组成的，在霍夫空间就是无数条正弦曲线，但这些正弦曲线会相交于一点（ρ₀，θ₀），把该点带入公式2和公式3就得了直线的斜率和截距，这样一条直线就被确定了下来。因此在用霍夫变换识别直线时，霍夫空间中的极大值就有可能对应一条直线。

下面给出函数HoughLines的霍夫变换直线检测的步骤：

1、对边缘图像进行霍夫空间变换；

2、在4邻域内找到霍夫空间变换的极大值；

3、对这些极大值安装由大到小顺序进行排序，极大值越大，越有可能是直线；

4、输出直线。

 

函数HoughLines的原型为：

void HoughLines(InputArray image,OutputArray lines, double rho, double theta, int threshold, double srn=0,double stn=0 )

image为输入图像，要求是单通道的二值图像

lines为输出直线向量，两个元素的向量（ρ，θ）代表一条直线，ρ是从原点（图像的左上角）的距离，θ是直线的角度（单位是弧度），0表示垂直线，π/2表示水平线

rho为距离分辨率

theta为角度分辨率

threshold为阈值，在步骤2中，只有大于该值的点才有可能被当作极大值，即至少有多少条正弦曲线交于一点才被认为是直线

srn和stn在多尺度霍夫变换的时候才会使用，在这里我们只研究经典霍夫变换

 

 

HoughLines函数是在sources/modules/imgproc/src/hough.cpp文件中定义的：

[cpp] view plain copy print ?

1. void cv::HoughLines( InputArray _image, OutputArray _lines,  
2.                      double rho, double theta, int threshold,  
3.                      double srn, double stn )  
4. {  
5.     //申请一段内存，用于存储霍夫变换后所检测到的直线  
6.     Ptr<CvMemStorage> storage = cvCreateMemStorage(STORAGE_SIZE);  
7.     //提取输入图像矩阵  
8.     Mat image = _image.getMat();  
9.     CvMat c_image = image;  
10.     //调用由C语音写出的霍夫变换的函数  
11.     CvSeq* seq = cvHoughLines2( &c_image, storage, srn == 0 && stn == 0 ?  
12.                     CV_HOUGH_STANDARD : CV_HOUGH_MULTI_SCALE,  
13.                     rho, theta, threshold, srn, stn );  
14.     //把由cvHoughLines2函数得到的霍夫变换直线序列转换为数组  
15.     seqToMat(seq, _lines);  
16. }  

void cv::HoughLines( InputArray _image, OutputArray _lines,
                     double rho, double theta, int threshold,
                     double srn, double stn )
{
    //申请一段内存，用于存储霍夫变换后所检测到的直线
    Ptr<CvMemStorage> storage = cvCreateMemStorage(STORAGE_SIZE);
    //提取输入图像矩阵
    Mat image = _image.getMat();
    CvMat c_image = image;
    //调用由C语音写出的霍夫变换的函数
    CvSeq* seq = cvHoughLines2( &c_image, storage, srn == 0 && stn == 0 ?
                    CV_HOUGH_STANDARD : CV_HOUGH_MULTI_SCALE,
                    rho, theta, threshold, srn, stn );
    //把由cvHoughLines2函数得到的霍夫变换直线序列转换为数组
    seqToMat(seq, _lines);
}

 

cvHoughLines2函数是霍夫变换直线检测的关键，函数的输出为所检测到的直线序列，它的第3个形参“srn == 0 && stn == 0 ? CV_HOUGH_STANDARD :CV_HOUGH_MULTI_SCALE”表示的是该霍夫变换是经典霍夫变换还是多尺度霍夫变换，它是由变量srn和stn决定的，只要这两个变量有一个不为0，就进行多尺度霍夫变换，否则为经典霍夫变换。另外cvHoughLines2函数不仅可以用于经典霍夫变换和多尺度霍夫变换，还可以用于概率霍夫变换。

[cpp] view plain copy print ?

1. CV_IMPL CvSeq*  
2. cvHoughLines2( CvArr* src_image, void* lineStorage, int method,  
3.                double rho, double theta, int threshold,  
4.                double param1, double param2 )  
5. {  
6.     CvSeq* result = 0;  
7.   
8.     CvMat stub, *img = (CvMat*)src_image;  
9.     CvMat* mat = 0;  
10.     CvSeq* lines = 0;  
11.     CvSeq lines_header;  
12.     CvSeqBlock lines_block;  
13.     int lineType, elemSize;  
14.     int linesMax = INT_MAX;    //输出最多直线的数量，设为无穷多  
15.     int iparam1, iparam2;  
16.   
17.     img = cvGetMat( img, &stub );  
18.     //确保输入图像是8位单通道  
19.     if( !CV_IS_MASK_ARR(img))  
20.         CV_Error( CV_StsBadArg, "The source image must be 8-bit, single-channel" );  
21.     //内存空间申请成功，以备输出之用  
22.     if( !lineStorage )  
23.         CV_Error( CV_StsNullPtr, "NULL destination" );  
24.     //确保rho，theta，threshold这三个参数大于0  
25.     if( rho <= 0 || theta <= 0 || threshold <= 0 )  
26.         CV_Error( CV_StsOutOfRange, "rho, theta and threshold must be positive" );  
27.   
28.     if( method != CV_HOUGH_PROBABILISTIC )    //经典霍夫变换和多尺度霍夫变换  
29.     {  
30.         //输出直线的类型为32位浮点双通道，即ρ和θ两个浮点型变量  
31.         lineType = CV_32FC2;  
32.         elemSize = sizeof(float)*2;  
33.     }  
34.     else    //概率霍夫变换  
35.     {  
36.         //输出直线的类型为32位有符号整型4通道，即两个像素点的4个坐标  
37.         lineType = CV_32SC4;  
38.         elemSize = sizeof(int)*4;  
39.     }  
40.     //判断lineStorage的类型，经分析lineStorage只可能是STORAGE类型  
41.     if( CV_IS_STORAGE( lineStorage ))  
42.     {  
43.         //在lineStorage内存中创建一个序列，用于存储霍夫变换的直线，lines为该序列的指针  
44.         lines = cvCreateSeq( lineType, sizeof(CvSeq), elemSize, (CvMemStorage*)lineStorage );  
45.     }  
46.     else if( CV_IS_MAT( lineStorage ))  
47.     {  
48.         mat = (CvMat*)lineStorage;  
49.   
50.         if( !CV_IS_MAT_CONT( mat->type ) || (mat->rows != 1 && mat->cols != 1) )  
51.             CV_Error( CV_StsBadArg,  
52.             "The destination matrix should be continuous and have a single row or a single column" );  
53.   
54.         if( CV_MAT_TYPE( mat->type ) != lineType )  
55.             CV_Error( CV_StsBadArg,  
56.             "The destination matrix data type is inappropriate, see the manual" );  
57.   
58.         lines = cvMakeSeqHeaderForArray( lineType, sizeof(CvSeq), elemSize, mat->data.ptr,  
59.                                          mat->rows + mat->cols - 1, &lines_header, &lines_block );  
60.         linesMax = lines->total;  
61.         cvClearSeq( lines );  
62.     }  
63.     else  
64.         CV_Error( CV_StsBadArg, "Destination is not CvMemStorage* nor CvMat*" );  
65.   
66.     iparam1 = cvRound(param1);  
67.     iparam2 = cvRound(param2);  
68.   
69.     switch( method )  
70.     {  
71.     case CV_HOUGH_STANDARD:    //经典霍夫变换  
72.           icvHoughLinesStandard( img, (float)rho,  
73.                 (float)theta, threshold, lines, linesMax );  
74.           break;  
75.     case CV_HOUGH_MULTI_SCALE:    //多尺度霍夫变换  
76.           icvHoughLinesSDiv( img, (float)rho, (float)theta,  
77.                 threshold, iparam1, iparam2, lines, linesMax );  
78.           break;  
79.     case CV_HOUGH_PROBABILISTIC:    //概率霍夫变换  
80.           icvHoughLinesProbabilistic( img, (float)rho, (float)theta,  
81.                 threshold, iparam1, iparam2, lines, linesMax );  
82.           break;  
83.     default:  
84.         CV_Error( CV_StsBadArg, "Unrecognized method id" );  
85.     }  
86.     //在前面判断lineStorage类型时，已经确定它是STORAGE类型，因此没有进入else if内，也就是没有对mat赋值，所以在这里进入的是else  
87.     if( mat )  
88.     {  
89.         if( mat->cols > mat->rows )  
90.             mat->cols = lines->total;  
91.         else  
92.             mat->rows = lines->total;  
93.     }  
94.     else  
95.         result = lines;  
96.     //返回lines序列的指针，该序列存储有霍夫变换的直线  
97.     return result;  
98. }  

CV_IMPL CvSeq*
cvHoughLines2( CvArr* src_image, void* lineStorage, int method,
               double rho, double theta, int threshold,
               double param1, double param2 )
{
    CvSeq* result = 0;

    CvMat stub, *img = (CvMat*)src_image;
    CvMat* mat = 0;
    CvSeq* lines = 0;
    CvSeq lines_header;
    CvSeqBlock lines_block;
    int lineType, elemSize;
    int linesMax = INT_MAX;    //输出最多直线的数量，设为无穷多
    int iparam1, iparam2;

    img = cvGetMat( img, &stub );
    //确保输入图像是8位单通道
    if( !CV_IS_MASK_ARR(img))
        CV_Error( CV_StsBadArg, "The source image must be 8-bit, single-channel" );
    //内存空间申请成功，以备输出之用
    if( !lineStorage )
        CV_Error( CV_StsNullPtr, "NULL destination" );
    //确保rho，theta，threshold这三个参数大于0
    if( rho <= 0 || theta <= 0 || threshold <= 0 )
        CV_Error( CV_StsOutOfRange, "rho, theta and threshold must be positive" );

    if( method != CV_HOUGH_PROBABILISTIC )    //经典霍夫变换和多尺度霍夫变换
    {
        //输出直线的类型为32位浮点双通道，即ρ和θ两个浮点型变量
        lineType = CV_32FC2;
        elemSize = sizeof(float)*2;
    }
    else    //概率霍夫变换
    {
        //输出直线的类型为32位有符号整型4通道，即两个像素点的4个坐标
        lineType = CV_32SC4;
        elemSize = sizeof(int)*4;
    }
    //判断lineStorage的类型，经分析lineStorage只可能是STORAGE类型
    if( CV_IS_STORAGE( lineStorage ))
    {
        //在lineStorage内存中创建一个序列，用于存储霍夫变换的直线，lines为该序列的指针
        lines = cvCreateSeq( lineType, sizeof(CvSeq), elemSize, (CvMemStorage*)lineStorage );
    }
    else if( CV_IS_MAT( lineStorage ))
    {
        mat = (CvMat*)lineStorage;

        if( !CV_IS_MAT_CONT( mat->type ) || (mat->rows != 1 && mat->cols != 1) )
            CV_Error( CV_StsBadArg,
            "The destination matrix should be continuous and have a single row or a single column" );

        if( CV_MAT_TYPE( mat->type ) != lineType )
            CV_Error( CV_StsBadArg,
            "The destination matrix data type is inappropriate, see the manual" );

        lines = cvMakeSeqHeaderForArray( lineType, sizeof(CvSeq), elemSize, mat->data.ptr,
                                         mat->rows + mat->cols - 1, &lines_header, &lines_block );
        linesMax = lines->total;
        cvClearSeq( lines );
    }
    else
        CV_Error( CV_StsBadArg, "Destination is not CvMemStorage* nor CvMat*" );

    iparam1 = cvRound(param1);
    iparam2 = cvRound(param2);

    switch( method )
    {
    case CV_HOUGH_STANDARD:    //经典霍夫变换
          icvHoughLinesStandard( img, (float)rho,
                (float)theta, threshold, lines, linesMax );
          break;
    case CV_HOUGH_MULTI_SCALE:    //多尺度霍夫变换
          icvHoughLinesSDiv( img, (float)rho, (float)theta,
                threshold, iparam1, iparam2, lines, linesMax );
          break;
    case CV_HOUGH_PROBABILISTIC:    //概率霍夫变换
          icvHoughLinesProbabilistic( img, (float)rho, (float)theta,
                threshold, iparam1, iparam2, lines, linesMax );
          break;
    default:
        CV_Error( CV_StsBadArg, "Unrecognized method id" );
    }
    //在前面判断lineStorage类型时，已经确定它是STORAGE类型，因此没有进入else if内，也就是没有对mat赋值，所以在这里进入的是else
    if( mat )
    {
        if( mat->cols > mat->rows )
            mat->cols = lines->total;
        else
            mat->rows = lines->total;
    }
    else
        result = lines;
    //返回lines序列的指针，该序列存储有霍夫变换的直线
    return result;
}

 

对于经典霍夫变换，调用的是icvHoughLinesStandard函数：

[cpp] view plain copy print ?

1. static void  
2. icvHoughLinesStandard( const CvMat* img, float rho, float theta,  
3.                        int threshold, CvSeq *lines, int linesMax )  
4. {  
5.     cv::AutoBuffer<int> _accum, _sort_buf;  
6.     cv::AutoBuffer<float> _tabSin, _tabCos;  
7.   
8.     const uchar* image;  
9.     int step, width, height;  
10.     int numangle, numrho;  
11.     int total = 0;  
12.     int i, j;  
13.     float irho = 1 / rho;  
14.     double scale;  
15.     //再次确保输入图像的正确性  
16.     CV_Assert( CV_IS_MAT(img) && CV_MAT_TYPE(img->type) == CV_8UC1 );  
17.   
18.     image = img->data.ptr;    //得到图像的指针  
19.     step = img->step;    //得到图像的步长  
20.     width = img->cols;    //得到图像的宽  
21.     height = img->rows;    //得到图像的高  
22.     //由角度和距离的分辨率得到角度和距离的数量，即霍夫变换后角度和距离的个数  
23.     numangle = cvRound(CV_PI / theta);  
24.     numrho = cvRound(((width + height) * 2 + 1) / rho);  
25.   
26.     //为累加器数组分配内存空间  
27.     //该累加器数组其实就是霍夫空间，它是用一维数组表示二维空间  
28.     _accum.allocate((numangle+2) * (numrho+2));  
29.     //为排序数组分配内存空间  
30.     _sort_buf.allocate(numangle * numrho);  
31.     //为正弦和余弦列表分配内存空间  
32.     _tabSin.allocate(numangle);  
33.     _tabCos.allocate(numangle);  
34.     //分别定义上述内存空间的地址指针  
35.     int *accum = _accum, *sort_buf = _sort_buf;  
36.     float *tabSin = _tabSin, *tabCos = _tabCos;  
37.     //累加器数组清零  
38.     memset( accum, 0, sizeof(accum[0]) * (numangle+2) * (numrho+2) );  
39.   
40.     float ang = 0;  
41.     //为避免重复运算，事先计算好sinθi/ρ和cosθi/ρ  
42.     for(int n = 0; n < numangle; ang += theta, n++ )  
43.     {  
44.         tabSin[n] = (float)(sin((double)ang) * irho);  
45.         tabCos[n] = (float)(cos((double)ang) * irho);  
46.     }  
47.   
48.     // stage 1. fill accumulator  
49.     //执行步骤1，逐点进行霍夫空间变换，并把结果放入累加器数组内  
50.     for( i = 0; i < height; i++ )  
51.         for( j = 0; j < width; j++ )  
52.         {  
53.             //只对图像的非零值处理，即只对图像的边缘像素进行霍夫变换  
54.             if( image[i * step + j] != 0 )  
55.                 for(int n = 0; n < numangle; n++ )  
56.                 {  
57.                     int r = cvRound( j * tabCos[n] + i * tabSin[n] );  
58.                     r += (numrho - 1) / 2;  
59.                     //r表示的是距离，n表示的是角点，在累加器内找到它们所对应的位置（即霍夫空间内的位置），其值加1  
60.                     accum[(n+1) * (numrho+2) + r+1]++;  
61.                 }  
62.         }  
63.   
64.     // stage 2. find local maximums  
65.     //执行步骤2，找到局部极大值，即非极大值抑制  
66.     for(int r = 0; r < numrho; r++ )  
67.         for(int n = 0; n < numangle; n++ )  
68.         {  
69.             //得到当前值在累加器数组的位置  
70.             int base = (n+1) * (numrho+2) + r+1;  
71.             if( accum[base] > threshold &&    //必须大于所设置的阈值  
72.                 //在4邻域内进行非极大值抑制  
73.                 accum[base] > accum[base - 1] && accum[base] >= accum[base + 1] &&  
74.                 accum[base] > accum[base - numrho - 2] && accum[base] >= accum[base + numrho + 2] )  
75.                 //把极大值位置存入排序数组内——sort_buf  
76.                 sort_buf[total++] = base;  
77.         }  
78.   
79.     // stage 3. sort the detected lines by accumulator value  
80.     //执行步骤3，对存储在sort_buf数组内的累加器的数据按由大到小的顺序进行排序  
81.     icvHoughSortDescent32s( sort_buf, total, accum );  
82.   
83.     // stage 4. store the first min(total,linesMax) lines to the output buffer  
84.     //执行步骤4，只输出前linesMax条直线，这里linesMax就等于total，即输出所有直线  
85.     linesMax = MIN(linesMax, total);  
86.     //事先定义一个尺度  
87.     scale = 1./(numrho+2);  
88.     for( i = 0; i < linesMax; i++ )  
89.     {  
90.         // CvLinePolar结构在该文件的前面被定义  
91.         CvLinePolar line;  
92.         //idx为极大值在累加器数组的位置  
93.         int idx = sort_buf[i];  
94.         //分离出该极大值在霍夫空间中的位置  
95.         int n = cvFloor(idx*scale) - 1;  
96.         int r = idx - (n+1)*(numrho+2) - 1;  
97.         //最终得到极大值所对应的角度和距离  
98.         line.rho = (r - (numrho - 1)*0.5f) * rho;  
99.         line.angle = n * theta;  
100.         //存储到序列内  
101.         cvSeqPush( lines, &line );  
102.     }  
103. }  

static void
icvHoughLinesStandard( const CvMat* img, float rho, float theta,
                       int threshold, CvSeq *lines, int linesMax )
{
    cv::AutoBuffer<int> _accum, _sort_buf;
    cv::AutoBuffer<float> _tabSin, _tabCos;

    const uchar* image;
    int step, width, height;
    int numangle, numrho;
    int total = 0;
    int i, j;
    float irho = 1 / rho;
    double scale;
    //再次确保输入图像的正确性
    CV_Assert( CV_IS_MAT(img) && CV_MAT_TYPE(img->type) == CV_8UC1 );

    image = img->data.ptr;    //得到图像的指针
    step = img->step;    //得到图像的步长
    width = img->cols;    //得到图像的宽
    height = img->rows;    //得到图像的高
    //由角度和距离的分辨率得到角度和距离的数量，即霍夫变换后角度和距离的个数
    numangle = cvRound(CV_PI / theta);
    numrho = cvRound(((width + height) * 2 + 1) / rho);

    //为累加器数组分配内存空间
    //该累加器数组其实就是霍夫空间，它是用一维数组表示二维空间
    _accum.allocate((numangle+2) * (numrho+2));
    //为排序数组分配内存空间
    _sort_buf.allocate(numangle * numrho);
    //为正弦和余弦列表分配内存空间
    _tabSin.allocate(numangle);
    _tabCos.allocate(numangle);
    //分别定义上述内存空间的地址指针
    int *accum = _accum, *sort_buf = _sort_buf;
    float *tabSin = _tabSin, *tabCos = _tabCos;
    //累加器数组清零
    memset( accum, 0, sizeof(accum[0]) * (numangle+2) * (numrho+2) );

    float ang = 0;
    //为避免重复运算，事先计算好sinθi/ρ和cosθi/ρ
    for(int n = 0; n < numangle; ang += theta, n++ )
    {
        tabSin[n] = (float)(sin((double)ang) * irho);
        tabCos[n] = (float)(cos((double)ang) * irho);
    }

    // stage 1. fill accumulator
    //执行步骤1，逐点进行霍夫空间变换，并把结果放入累加器数组内
    for( i = 0; i < height; i++ )
        for( j = 0; j < width; j++ )
        {
            //只对图像的非零值处理，即只对图像的边缘像素进行霍夫变换
            if( image[i * step + j] != 0 )
                for(int n = 0; n < numangle; n++ )
                {
                    int r = cvRound( j * tabCos[n] + i * tabSin[n] );
                    r += (numrho - 1) / 2;
                    //r表示的是距离，n表示的是角点，在累加器内找到它们所对应的位置（即霍夫空间内的位置），其值加1
                    accum[(n+1) * (numrho+2) + r+1]++;
                }
        }

    // stage 2. find local maximums
    //执行步骤2，找到局部极大值，即非极大值抑制
    for(int r = 0; r < numrho; r++ )
        for(int n = 0; n < numangle; n++ )
        {
            //得到当前值在累加器数组的位置
            int base = (n+1) * (numrho+2) + r+1;
            if( accum[base] > threshold &&    //必须大于所设置的阈值
                //在4邻域内进行非极大值抑制
                accum[base] > accum[base - 1] && accum[base] >= accum[base + 1] &&
                accum[base] > accum[base - numrho - 2] && accum[base] >= accum[base + numrho + 2] )
                //把极大值位置存入排序数组内——sort_buf
                sort_buf[total++] = base;
        }

    // stage 3. sort the detected lines by accumulator value
    //执行步骤3，对存储在sort_buf数组内的累加器的数据按由大到小的顺序进行排序
    icvHoughSortDescent32s( sort_buf, total, accum );

    // stage 4. store the first min(total,linesMax) lines to the output buffer
    //执行步骤4，只输出前linesMax条直线，这里linesMax就等于total，即输出所有直线
    linesMax = MIN(linesMax, total);
    //事先定义一个尺度
    scale = 1./(numrho+2);
    for( i = 0; i < linesMax; i++ )
    {
        // CvLinePolar结构在该文件的前面被定义
        CvLinePolar line;
        //idx为极大值在累加器数组的位置
        int idx = sort_buf[i];
        //分离出该极大值在霍夫空间中的位置
        int n = cvFloor(idx*scale) - 1;
        int r = idx - (n+1)*(numrho+2) - 1;
        //最终得到极大值所对应的角度和距离
        line.rho = (r - (numrho - 1)*0.5f) * rho;
        line.angle = n * theta;
        //存储到序列内
        cvSeqPush( lines, &line );
    }
}

由于函数HoughLines只输出直线所对应的角度和距离，所以在进行直线检测的时候还要把其转换为直角坐标系下的数据，另外输入图像还必须是边缘图像，下面就是具体的实例：

[cpp] view plain copy print ?

1. #include "opencv2/core/core.hpp"  
2. #include "opencv2/highgui/highgui.hpp"  
3. #include "opencv2/imgproc/imgproc.hpp"  
4. #include <iostream>  
5. using namespace cv;  
6. using namespace std;  
7.   
8. int main( int argc, char** argv )  
9. {  
10.     Mat src, edge;  
11.     src=imread("building.jpg");  
12.     if( !src.data )    
13.         return -1;    
14.     //Canny边缘检测  
15.     Canny(src,edge,50,200,3);  
16.     //定义输出数组，用于存储直线的角度和距离这两个变量  
17.     vector<Vec2f> lines;  
18.     //距离分辨率为1，角度分辨率为π/180，阈值为215  
19.     //阈值的选取直接影响到输出直线的数量  
20.     HoughLines(edge,lines,1,CV_PI/180,215,0,0);  
21.     //画直线  
22.     for( size_t i = 0; i < lines.size(); i++ )  
23.     {  
24.         //提取出距离和角度  
25.         float rho = lines[i][0], theta = lines[i][1];  
26.         //定义两个点，两点确定一条直线  
27.         //计算得到的两点的坐标为（ρcosθ-1000sinθ，ρsinθ+1000cosθ），（ρcosθ+1000sinθ，ρsinθ-1000cosθ）  
28.         Point pt1, pt2;  
29.         double a = cos(theta), b = sin(theta);  
30.         double x0 = a*rho, y0 = b*rho;  
31.         pt1.x = cvRound(x0 + 1000*(-b));  
32.         pt1.y = cvRound(y0 + 1000*(a));  
33.         pt2.x = cvRound(x0 - 1000*(-b));  
34.         pt2.y = cvRound(y0 - 1000*(a));  
35.         //在原图上画宽带为2的红线  
36.         line( src, pt1, pt2, Scalar(0,0,255),2);  
37.     }  
38.    
39.     namedWindow( "lines", CV_WINDOW_AUTOSIZE );  
40.     imshow( "lines", src );  
41.     waitKey(0);  
42.   
43.     return 0;  
44. }  

#include "opencv2/core/core.hpp"
#include "opencv2/highgui/highgui.hpp"
#include "opencv2/imgproc/imgproc.hpp"
#include <iostream>
using namespace cv;
using namespace std;

int main( int argc, char** argv )
{
	Mat src, edge;
	src=imread("building.jpg");
	if( !src.data )  
		return -1;  
    //Canny边缘检测
	Canny(src,edge,50,200,3);
    //定义输出数组，用于存储直线的角度和距离这两个变量
	vector<Vec2f> lines;
    //距离分辨率为1，角度分辨率为π/180，阈值为215
    //阈值的选取直接影响到输出直线的数量
	HoughLines(edge,lines,1,CV_PI/180,215,0,0);
    //画直线
    for( size_t i = 0; i < lines.size(); i++ )
    {
        //提取出距离和角度
        float rho = lines[i][0], theta = lines[i][1];
        //定义两个点，两点确定一条直线
        //计算得到的两点的坐标为（ρcosθ-1000sinθ，ρsinθ+1000cosθ），（ρcosθ+1000sinθ，ρsinθ-1000cosθ）
        Point pt1, pt2;
        double a = cos(theta), b = sin(theta);
        double x0 = a*rho, y0 = b*rho;
        pt1.x = cvRound(x0 + 1000*(-b));
        pt1.y = cvRound(y0 + 1000*(a));
        pt2.x = cvRound(x0 - 1000*(-b));
        pt2.y = cvRound(y0 - 1000*(a));
        //在原图上画宽带为2的红线
        line( src, pt1, pt2, Scalar(0,0,255),2);
    }
 
    namedWindow( "lines", CV_WINDOW_AUTOSIZE );
    imshow( "lines", src );
    waitKey(0);

    return 0;
}

下图霍夫变换直线检测的结果：