已知空间三点的三维坐标，求这三个点所确定的空间圆的圆心坐标和半径。
这是测绘领域经常会碰到的问题，今天示例关于这个空间三点求圆心程序的实现过程。
1、基础理论——空间三点确定圆心
已知空间三点的坐标为（x1,y1,z1），(x2,y2,z2)，(x3,y3,z3)，求这三个点所确定的空间圆的圆心坐标和半径。
根据圆心到三点的距离都为半径可列出下列三式：
由（1）=（2）得：
由（1）=（3）得：
联立（4）、（5）、（6）可得：
2、界面设计
所用Visual Studio为2015版本，调用控件构建如下界面，
主要通过手输坐标进行单次计算，如果想通过读取文本数据批量计算，可以结合我之前的文章：txt文本读取、保存的代码实现【vb.Net】。
用一组数据测试的程序运行结果如下：
3、程序实现
本程序的整个工程文件资源已经上传，想用的话可以在这里下载（免费）：
重点主要有两个，
（1）理论公式的实现，这部分的代码如下：
Public Function qiuyuanxin(ByVal p1(,) As Double, ByVal p2(,) As Double, ByVal p3(,) As Double, ByRef yuanxin(,) As Double) As Double
'------------------空间三点求圆心
Dim x1 As Double, x2 As Double, x3 As Double
Dim y1 As Double, y2 As Double, y3 As Double
Dim z1 As Double, z2 As Double, z3 As Double
x1 = p1(1, 1)
y1 = p1(1, 2)
z1 = p1(1, 3)
x2 = p2(1, 1)
y2 = p2(1, 2)
z2 = p2(1, 3)
x3 = p3(1, 1)
y3 = p3(1, 2)
z3 = p3(1, 3)
Dim A1 As Double, B1 As Double, C1 As Double, D1 As Double
A1 = y1 * z2 - y1 * z3 - z1 * y2 + z1 * y3 + y2 * z3 - y3 * z2
B1 = x1 * z3 - x1 * z2 + z1 * x2 - z1 * x3 - x2 * z3 + x3 * z2
C1 = x1 * y2 - x1 * y3 - y1 * x2 + y1 * x3 + x2 * y3 - x3 * y2
D1 = x1 * y3 * z2 - x1 * y2 * z3 + x2 * y1 * z3 - x3 * y1 * z2 - x2 * y3 * z1 + x3 * y2 * z1
Dim A2 As Double
Dim B2 As Double
Dim C2 As Double
Dim D2 As Double
A2 = (x2 - x1) * 2
B2 = (y2 - y1) * 2
C2 = (z2 - z1) * 2
D2 = x1 ^ 2 + y1 ^ 2 + z1 ^ 2 - x2 ^ 2 - y2 ^ 2 - z2 ^ 2
Dim A3 As Double
Dim B3 As Double
Dim C3 As Double
Dim D3 As Double
A3 = (x3 - x1) * 2
B3 = (y3 - y1) * 2
C3 = (z3 - z1) * 2
D3 = x1 ^ 2 + y1 ^ 2 + z1 ^ 2 - x3 ^ 2 - y3 ^ 2 - z3 ^ 2
Dim a(3, 3) As Double
Dim B(,) As Double
Dim C(3, 3) As Double
Dim D(3, 1) As Double
Dim m(3, 1) As Double
D(1, 1) = D1
D(2, 1) = D2
D(3, 1) = D3
a(1, 1) = A1
a(1, 2) = B1
a(1, 3) = C1
a(2, 1) = A2
a(2, 2) = B2
a(2, 3) = C2
a(3, 1) = A3
a(3, 2) = B3
a(3, 3) = C3
Call jzqn(a, B)
For i = 1 To 3
For j = 1 To 3
C(i, j) = B(i, j) * (-1)
Next j
Next i
For i = 1 To 3
For k = 1 To 3
m(i, 1) = m(i, 1) + C(i, k) * D(k, 1) '矩阵相乘
Next k
Next i
For i = 1 To 3
yuanxin(1, i) = m(i, 1)
Next i
'求半径
yuanxin(1, 4) = Math.Sqrt((x1 - yuanxin(1, 1)) ^ 2 + (y1 - yuanxin(1, 2)) ^ 2 + (z1 - yuanxin(1, 3)) ^ 2)
End Function（2）其中调用的矩阵求逆函数jzqn( )，这部分的代码如下：
Public Function jzqn(ByVal qa(，) As Double, ByRef na(，) As Double) As Double
'-------------------------------------------矩阵求逆
Dim a(，)
Dim zj As Double, xxx As Double, q As Integer
Dim n = UBound(qa, 1)
ReDim na(n, n)
ReDim a(n, 2 * n)
For i = 1 To n
For j = 1 To n
a(i, j) = qa(i, j)
Next j
Next i
For i = 1 To n
For j = n + 1 To 2 * n
If j - i = n Then
a(i, j) = 1
Else
a(i, j) = 0
End If
Next j
Next i
For i = 1 To n
If a(i, i) = 0 Then
For q = i To n
If a(q, i) <> 0 Then
For w = i To 2 * n
zj = a(i, w)
a(i, w) = a(q, w)
a(q, w) = zj
Next w
Exit For
End If
Next q
If q > n Then
MsgBox("此矩阵不可逆")
Exit Function
End If
End If
For k = 2 * n To i Step -1
a(i, k) = a(i, k) / a(i, i)
Next k
For j = i + 1 To n
If a(j, i) <> 0 Then
For k = 2 * n To i Step -1
a(j, k) = a(j, k) / a(j, i) - a(i, k)
Next k
End If
Next j
Next i
For i = n To 1 Step -1
If a(i, i) = 0 Then
For q = i - 1 To 1 Step -1
If a(q, i) <> 0 Then
For w = i To 2 * n
zj = a(i, w)
a(i, w) = a(q, w)
a(q, w) = zj
Next w
Exit For
End If
Next q
End If
For k = 2 * n To i Step -1
a(i, k) = a(i, k) / a(i, i)
Next k
For j = i - 1 To 1 Step -1
If a(j, i) <> 0 Then
xxx = a(j, i)
For k = 2 * n To 1 Step -1
a(j, k) = a(j, k) / xxx - a(i, k)
Next k
End If
Next j
Next i
For i = 1 To n
For j = 1 To n
na(i, j) = a(i, j + n)
Next j
Next i
End Function如果是习惯C++的朋友，矩阵求逆可以参考我之前写的这篇文章：N阶矩阵求逆的方法及代码实现【C++】
欢迎交流！
标签：Dim,vb,Double,Next,圆心,y1,Net,x1,z1 From： https://blog.csdn.net/qq_34451885/article/details/144454052