已知空间N点坐标求圆心坐标,半径 文档

已知空间N点坐标求圆心坐标，半径

注意哦这里是求圆心不是球心哦

条件：已知空间N点坐标，格式如下求圆心坐标，半径

-33.386698 -12.312448 -2301.396442

-33.668120 -12.571431 -2300.390996

-33.838611 -12.774933 -2299.691688

-34.079235 -13.616660 -2298.326277

-34.254998 -13.441280 -2298.192657

-34.356542 -13.755525 -2297.796473 ........................................................

实现方法：用最小二乘法拟合出球面方程和平面方程，两者相交即为所求圆曲线球面方程：(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组数据如下形式

现在就是求A B C D 的问题了

具体步骤如下

平面方程 A' x + B' y + C' z + D' = 0

可化为 A x + B y + C z -1 =0

用矩阵表示如下

同样可以求出A B C

这样就有了球面方程球心坐标球半径和平面方程了如图所示（网络找的图，意思着看囧）

圆心为球心到平面的垂足（也就是平面外一点到平面的坐标问题）

半径为sqrt(球半径^2 - 球心到平面的距离^2)

设平面方程为 A X + B Y + C Z +D = 0 ①球心坐标(a,b,c)

平面外一点到平面的坐标问题：

则平面的法向量为（A ,B ,C ）

设垂足即圆心(x' y' z') 球心到圆心的连线与法向量是平行的

可以得到如下(x' -a )/A = (y' -b)/B = (z' - c)/C = t ②

由②得x' = At + a; y' = B t + b ; z' = C t + c;

带入①可以得到(x' , y' , z')

平面外一点到平面的距离问题

d = abs(Ax' + By' + C z' +D) / sqrt(A^2 + B^2 +C^2)

到此为止已经有了足够的理论知识，下面是代码分别用OPENCV 和C 实现

1// 最小二乘法拟合球.cpp : 定义控制台应用程序的入口点。

2//

3

4#include "stdafx.h"

5#include

6#include

7#include

8usingnamespace std;

9

10#include

11#include

12#include

13

14#ifdef DEBUG

15#pragma comment(lib,"cxcore200d.lib")

16#pragma comment(lib,"cv200d.lib")

17#pragma comment(lib,"highgui200d.lib")

18#else

19#pragma comment(lib,"cxcore200.lib")

20#pragma comment(lib,"cv200.lib")

21#pragma comment(lib,"highgui200.lib")

22

23#endif

24

25void fitRound(char*filepath,int n)

26{

27 ifstream file(filepath);

28//int n = 0; //数据组数

29

30//file>>n;

31/*

32 x ~ NX4

33 | x1 y1 z1 -1 |

34 |. .. .... .......... |

35 | .. ...... ..........|

36 |xn yn zn -1|

37*/

38int xsize =4*n*sizeof(double);

39double*x = (double*)malloc(xsize);

40/*

41 y ~ NX1

42 |x1^2 + y1^2 +z1^2 |

43 |.......................................|

44 |.......................................|

45 |xn^2+yn^2+zn^2 |

46*/

47int ysize = n *sizeof(double);

48double*y = (double*)malloc(ysize);

49

50/*