利用三个平面或者空间的点位坐标可以求出该平面圆或者球的圆心或者球心和半径
import numpy as np
from numpy.linalg import det
# 由圆上三点确定圆心和半径# INPUT
# p1 : - 第一个点坐标, list或者array 1x3
# p2 : - 第二个点坐标, list或者array 1x3
# p3 : - 第三个点坐标, list或者array 1x3
# 若输入1x2的行向量, 末位自动补0, 变为1x3的行向量
def points2circle(p1, p2, p3):p1 = np.array(p1)p2 = np.array(p2)p3 = np.array(p3)num1 = len(p1)num2 = len(p2)num3 = len(p3)# 输入检查if (num1 == num2) and (num2 == num3):if num1 == 2:p1 = np.append(p1, 0)p2 = np.append(p2, 0)p3 = np.append(p3, 0)elif num1 != 3:print('\t仅支持二维或三维坐标输入')return Noneelse:print('\t输入坐标的维数不一致')return None# 共线检查temp01 = p1 - p2temp02 = p3 - p2temp03 = np.cross(temp01, temp02)#计算两个向量（向量数组）的叉乘。叉乘返回的数组既垂直于a，又垂直于b。# 如果a,b是向量数组，则向量在最后一维定义。该维度可以为2，也可以为3. 为2的时候会自动将第三个分量视作0补充进去计算。temp = (temp03 @ temp03) / (temp01 @ temp01) / (temp02 @ temp02) # @装饰器的格式来写的目的就是为了书写简单方便# temp03 @ temp03中的@ 含义是数组中每个元素的平方之和if temp < 10**-6:print('\t三点共线, 无法确定圆')return Nonetemp1 = np.vstack((p1, p2, p3)) # 行拼接temp2 = np.ones(3).reshape(3, 1) # 以a行b列的数组形式显示mat1 = np.hstack((temp1, temp2)) # size = 3x4m = +det(mat1[:, 1:])n = -det(np.delete(mat1, 1, axis=1)) # axis=1相对于把每一行当做列来排列p = +det(np.delete(mat1, 2, axis=1))q = -det(temp1)temp3 = np.array([p1 @ p1, p2 @ p2, p3 @ p3]).reshape(3, 1)temp4 = np.hstack((temp3, mat1))# 使用 stack，可以将一个列表转换为一个numpy数组，当axis=0的时候，和 使用 np.array() 没有什么区别，# 但是当 axis=1的时候，那么就是对每一行进行在列方向上进行运算，也就是列方向结合，# 此时矩阵的维度也从（2,3）变成了（3,2）# hstack(tup) ，参数tup可以是元组，列表，或者numpy数组，返回结果为numpy的数组temp5 = np.array([2 * q, -m, -n, -p, 0])mat2 = np.vstack((temp4, temp5)) # size = 4x5A = +det(mat2[:, 1:])B = -det(np.delete(mat2, 1, axis=1))C = +det(np.delete(mat2, 2, axis=1))D = -det(np.delete(mat2, 3, axis=1))E = +det(mat2[:, :-1])pc = -np.array([B, C, D]) / 2 / Ar = np.sqrt(B * B + C * C + D * D - 4 * A * E) / 2 / abs(A)return pc, r# pc : - 圆心坐标, array 1x3# r : - 半径, 标量pc2, r2 = points2circle([1, 2, -1], [-2, 1, 2], [0, -3, -3])
print(pc2, r2)
# 调用示例1 - 平面上三个点
# pc1, r1 = points2circle([1, 2], [-2, 1], [0, -3])
# 调用示例2 - 空间中三个点
# pc2, r2 = points2circle([1, 2, -1], [-2, 1, 2], [0, -3, -3])运行结果如下：
文中有部分内容参考了 http://t.csdn.cn/ec8Ke
大家想了解更多详细的知识可以参考：http://t.csdn.cn/ec8Ke
