對於網上常用的擬合圓程式碼(經過修改, 因為除數可能為0)

/*
* 參考: http://blog.csdn.net/liyuanbhu/article/details/50889951
* 通過最小二乘法來擬合圓的資訊
* pts: 所有點座標
* center: 得到的圓心座標
* radius: 圓的半徑
*/
static bool CircleInfo2(std::vector<cv::Point2d>& pts, cv::Point2d& center, double& radius)
{
    center = cv::Point2d(0, 0);
    radius = 0.0;
    if (pts.size() < 3) return false;;

    double sumX = 0.0;
    double sumY = 0.0;
    double sumX2 = 0.0;
    double sumY2 = 0.0;
    double sumX3 = 0.0;
    double sumY3 = 0.0;
    double sumXY = 0.0;
    double sumX1Y2 = 0.0;
    double sumX2Y1 = 0.0;
    const double N = (double)pts.size();
    for (int i = 0; i < pts.size(); ++i)
    {
        double x = pts.at(i).x;
        double y = pts.at(i).y;
        double x2 = x * x;
        double y2 = y * y;
        double x3 = x2 *x;
        double y3 = y2 *y;
        double xy = x * y;
        double x1y2 = x * y2;
        double x2y1 = x2 * y;

        sumX += x;
        sumY += y;
        sumX2 += x2;
        sumY2 += y2;
        sumX3 += x3;
        sumY3 += y3;
        sumXY += xy;
        sumX1Y2 += x1y2;
        sumX2Y1 += x2y1;
    }
    double C = N * sumX2 - sumX * sumX;
    double D = N * sumXY - sumX * sumY;
    double E = N * sumX3 + N * sumX1Y2 - (sumX2 + sumY2) * sumX;
    double G = N * sumY2 - sumY * sumY;
    double H = N * sumX2Y1 + N * sumY3 - (sumX2 + sumY2) * sumY;

    double denominator = C * G - D * D;
    if (std::abs(denominator) < DBL_EPSILON) return false;
    double a = (H * D - E * G) / (denominator);
    denominator = D * D - G * C;
    if (std::abs(denominator) < DBL_EPSILON) return false;
    double b = (H * C - E * D) / (denominator);
    double c = -(a * sumX + b * sumY + sumX2 + sumY2) / N;

    center.x = a / (-2);
    center.y = b / (-2);
    radius = std::sqrt(a * a + b * b - 4 * c) / 2;
    return true;
}

---------------------

本文來自 wfh2015 的CSDN 部落格 ,全文地址請點選:https://blog.csdn.net/wfh2015/article/details/79626581?utm_source=copy 

然而,很多人對這部分程式碼表示有疑問,所以我根據網上的資料另外寫了一個部分,但是速度差很多,僅僅用於學習程式碼公式是按照這個部落格來的。個人認為,這份部落格的大部分是正確的,比如求圓的圓心是正確的。但是最後面的計算圓的半徑部落格公式可能有問題,我根據下面網友評論修改了一下,得到正確的結果。下面是我寫的程式碼

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

/////////////////////////////////////根據公式編寫程式碼//////////////////////////////////////////
/**
計算所有點的X座標均值
pts: 點座標
return: 均值
*/
static double meanX(std::vector<cv::Point2d>& pts)
{
    if (pts.size() <= 0)    return 0.0;

    const double N = (double)pts.size();
    double sum = 0.0;
    for (int i = 0; i < pts.size(); ++i)
    {
        cv::Point2d pt = pts.at(i);
        sum += (pt.x);
    }

    double avg = sum / N;
    return avg;
}

/**
計算所有點的Y座標均值
pts: 點座標
return: 均值
*/
static double meanY(std::vector<cv::Point2d>& pts)
{
    if (pts.size() <= 0)    return 0.0;
    const double N = (double)pts.size();

    double sum = 0.0;
    for (int i = 0; i < pts.size(); ++i)
    {
        cv::Point2d pt = pts.at(i);
        sum += (pt.y);
    }

    double avg = sum / N;
    return avg;
}

static double Ui(std::vector<cv::Point2d>& pts, int index, double meanXValue)
{
    double xi = pts.at(index).x;
    return xi - meanXValue;
}

static double Vi(std::vector<cv::Point2d>& pts, int index, double meanYValue)
{
    double yi = pts.at(index).y;
    return yi - meanYValue;
}

