OpenCV图像处理：霍夫圆检测与最小二乘法拟合圆详解

OpenCV图像处理：霍夫圆检测与最小二乘法拟合圆详解

引用

CSDN

1.

https://m.blog.csdn.net/matt45m/article/details/140417972

在图像处理领域，圆检测是一个常见的任务，广泛应用于各种场景，如工业检测、医学影像分析等。本文将详细介绍两种常用的圆检测方法：霍夫圆检测和最小二乘法拟合圆，并提供具体的实现代码。

霍夫圆检测

霍夫圆检测是基于霍夫变换的一种圆检测算法，通过在参数空间中寻找圆的参数来检测图像中的圆。OpenCV提供了HoughCircles函数来实现这一功能。

函数原型

void HoughCircles( InputArray image, OutputArray circles, int method, double dp, double minDist, double param1 = 100, double param2 = 100, int minRadius = 0, int maxRadius = 0 );

参数说明

• image：输入的8-bit单通道灰度图
• circles：输出的圆的向量，每个圆包含3或4个元素（x坐标、y坐标、半径、投票数）
• method：检测方法，可选值包括：
• HOUGH_STANDARD：标准霍夫变换
• HOUGH_PROBABILISTIC：基于概率的霍夫变换
• HOUGH_MULTI_SCALE：标准霍夫变换的多尺度变种
• HOUGH_GRADIENT：霍夫梯度
• dp：累加器分辨率与图像分辨率的反比
• minDist：检测到的圆的中心之间的最小距离
• param1：如果检测方法为HOUGH_GRADIENT，代表Canny边缘提取的高阈值
• param2：如果检测方法为HOUGH_GRADIENT，代表检测圆心的累加器的阈值
• minRadius：最小半径
• maxRadius：最大半径

使用步骤

1. 转换为灰度图像：如果图像不是灰度的，需要先转换为灰度图像。
2. 使用高斯模糊：对图像应用高斯模糊，以减少噪声和细节。
3. 边缘检测：使用Canny算法或其他边缘检测方法来获取图像的边缘。
4. 霍夫圆检测：使用HoughCircles函数检测边缘图像中的圆形。
5. 绘制圆：在原始图像上绘制检测到的圆。

示例代码

#include <opencv2/opencv.hpp>
#include <iostream>

using namespace cv;
using namespace std;

int main() {
    // 读取图像
    Mat image = imread("path_to_your_image.jpg", IMREAD_COLOR);
    if (image.empty()) {
        cerr << "Could not read the image" << endl;
        return 1;
    }

    // 转换为灰度图像
    Mat gray;
    cvtColor(image, gray, COLOR_BGR2GRAY);

    // 使用高斯模糊
    Mat blurred;
    GaussianBlur(gray, blurred, Size(5, 5), 0);

    // Canny 边缘检测
    Mat edges;
    Canny(blurred, edges, 100, 200);

    // 霍夫圆检测
    vector<Vec3f> circles;
    HoughCircles(edges, circles,
                 3, // 霍夫变换的分辨率
                 blurred.rows / 8, // 霍夫变换的灵敏度，较小的值更敏感
                 100, // 阈值，确定何时一个圆被认为是检测到的
                 30, // 最小圆半径
                 200); // 最大圆半径

    // 绘制圆
    for (size_t i = 0; i < circles.size(); i++) {
        Vec3f circle = circles[i];
        Point center(cvRound(circle[0]), cvRound(circle[1]));
        int radius = cvRound(circle[2]);
        circle(image, center, radius, Scalar(0, 255, 0), 3, 8, 0);
    }

    // 显示结果
    imshow("Detected Circles", image);
    waitKey(0);
    destroyAllWindows();

    return 0;
}

最小二乘法拟合圆

最小二乘法是一种数学优化技术，通过最小化误差的平方和来寻找数据的最佳函数匹配。在计算机视觉和图像处理中，最小二乘法常用于几何形状的拟合，比如圆。

对于圆的拟合，我们通常有一组二维点 ((x_i, y_i))，并希望找到一个圆，使得这些点到圆周的垂直距离之和的平方最小。圆的方程可以表示为：

[
(x - a)^2 + (y - b)^2 = r^2
]

其中 ((a, b)) 是圆心坐标，(r) 是半径。

使用步骤

1. 建立目标函数：定义目标函数为所有点到圆周的垂直距离平方和。
2. 建立法线方程：对于圆上的每个点 ((x_i, y_i))，从圆心 ((a, b)) 到该点的线段的斜率 (m) 为 ((y_i - b) / (x_i - a))。法线与此线段垂直，所以法线的斜率为 (-1/m)。
3. 最小化误差：将法线方程代入圆的方程，并最小化目标函数。这通常涉及到对目标函数求导并令导数为零求解。
4. 求解方程组：将得到的方程组求解，得到圆心 ((a, b)) 和半径 (r)。
5. 验证和优化：检查解的合理性，并对结果进行必要的优化或迭代。

C++实现代码

