最小二乘法拟合圆c语言,最小二乘法拟合圆公式推导及其实现

最小二乘法拟合圆是一种常用的曲线拟合方法，可以用于将一组点拟合成一个圆。其原理是将点到圆心的距离的平方作为误差，通过最小化误差来求解圆心坐标和半径。

以下是最小二乘法拟合圆的公式推导及其在C语言中的实现：

1. 公式推导

设圆的方程为：(x-a)^2 + (y-b)^2 = r^2

将每个点代入圆的方程中，得到如下方程组：

(x1-a)^2 + (y1-b)^2 = r^2
(x2-a)^2 + (y2-b)^2 = r^2
...
(xn-a)^2 + (yn-b)^2 = r^2

将圆心坐标和半径表示为未知数，即：

A = [a,b,r]

则上述方程组可以表示为：

f(A) = [ (x1-a)^2 + (y1-b)^2 - r^2,
(x2-a)^2 + (y2-b)^2 - r^2,
...
(xn-a)^2 + (yn-b)^2 - r^2 ]

误差函数为：

E(A) = f(A)^T * f(A)

对误差函数求导，得到：

dE/dA = -2 * J^T * f(A)

其中，J为雅可比矩阵，表示误差函数对未知数的偏导数，即：

J = [ -2*(x1-a), -2*(y1-b), -2*r;
-2*(x2-a), -2*(y2-b), -2*r;
...
-2*(xn-a), -2*(yn-b), -2*r ]

将雅可比矩阵和方程组代入上式，得到：

dE/dA = -2 * ( J^T * J * A + J^T * f(A) )

令dE/dA = 0，解得：

A = ( J^T * J )^-1 * J^T * f(A)

即可求出圆心坐标和半径。

2. C语言实现

在C语言中，可以使用矩阵运算库（如Eigen）来实现矩阵的乘法、转置和求逆等操作。以下是一个简单的C语言实现：

#include <Eigen/Dense>
#include <vector>

using namespace Eigen;
using namespace std;

// 最小二乘法拟合圆
// 输入：points-点集，格式为vector<Vector2d>
// 输出：圆心坐标和半径，格式为Vector3d
Vector3d fitCircle(const vector<Vector2d>& points)
{
    int n = points.size();

    // 构建J和f(A)矩阵
    MatrixXd J(n, 3);
    VectorXd f(n);
    for (int i = 0; i < n; i++) {
        double x = points[i](0);
        double y = points[i](1);
        J(i, 0) = -2 * x;
        J(i, 1) = -2 * y;
        J(i, 2) = -2;
        f(i) = -(x * x + y * y);
    }

    // 求解圆心坐标和半径
    Vector3d A = (J.transpose() * J).inverse() * J.transpose() * f;
    A(0) = -A(0) / 2;
    A(1) = -A(1) / 2;
    A(2) = sqrt(A(0) * A(0) + A(1) * A(1) - A(2));

    return A;
}

使用示例：

vector<Vector2d> points;
points.push_back(Vector2d(0, 1));
points.push_back(Vector2d(1, 0));
points.push_back(Vector2d(0, -1));
points.push_back(Vector2d(-1, 0));

Vector3d A = fitCircle(points);
cout << "圆心坐标：" << A.head<2>().transpose() << endl;
cout << "半径：" << A(2) << endl;