空间圆上N点拟合圆的算法

1、实际问题

在空间圆上测量出N个点（如下图）的三维坐标（x,y,z），即平面x，y坐标加上高程z，根据这N组点的坐标求出圆心的坐标和半径。

2、实现方法

利用用最小二乘法拟合出球面方程和平面方程，二者相交即为所求圆曲线。

2.1 球面方程解算

球面方程为：

(x - a)^2  + (y - b)^2 + (z - c)^2 = R^2

即： x^2 + y^2 + z ^2 + a^2 +b^2 +c^2 - 2ax - 2by -2cz = R ^2

令： A = 2a  ; B =2b ; C =2c  D=a^2 +b^2 +c^2 -R^2 

得： x^2 + y^2 + z ^2 - Ax -By - Cz + D =0

即：

Ax +By +Cz -D = x^2 +y^2 +z^2

设N组点的坐标数据为(x1,y1,z1)、(x2,y2,z2)、...(xn,yn,zn)

则用矩阵表示如下：

令：

   ，  ，

则转成：

得：

解算出A、B、C、D之后，可得：

2.2 平面方程解算

平面方程为： 

Ex + Fy + Gz = 1

用矩阵表示如下：

解算过程同球面方程，可求出E、F、G。

如此便可得到球面方程、球心坐标、球半径和平面方程了。

2.3 拟合圆圆心及半径解算

如下图所示（仅示意）：

圆心A为球心O到平面E的垂足

设平面方程为 ：

          Ex + Fy + Gz - 1  =0                   ①   

则平面的法向量为 （E,F,G）

球心O的坐标为(a,b,c) ，设垂足即圆心A坐标为 (x', y', z')，因为球心O到圆心A的连线与法向量是平行的，可以得到如下：

(x' -a )/E = (y' -b)/F = (z' - c)/G  = t            ②

由 ②得 x' = Et + a； y' = Ft + b；z' = Gt + c；

带入① 可以得到：

即可求出最终拟合圆的圆心坐标(x' , y' , z')。

球心O到圆心A的距离可由平面外一点到平面的距离公式求得：

d=abs(Ex' + Fy' + Gz' -1) / sqrt(E^2 + F^2 +G^2)

则拟合圆的半径r为:

r = sqrt(R^2 -d^2)

以上便是解决问题的理论基础，后续便可以用代码来实现了。