绘制平面若干点的最小覆盖圆

绘制平面若干点的最小覆盖圆

编程思路：
1、随机生成10个点，x,y坐标取值在[1,9]区间，求最小覆盖圆的半径，参见求平面中若干点的覆盖圆的最小半径
2、判断任意三点的关系，锐角三角形时，求外接圆的圆心坐标和半径，其它情况下求最远两个点的中点坐标，半径是最远两点距离之半；找出所有这些半径的最大值，即所有这些点的最小覆盖圆半径。
3、绘制散点图，绘制最小覆盖圆。
说明：三点构成锐角三角形时，覆盖圆圆心坐标采用两种方法进行计算（通过circle和cir两个函数），一种是通过sympy解方程的方法求得，另外一种是直接套公式计算，两者相互印证，如果出现偏差，程序将直接报错退出。

def circle(z):
    p1, p2, p3 = z
    x1, y1 = p1
    x2, y2 = p2
    x3, y3 = p3
    f1 = (x - x1) ** 2 + (y - y1) ** 2 - (x - x2) ** 2 - (y - y2) ** 2
    f2 = (x - x1) ** 2 + (y - y1) ** 2 - (x - x3) ** 2 - (y - y3) ** 2
    return sympy.solve([f1, f2], [x, y])


def cir(z):
    p1, p2, p3 = z
    x1, y1 = p1
    x2, y2 = p2
    x3, y3 = p3
    a = x1 ** 2 + y1 ** 2
    b = x2 ** 2 + y2 ** 2
    c = x3 ** 2 + y3 ** 2
    g = (y3 - y2) * x1 + (y1 - y3) * x2 + (y2 - y1) * x3
    x0 = ((b - c) * y1 + (c - a) * y2 + (a - b) * y3) / (2 * g)
    y0 = ((c - b) * x1 + (a - c) * x2 + (b - a) * x3) / (2 * g)
    return [x0, y0]


def judging_triangles(p1, p2, p3):
    x1, y1 = p1
    x2, y2 = p2
    x3, y3 = p3
    # 判断是否共线
    mark = 0  # 不共线
    if x1 == x2 and x1 == x3:
        mark = 1  # 共线
    elif x1 == x2 or x1 == x3:
        mark = 0
    elif (y1 - y2) / (x1 - x2) == (y1 - y3) / (x1 - x3):
        '''斜率相同共线'''
        mark = 1
    dis_squared12 = (x1 - x2) ** 2 + (y1 - y2) ** 2
    dis_squared13 = (x1 - x3) ** 2 + (y1 - y3) ** 2
    dis_squared23 = (x2 - x3) ** 2 + (y2 - y3) ** 2
    dis_squared = [dis_squared12, dis_squared13, dis_squared23]
    dis_squared = sorted(dis_squared)
    if mark == 1 or dis_squared[2] >= dis_squared[0] + dis_squared[1]:
        '''满足三点共线或直角三角形或钝角三角形条件'''
        radius = math.sqrt(dis_squared[-1]) / 2
        if dis_squared[-1] == dis_squared12:
            x0, y0 = (x1 + x2) / 2, (y1 + y2) / 2
        elif dis_squared[-1] == dis_squared13:
            x0, y0 = (x1 + x3) / 2, (y1 + y3) / 2
        else:
            x0, y0 = (x2 + x3) / 2, (y2 + y3) / 2
        return [radius, x0, y0]
    else:
        '''满足锐角三角形条件'''
        a = math.sqrt(dis_squared[0])
        b = math.sqrt(dis_squared[1])
        c = math.sqrt(dis_squared[2])
        p = (a + b + c) / 2
        radius = (a * b * c) / (4 * math.sqrt(p * (p - a) * (p - b) * (p - c)))
        p = [p1, p2, p3]
        xy = circle(p)
        verification = cir(p)
        x0, y0 = xy[x], xy[y]
        xv0, yv0 = verification
        if float(x0) != xv0 or float(y0) != yv0:
            print("Validation error!")
            exit()
        return [radius, x0, y0]


if __name__ == '__main__':
    import math
    import random
    import sympy
    import numpy as np
    import matplotlib.pyplot as plt

    x = sympy.Symbol('x')
    y = sympy.Symbol('y')
    point_list = [[random.randint(1, 9), random.randint(1, 9)] for _ in range(10)]
    print(point_list)
    result = []
    length = len(point_list)
    for i in range(length - 2):
        for j in range(i + 1, length - 1):
            for k in range(j + 1, length):
                result.append(judging_triangles(point_list[i], point_list[j], point_list[k]))
    result = sorted(result, key=lambda xx: xx[0])
    print(result[-1])
    x, y = [], []
    for p in point_list:
        x.append(p[0])
        y.append(p[1])
    plt.figure(figsize=(5, 5))
    plt.scatter(x, y) 
    ax = plt.gca()
    ax.set_xlim((-2, 12))
    ax.set_ylim((-2, 12))
    disk1 = plt.Circle((result[-1][1], result[-1][2]), result[-1][0], color='b', fill=False)
    ax.add_artist(disk1)
    plt.show()