static double Suvv(std::vector<cv::Point2d>& pts, double mean_x, double mean_y)
{
    double sum = 0.0;
    for (int i = 0; i < pts.size(); ++i)
    {
        double u = Ui(pts, i, mean_x);
        double v = Vi(pts, i, mean_y);
        double cur = u * v *v;
        sum += cur;
    }
    return sum;
}

static double Suuv(std::vector<cv::Point2d>& pts, double mean_x, double mean_y)
{
    double sum = 0.0;
    for (int i = 0; i < pts.size(); ++i)
    {
        double u = Ui(pts, i, mean_x);
        double v = Vi(pts, i, mean_y);
        double cur = u * u * v;
        sum += cur;
    }
    return sum;
}

static double Suv(std::vector<cv::Point2d>& pts, double mean_x, double mean_y)
{
    double sum = 0.0;
    for (int i = 0; i < pts.size(); ++i)
    {
        double u = Ui(pts, i, mean_x);
        double v = Vi(pts, i, mean_y);
        double cur = u * v;
        sum += cur;
    }
    return sum;
}

static double Svv(std::vector<cv::Point2d>& pts, double mean_y)
{
    double sum = 0.0;
    for (int i = 0; i < pts.size(); ++i)
    {
        double value = Vi(pts, i, mean_y);
        double cur = value * value;
        sum += cur;
    }
    return sum;
}

static double Suu(std::vector<cv::Point2d>& pts, double mean_x)
{
    double sum = 0;
    for (int i = 0; i < pts.size(); ++i)
    {
        double value = Ui(pts, i, mean_x);
        double cur = value * value;
        sum += cur;
    }
    return sum;
}

static double Svvv(std::vector<cv::Point2d>& pts, double mean_y)
{
    double sum = 0.0;
    for (int i = 0; i < pts.size(); ++i)
    {
        double v = Vi(pts, i, mean_y);
        double cur = v * v * v;
        sum += cur;
    }
    return sum;
}

static double Suuu(std::vector<cv::Point2d>& pts, double mean_x)
{
    double sum = 0.0;
    for (int i = 0; i < pts.size(); ++i)
    {
        double u = Ui(pts, i, mean_x);
        double cur = u * u * u;
        sum += cur;
    }
    return sum;
}

static double Uc(double A, double B, double C, double D, double E, double F, double G)
{
    double numerator = A * B - C * D - E * D + B * F;
    double denominator = 2 * (B * B - G * D);
    double result = numerator / denominator;
    return result;
}

static double Vc(double A, double B, double C, double D, double E, double F, double G)
{
    double numerator = -1 * G * A + C * B + B * E - G * F;
    double denominator = 2 * (B * B - G * D);
    double result = numerator / denominator;
    return result;
}

static double GetRadius(std::vector<cv::Point2d>& pts, cv::Point2d center)
{
    double sum = 0.0;
    for (int i = 0; i < pts.size(); ++i)
    {
        double p1 = pts.at(i).x - center.x;
        double p2 = pts.at(i).y - center.y;
        double cur = p1*p1 + p2*p2;
        sum += cur;
    }
    double radius = std::sqrt(sum / (pts.size()));
    return radius;
}

static void CircleInfo(std::vector<cv::Point2d>& pts, cv::Point2d& center, double& radius)
{
    double mean_x = meanX(pts);
    double mean_y = meanY(pts);

    double A = Suuv(pts, mean_x, mean_y);
    double B = Suv(pts, mean_x, mean_y);
    double C = Suuu(pts, mean_x);
    double D = Svv(pts, mean_y);
    double E = Suvv(pts, mean_x, mean_y);
    double F = Svvv(pts, mean_y);
    double G = Suu(pts, mean_x);

    double uc = Uc(A, B, C, D, E, F, G);
    double vc = Vc(A, B, C, D, E, F, G);
    center.x = uc + mean_x;
    center.y = vc + mean_y;

    radius = GetRadius(pts, center);
}

---------------------

本文來自 wfh2015 的CSDN 部落格 ,全文地址請點選:https://blog.csdn.net/wfh2015/article/details/79626581?utm_source=copy 

另外,附上所有的測試程式碼下載。  最後,感謝我參考的部落格作者的公式推導。謝謝!