圆拟合理论推导和Matlab实现

最新推荐文章于 2025-10-19 15:11:09 发布

原创最新推荐文章于 2025-10-19 15:11:09 发布 · 置顶 · 3.1k 阅读

45 ·

CC 4.0 BY-SA版权

文章标签：

#matlab #线性代数 #算法

优化方法专栏收录该内容

2 篇文章

订阅专栏

这篇博客介绍了如何通过矩阵论的方法推导最佳二乘圆的解析解，并提供了Matlab实现代码。文章首先阐述了问题背景和理论推导，通过非线性变换将问题转化为线性方程组，然后利用广义逆求解。最后，给出了Matlab代码示例，展示了如何应用推导结果进行实际的圆拟合操作。

最佳二乘圆的推导和Matlab实现

1 理论推导

最近在学习的过程中需要使用最佳二乘圆和最佳二乘椭圆以及其他一些最佳二乘的二次曲线，恰巧又刚好学过一点矩阵论。刚好拿来练习一下。

首先搞清楚问题，最佳二乘拟合圆，即输入为散点集 ${(xi,yi)∣xi∈R,yi∈R}\{(x_i,y_i)|x_i\in \R,y_i\in \R\}$ ,定义新的向量 $x2...xn]T∈Rn}\{X|X=[x_1\ x_2...x_n]^T\in \R^n\}$ 和 $y2...yn]T∈Rn}\{Y|Y=[y_1\ y_2...y_n]^T\in \R^n\}$ 。我们都知道不共线的三点能够唯一确定一个圆（共线的三点确定的圆圈半径无穷大，圆心在无穷远点）。那么当上面的 $n = 3$ 时，我们能够唯一的确定圆圈（把直线归纳为特殊的圆）；当上面的 $n > 3$ 时，就存在可能我们新加入的点并不在原来的圆上。最小二乘的思想是，多个点求出的圆与输入散点集合的距离尽可能的近。这里的近一般是使用向量间的欧式距离进行描述。我们知道在二维平面，唯一确定一个圆需要三个参数即圆心横坐标，圆心纵坐标以及半径三个参数，即我们需要求出这三个参数即可。

然后开始我们正式的推导，我们在中学时代就已经知道，圆方程的一般形式是：
$x-x_0)^2+(y-y_0)^2=R^2$
将其展开为：
$x^2+y^2-2xx_0-2yy_0+x_0^2+y_0^2=R^2$
[数学从来没有简单或复杂，只有已知与未知。]接下来这一步是关键，目的在于把关于 $x_0,y_0,R$ 的非线性方程转换为线性方程，然后利用矩阵的工具进行求解。变换后的方程为:
$x^2+y^2+ax+by+c=0$
其中：
$\begin{cases} \begin{aligned} a &= -2x_0\\ b &=-2y_0\\ c &=x_0^2+y_0^2-R^2 \end{aligned} \end{cases}$
我们使用一个非线性变换，把关于 $x_0,y_0,R$ 的非线性方程变成了关于 $a, b, c$ 的线性方程。既然是线性方程我们就可以写出矩阵形式：
$\begin{bmatrix} x&y&1\\ \end{bmatrix}\times \begin{bmatrix} a\\ b\\ c \end{bmatrix}= -x^2-y^2$
这样对于每个点，我们都能得到一个线性方程。对于 $n$ 个点，我们将得到的方程写成矩阵的形式：
$\begin{bmatrix} x_1 & y_1 & 1 \\ x_2 & y_2 & 1 \\ \vdots & \vdots &\vdots\\ x_n & y_n & 1 \\ \end{bmatrix}\times \begin{bmatrix} a\\ b\\ c\\ \end{bmatrix}= \begin{bmatrix} -x_1^2-y_1^2\\ -x_2^2-y_2^2\\ \vdots\\ -x_n^2-y_n^2\\ \end{bmatrix}$
这个问题中，未知数的个数为3，独立方程的个数为 $n$ 。矩阵论中给出了形如 $A X = b$ 的方程组在独立方程个数多于未知数个数时的最佳二乘解为 $X=A^+b$ ，其中 $A^+$ 表示的是矩阵 $A$ 的M-P广义逆。

