在过去的一周时间里一直在计算师兄交代的任务

，啥都没学，就一直在干这活。其中呢，我想的一个算法中就用到了本文中的这个小知识(主要就是线性代数和高中的一些知识结合起来)。师兄说做的这个东西的内容要保密，但是我想了一下，本篇的这个内容还是可以说的

，这不涉及我做的那些东西的内容。

我想大家应该都知道圆的一般方程吧：

这个平方项前面的系数可以不写A，写1也是一样的(除一下就行)。

比如说现在有圆上的三个点：

(x1, y1)

(x2, y2)

(x3, y3)

和圆上的任意点：

(x, y)

这个四个点满足上述的圆的一般方程，因此可以把这四个点带入方程中得到如下方程组：

学过线性代数的应该都有印象：齐次线性方程组要有非零解(这里要解 圆的一般方程中的系数ABCD)，就要求系数行列式等于0。(这里是求解ABCD这4个系数)所以有下式：

把这个左边的行列式按照第一行展开：

这不就是上面的那个圆的一般方程嘛

。

所以A, B, C, D 这四个系数不就求出来了嘛

，如下：

那么现在就已经知道了圆的一般方程了，怎么求圆心和半径呢

？当然是配方啊！配方结果如下：

好了，现在圆心和半径就知道了，如下：

公式都推导出来了，下面该写代码了！

计算圆心半径的函数

function [x_center, y_center, r] = Cal_Radius_Center(x, y)

%% 本函数根据三点坐标来计算圆心和半径

%% x和y都是1x3的列向量

%   实际上是x值和z值

one = ones(size(x));

% 圆的一般方程系数矩阵的行列式

% A*x^2 + A*y^2 + B*x + C*y + D = 0

A =  det([x, y, one]);

B = -det([x.^2 + y.^2, y, one]);

C =  det([x.^2 + y.^2, x, one]);

D = -det([x.^2 + y.^2, x, y]);

x_center = -B/(2*A);

y_center = -C/(2*A);

r = sqrt( (B^2+C^2-4*A*D) / (4*A^2) );

end

应用一下：

现在有-0.9~0.9范围内的若干个单位圆上的数据点，但是这些点加上了噪声

，我们来利用这些点来看一下平均圆心和半径与实际的圆心和半径之间的误差多大。

clc

clear

close all

% 我们来若干个点

x = -0.9:0.01:0.9;

y = sqrt(1 - x.^2);

% 再加点噪声

y = y + 0.1 .* sin(rand(size(x)));

% 画出这些散点看看

figure

scatter(x, y, 's', 'SizeData', 100)

% 拟合一下

p = polyfit(x, y, 6); % 6次(自己调整)

x_p = linspace(-0.9, 0.9, 600);

y_p = polyval(p, x_p);

hold on

plot(x_p, y_p, 'r', 'LineWidth', 3)

legend('原数据', '拟合后')

daspect([1 1 1])  % 让坐标轴采用相同的标尺

%% 开始计算这些散点所在的圆的圆心和半径(结果取平均值)

% 把x_p和y_p都分成三段，每次从其中各随机取出一个值凑成三组坐标用于计算

len = fix(length(x_p)/3);

x_selected_zone1 = x_p(1:len);

x_selected_zone2 = x_p(len+1:2*len);

x_selected_zone3 = x_p(2*len+1:3*len);

y_selected_zone1 = y_p(1:len);

y_selected_zone2 = y_p(len+1:2*len);

y_selected_zone3 = y_p(2*len+1:3*len);

x_centers = [];

y_centers = [];

radius = [];

figure

scatter(x, y, 's', 'SizeData', 100)

hold on

plot(x_p, y_p, 'r', 'LineWidth', 3)

while ~isempty(x_selected_zone1)

% 每从中抽取一个值后，把抽出来的这个值从数组中删除，取出x对应的y，并把取出来的y也赋空

x1 = x_selected_zone1( unidrnd( length(x_selected_zone1) ) );

y1 = y_selected_zone1( x_selected_zone1 == x1 );

y_selected_zone1( x_selected_zone1 == x1 ) = [];

x_selected_zone1( x_selected_zone1 == x1 ) = []; % x一定要放在后面赋空，因为z是要依据x来取的

x2 = x_selected_zone2( unidrnd( length(x_selected_zone2) ) );

y2 = y_selected_zone2( x_selected_zone2 == x2 );

y_selected_zone2( x_selected_zone2 == x2 ) = [];

x_selected_zone2( x_selected_zone2 == x2 ) = [];

x3 = x_selected_zone3( unidrnd( length(x_selected_zone3) ) );

y3 = y_selected_zone3( x_selected_zone3 == x3 );

y_selected_zone3( x_selected_zone3 == x3 ) = [];

x_selected_zone3( x_selected_zone3 == x3 ) = [];

x = [x1, x2, x3]';

y = [y1, y2, y3]';

% 计算半径和圆心

[x_center, y_center, r] = Cal_Radius_Center(x, y);

x_centers = [x_centers, x_center];

y_centers = [y_centers, y_center];

radius = [radius, r];

hold on

line([x_center, x1], [y_center, y1])

hold on

line([x_center, x2], [y_center, y2])

hold on

line([x_center, x3], [y_center, y3])

end

daspect([1 1 1])

%% 可以计算圆心中值点了和半径均值

x_centers_mean = mean(x_centers);

y_centers_mean = mean(y_centers);

radius_mean = mean(radius);

disp('最终平均圆心坐标和平均半径')

disp('===================================')

disp(['平均圆心坐标：(', num2str(x_centers_mean), ', ', num2str(y_centers_mean), ')'])

disp('')

disp(['平均半径：', num2str(radius_mean)])

disp('===================================')

结果：

(和圆心(0, 0) 以及半径1是很接近的)

每次运行的结果是不会相同的(因为是随机从三段中各取一个点来算的，最后再平均一下)，但是误差也很小。

--END--