空间三点共圆的计算代码（Python版本）

• 已知空间三点坐标，计算圆心和半径

def calc_circle_center_and_radius(p1, p2, p3):

    x1 = p1.x
    y1 = p1.y
    z1 = p1.z

    x2 = p2.x
    y2 = p2.y
    z2 = p2.z

    x3 = p3.x
    y3 = p3.y
    z3 = p3.z

    a1 = (y1*z2 - y2*z1 - y1*z3 + y3*z1 + y2*z3 - y3*z2)
    b1 = -(x1*z2 - x2*z1 - x1*z3 + x3*z1 + x2*z3 - x3*z2)
    c1 = (x1*y2 - x2*y1 - x1*y3 + x3*y1 + x2*y3 - x3*y2)
    d1 = -(x1*y2*z3 - x1*y3*z2 - x2*y1*z3 + x2*y3*z1 + x3*y1*z2 - x3*y2*z1)

    a2 = 2 * (x2 - x1)
    b2 = 2 * (y2 - y1)
    c2 = 2 * (z2 - z1)
    d2 = x1*x1 + y1*y1 + z1*z1 - x2*x2 - y2*y2 - z2*z2

    a3 = 2 * (x3 - x1)
    b3 = 2 * (y3 - y1)
    c3 = 2 * (z3 - z1)
    d3 = x1*x1 + y1*y1 + z1*z1 - x3*x3 - y3*y3 - z3*z3

    x = -(b1*c2*d3 - b1*c3*d2 - b2*c1*d3 + b2*c3*d1 + b3*c1*d2 - b3*c2*d1) /\
        (a1*b2*c3 - a1*b3*c2 - a2*b1*c3 + a2*b3*c1 + a3*b1*c2 - a3*b2*c1)
    y = (a1*c2*d3 - a1*c3*d2 - a2*c1*d3 + a2*c3*d1 + a3*c1*d2 - a3*c2*d1) /\
        (a1*b2*c3 - a1*b3*c2 - a2*b1*c3 + a2*b3*c1 + a3*b1*c2 - a3*b2*c1)
    z = -(a1*b2*d3 - a1*b3*d2 - a2*b1*d3 + a2*b3*d1 + a3*b1*d2 - a3*b2*d1) /\
        (a1*b2*c3 - a1*b3*c2 - a2*b1*c3 + a2*b3*c1 + a3*b1*c2 - a3*b2*c1)

    ptCenter = Point(x, y, z)
    radius = math.sqrt((x1 - x)*(x1 - x) + (y1 - y)
                       * (y1 - y) + (z1 - z)*(z1 - z))

    return ptCenter, radius

• 计算圆的单位法向量

def calc_norm_3d(p1, p2, p3):

    n0 = (p2.y-p1.y)*(p3.z-p1.z)-(p3.y-p1.y)*(p2.z-p1.z)
    n1 = (p2.z-p1.z)*(p3.x-p1.x)-(p3.z-p1.z)*(p2.x-p1.x)
    n2 = (p2.x-p1.x)*(p3.y-p1.y)-(p3.x-p1.x)*(p2.y-p1.y)

    n = np.array([n0, n1, n2])
    n_norm = np.linalg.norm(n)
    n = n/n_norm
    return n

• 获取圆的参数表示

    n = calc_norm_3d(p1, p2, p3)

    c, r = calc_circle_center_and_radius(p1, p2, p3)

    d1 = np.array([p1.x-c.x, p1.y-c.y, p1.z-c.z])
    d2 = np.array([p2.x-c.x, p2.y-c.y, p2.z-c.z])
    d3 = np.array([p3.x-c.x, p3.y-c.y, p3.z-c.z])
    n1 = np.array([p1.x-c.x, p1.y-c.y, p1.z-c.z])
    n1_norm = np.linalg.norm(n1)
    n1 = n1/n1_norm

    n2 = np.cross(n, n1)

    cos_angle = d1.dot(d3)/(np.linalg.norm(d1)*np.linalg.norm(d3))
    theta = np.arccos(cos_angle)

    theta_list = np.linspace(0, theta, 10)

    # Points on the arc from p1 to p3
    for i in range(0, 10):
        x = c.x + r*n1[0]*np.cos(theta_list[i]) + r*n2[0]*np.sin(theta_list[i])
        y = c.y + r*n1[1]*np.cos(theta_list[i]) + r*n2[1]*np.sin(theta_list[i])
        z = c.z + r*n1[2]*np.cos(theta_list[i]) + r*n2[2]*np.sin(theta_list[i])