Circle3f circle_least_squares(std::vector<cv::Point2f> pts)
{
    Circle3f lsc;

    int adj_count = 0;

LABEL_ADJ:

    const int pt_len = (int)pts.size();

    if (pt_len < 3)
    {
        lsc.center.x = 0;
        lsc.center.y = 0;
        lsc.radius = 0;

        return lsc;
    }

    cv::Mat A(pt_len, 3, CV_32FC1);
    cv::Mat b(pt_len, 1, CV_32FC1);

    for (int r = 0; r < pt_len; r++)
    {
        float* pData = A.ptr<float>(r);

        pData[0] = pts[r].x * 2.0f;
        pData[1] = pts[r].y * 2.0f;
        pData[2] = 1.0f;
    }

    float* pb = (float*)b.data;

    for (int i = 0; i < pt_len; i++)
    {
        pb[i] = (float)(pts[i].x * pts[i].x + pts[i].y * pts[i].y);
    }

    cv::Mat A_Trans;
    transpose(A, A_Trans);

    cv::Mat Inv_A;
    invert(A_Trans * A, Inv_A);

    const cv::Mat x = Inv_A * A_Trans * b;

    lsc.center.x = x.at<float>(0, 0);
    lsc.center.y = x.at<float>(1, 0);
    lsc.radius = (float)sqrt(lsc.center.x * lsc.center.x + lsc.center.y * lsc.center.y + x.at<float>(2, 0));

    const int lr = 1;				// 学习率, 一般拟合的点都比较多, 所以我设置的比较大, 可以根据你的情况来设置
    const int iters = pt_len;		// 迭代次数, 我设置成了有多少个点就迭代多少次, 也可以根据实际情况设置

    std::vector<float> losses(pt_len);	// 每次迭代后的 loss 值
    std::vector<float> min_loss(pt_len);	// 每次迭代后的最小 loss
    std::vector<float> root_val(pt_len);	// 每次迭代中的开平方值, 方便以后使用

    for (int i = 0; i < iters; i++)
    {
        float loop_loss = 0;

        for (int j = 0; j < pt_len; j++)
        {
            // 这里第一次迭代的 x, y, r 是最小二乘的结果, 第二次迭代开始就是修正后的结果
            root_val[j] = sqrt((pts[j].x - lsc.center.x) * (pts[j].x - lsc.center.x) +
                (pts[j].y - lsc.center.y) * (pts[j].y - lsc.center.y));

            const float loss = root_val[j] - lsc.radius;

            losses[j] = loss;
            loop_loss += fabs(loss);
        }

        min_loss[i] = loop_loss;

        // 如果 loss 值不再减小, 就提前结束
        if (i > 0 && min_loss[i] > min_loss[i - 1])
        {
            break;
        }

        // 下面三个是梯度值
        float gx = 0;
        float gy = 0;
        float gr = 0;

        for (int j = 0; j < pt_len; j++)
        {
            // 在计算梯度时要先计算偏导数, 再将 x 代数公式得到
            float gxi = (lsc.center.x - pts[j].x) / root_val[j];

            if (losses[j] < 0)
            {
                gxi *= (-1);
            }

            float gyi = (lsc.center.y - pts[j].y) / root_val[j];

            if (losses[j] < 0)
            {
                gyi *= (-1);
            }

            float gri = -1;

            if (losses[j] < 0)
            {
                gri = 1;
            }

            gx += gxi;
            gy += gyi;
            gr += gri;
        }

        gx /= pt_len;
        gy /= pt_len;
        gr /= pt_len;

        lsc.center.x -= (lr * gx);
        lsc.center.y -= (lr * gy);
        lsc.radius -= (lr * gr);
    }

    if (adj_count < 1)
    {
        adj_count++;

        const cv::Point2f pt(lsc.center.x, lsc.center.y);

        for (auto it = pts.begin(); it != pts.end();)
        {
            const double dist = norm(*it - pt);

            if (dist < lsc.radius * 0.90 || dist > lsc.radius * 1.10)
            {
                it = pts.erase(it);
            }
            else
            {
                it++;
            }
        }

        goto LABEL_ADJ;
    }

    return lsc;
}

调用示例

int main()
{
    cv::Mat cv_src = cv::imread("58.jpg");

    cv::Mat cv_gray;
    cv::cvtColor(cv_src, cv_gray, cv::COLOR_BGR2GRAY);

    cv::Mat cv_canny;
    cv::Canny(cv_gray, cv_canny, 50, 150);

    cv::imshow("src", cv_src);
    cv::imshow("canny", cv_canny);

    std::vector<cv::Point> largest_contour;
    std::vector<std::vector<cv::Point>> contours;
    cv::findContours(cv_canny, contours, cv::RETR_LIST, cv::CHAIN_APPROX_SIMPLE);

    double max_area = 0.0;

    for (size_t i = 0; i < contours.size(); ++i)
    {
        double area = cv::contourArea(contours[i]);
        if (area > max_area)
        {
            max_area = area;
            largest_contour = contours[i];
        }
    }

    std::vector<cv::Point2f> contours_float;
    for (const cv::Point& p : largest_contour)
    {
        contours_float.push_back(cv::Point2f(p.x, p.y));
    }

    Circle3f c = circle_least_squares(contours_float);

    cv::circle(cv_src, c.center, 2, cv::Scalar(255, 0, 255), 4, cv::LINE_AA);
    cv::circle(cv_src, c.center, c.radius, cv::Scalar(0, 0, 255), 4, cv::LINE_AA);
    cv::imshow("circle", cv_src);
    cv::waitKey();
}