D:\Python\study\venv\Scripts\python.exe D:/Python/study/20200725/外接圆坐标.py
[[7, 5], [9, 3], [7, 4], [7, 6], [1, 8], [7, 1], [6, 3], [8, 4], [5, 2], [3, 6]]
[4.730925575903064, 125/26, 135/26]

Process finished with exit code 0

程序修改：

from math import sqrt
from random import random
from time import time
import matplotlib.pyplot as plt    # 导入绘图库


def circle(z):   # 计算锐角三角形外接圆圆心坐标
    p1, p2, p3 = z
    [x1, y1], [x2, y2], [x3, y3] = p1, p2, p3
    a, b, c = x1 ** 2 + y1 ** 2, x2 ** 2 + y2 ** 2, x3 ** 2 + y3 ** 2
    g = (y3 - y2) * x1 + (y1 - y3) * x2 + (y2 - y1) * x3
    x0 = ((b - c) * y1 + (c - a) * y2 + (a - b) * y3) / (2 * g)
    y0 = ((c - b) * x1 + (a - c) * x2 + (b - a) * x3) / (2 * g)
    return [x0, y0]


def judge_triangle(p1, p2, p3):   # 判断任意三个点的位置关系，并返回最小覆盖圆半径和坐标
    [x1, y1], [x2, y2], [x3, y3] = p1, p2, p3
    # 判断是否共线
    mark = 0  # 不共线
    if x1 == x2 and x1 == x3:
        mark = 1  # 共线
    elif x1 == x2 or x1 == x3:
        mark = 0
    elif (y1 - y2) / (x1 - x2) == (y1 - y3) / (x1 - x3):
        # 斜率相同共线
        mark = 1
    dis_squared12 = (x1 - x2) ** 2 + (y1 - y2) ** 2
    dis_squared13 = (x1 - x3) ** 2 + (y1 - y3) ** 2
    dis_squared23 = (x2 - x3) ** 2 + (y2 - y3) ** 2
    dis_squared = [dis_squared12, dis_squared13, dis_squared23]
    dis_squared = sorted(dis_squared)
    if mark == 1 or dis_squared[2] >= dis_squared[0] + dis_squared[1]:
        # 满足三点共线或直角三角形或钝角三角形条件
        radius = sqrt(dis_squared[-1]) / 2
        if dis_squared[-1] == dis_squared12:
            x0, y0 = (x1 + x2) / 2, (y1 + y2) / 2
        elif dis_squared[-1] == dis_squared13:
            x0, y0 = (x1 + x3) / 2, (y1 + y3) / 2
        else:
            x0, y0 = (x2 + x3) / 2, (y2 + y3) / 2
        return [radius, x0, y0]
    else:
        # 满足锐角三角形条件
        a, b, c = sqrt(dis_squared[0]), sqrt(dis_squared[1]), sqrt(dis_squared[2])
        p = (a + b + c) / 2
        radius = (a * b * c) / (4 * sqrt(p * (p - a) * (p - b) * (p - c)))
        p = [p1, p2, p3]
        return [radius, circle(p)[0], circle(p)[1]]


if __name__ == '__main__':
    n = int(input("请输入点数：\n"))
    print("计算中......")
    t0 = time()
    point_list, Ax, Ay = [], [], []    # 把随机产生的所有点的x,y坐标分别保存到列表中，为后续绘制散点图作准备
    while True:   # 随机产生n个不重复的点，坐标取值范围0~10
        ran_x, ran_y = random()*10, random()*10
        if [ran_x, ran_y] not in point_list:   # 确保随机数不重复
            point_list.append([ran_x, ran_y])
            Ax.append(ran_x), Ay.append(ran_y)
        if len(point_list) == n:
            break
    length, pp = len(point_list), [0, 0, 0]
    for i in range(length - 2):
        for j in range(i + 1, length - 1):
            for k in range(j + 1, length):
                pp1 = judge_triangle(point_list[i], point_list[j], point_list[k])
                if pp1[0] > pp[0]:
                    pp = pp1
    plt.rcParams['font.sans-serif'] = ['SimSun']    # 显示中文必须要有此两行
    plt.rcParams['axes.unicode_minus'] = False
    plt.figure("绘制平面上若干点最小覆盖圆", figsize=(7, 7))   # 设定画布（绘制区域）尺寸
    plt.scatter(Ax, Ay), plt.scatter(pp[1], pp[2], c='r')   # 绘制散点图，最小覆盖圆圆心点用红色表示
    plt.title("最小覆盖圆圆心（红色点）坐标[{:.5f},{:.5f}]，半径{:.5f}".format(pp[1], pp[2], pp[0]))
    ax = plt.gca()  # 创建一个移动坐标轴的实例
    ax.set_xlim((-2, 12)), ax.set_ylim((-2, 12))   # 设置坐标轴的起始和终止刻度
    disk1 = plt.Circle((pp[1], pp[2]), pp[0], color='b', fill=False)   # 绘制圆
    ax.add_artist(disk1)  # 添加到轴
    print("耗时{:.3f}s".format(time() - t0))
    plt.show()   # 显示