对于我们这个问题我们需要求下面这个矩阵的广义逆:
$A=\begin{bmatrix} x_1 & y_1 & 1 \\ x_2 & y_2 & 1 \\ \vdots & \vdots &\vdots\\ x_n & y_n & 1 \\ \end{bmatrix}$
对于拟合问题，即 $n > 3$ 时，矩阵 $A$ 是一个列满秩矩阵。矩阵论中给出了列满秩矩阵的计算公式为：
$A^+=(A^HA)^-A^H$
对于实矩阵， $A^H$ 表示其对称矩阵，则未知数构成的向量的解为：
$\begin{bmatrix} a\\ b\\ c \end{bmatrix}=(A^HA)^-A^H \begin{bmatrix} -x_1^2-y_1^2\\ -x_2^2-y_2^2\\ \vdots\\ -x_n^2-y_n^2\\ \end{bmatrix}=B^-C$
计算 $B$ 矩阵:
$B=A^HA=\begin{bmatrix} x_1 &x_2&\cdots&x_n \\ y_1 &y_2&\cdots&y_n \\ 1 &1&\cdots&1\\ \end{bmatrix}\times \begin{bmatrix} x_1 & y_1 & 1 \\ x_2 & y_2 & 1 \\ \vdots & \vdots &\vdots\\ x_n & y_n & 1 \\ \end{bmatrix}\\ =\begin{bmatrix} \sum_{i=1}^{n}x_i^2&\sum_{i=1}^nx_iy_i&\sum_{i=1}^nx_i\\ \sum_{i=1}^nx_iy_i&\sum_{i=1}^{n}y_i^2&\sum_{i=1}^ny_i\\ \sum_{i=1}^nx_i&\sum_{i=1}^ny_i&n \end{bmatrix}$
计算 $C$ 矩阵：
$C=A^H*\begin{bmatrix} -x_1^2-y_1^2\\ -x_2^2-y_2^2\\ \vdots\\ -x_n^2-y_n^2\\ \end{bmatrix}= \begin{bmatrix} x_1 &x_2&\cdots&x_n \\ y_1 &y_2&\cdots&y_n \\ 1 &1&\cdots&1\\ \end{bmatrix}\times \begin{bmatrix} -x_1^2-y_1^2\\ -x_2^2-y_2^2\\ \vdots\\ -x_n^2-y_n^2\\ \end{bmatrix} =-\begin{bmatrix} \sum_{i=1}^n(x_i^3+x_iy_i^2)\\ \sum_{i=1}^n(x_i^2y_i+y_i^3)\\ \sum_{i=1}^n(x_i^2+y_i^2)\\ \end{bmatrix}$

至此我们在具备矩阵论知识的基础上，推导了最佳拟合圆的解析解。

2 Matlab 实现

[我说你的理论没用，把代码拿出来试试。]在有了理论基础的前提下，代码就可以四两拨千斤。

最小二乘圆拟合函数：

%%020406081012141618202224262830323436384042444648505254565860626466687072747678%%
%@func:Fit the scatters by least-square circle;
%@inpt:Coordinate set (x,y);
%@oupt:The circle struct contains center and radius;
%@note:
%@date:2021/01/23 by xiyin,li in hust;
function [Circle] = circleFit(x,y)
xlength = length(x);
if(xlength ~= length(y) | xlength  < 3)
    warning('the length of x and y must be same and greater than 3');
else
    xx = x.*x;xy = x.*y;yy = y.*y;
    B = [sum(xx) sum(xy) sum(x);sum(xy) sum(yy) sum(y);sum(x) sum(y) xlength];
    C = [sum(xx.*x+yy.*x);sum(xx.*y+yy.*y);sum(xx+yy)];
    D = inv(B)*C;
    Circle.x0 = D(1)/2;
    Circle.y0 = D(2)/2;
    Circle.R = sqrt(Circle.x0*Circle.x0+Circle.y0*Circle.y0+D(3));
    theta = linspace(0.1*pi,1.9*pi,100)
    Cirlcex = Circle.x0 + Circle.R*cos(theta);
    Cirlcey = Circle.y0 + Circle.R*sin(theta);
    plot(x,y,'o');
    hold on
    plot(Cirlcex ,Cirlcey,'r','linewidth',1.6);
    legend('输入散点','最小二乘圆')
    axis equal
end
end

使用如下的例子进行测试：

clc,clear;
close all
Sensor012 = csvread('CircleSensor012.csv');
Sensor.num0.x = Sensor012(1,:);
Sensor.num0.y = Sensor012(2,:);
Circle = circleFit(Sensor.num0.x,Sensor.num0.y);