最小覆盖圆问题·三点定圆

问题描述

        给定一组点Point[n]，n > 0，求覆盖这组点的最小圆O(x,y,R);

定理1（三点定圆）：一组点Point[n]的最小覆盖圆最多由3不相同的点确定.

论证：

（1）当 n = 1，最小覆盖圆O(x，y， R) = (Point[0].x， Point[0].y，0)；

（2）当 n = 2，最小覆盖圆O(x，y， R) = ((Point[0].x + Point[1].x)/2 ， (Point[0].y + Point[1].y)/2，dis(Point[0], Point[1])/2)；

（3）当 n = 3，分类讨论：

        a.如果Point[0]、Point[1]、Point[3]三点构成钝角三角形，则最小覆盖原由构成最大值直线距离的两点确定；

证明：

1.必须满足：最小覆盖圆的直径L >= 钝角三角形的斜边a, 这样才能保证圆能包住三点

2 若将最小覆盖圆的直径与钝角三角形的斜边重合，依据定律1，斜边所对应的顶点一定在圆内，找到一个覆盖圆。

联合1，2得证。

b. 如果Point[0]、Point[1]、Point[3]三点构成锐角三角形，最小覆盖原由三点确定。（这里三个点不相同，否则属于（1）、（2）范畴）

1.若以最长边与圆直径重合，依据定律2，必然存在一点在圆外，故：最小覆盖圆直径L > 三角形最大边的边长a

2.以三点带入圆方差并联立解方程组必有唯一解，所确定的圆记为C

3.尝试将C缩小一点，若圆心不变得到C1，3点必然在圆外，如果为了包住3点向某个方向移动圆心，必然会远离某个点，该点必然在圆外

故C为最小覆盖圆

（4）当 n > 3, 如果处于最小覆盖圆的边界上的点超过3个，从中挑选3个点即可利用（1）（2）（3）确定一个最小覆盖圆。

定理2（边界点判定法制）：如果Point[i]在Point[1]...Point[i - 1]确定的最小覆盖圆O[i-1]的外部（且即O[i] != O[i-1]），那么Point[i]一定在Point[1]...Point[i]确定的最小覆盖圆O[i]的边界上。

论证(反证法)：

  条件A：O[i-1]是Point[1]...Point[i - 1]确定的最小覆盖圆，O[i]是Point[1]...Point[i]确定的最小覆盖圆，Point[i]在O[i-1]的外部（且O[i] != O[i-1]）。

  条件B：Point[i]一定在O[i]的边界上。

        假设Point[i]不在O[i]的边界上，那么Point[i]在O[i]的外部或内部（非B）。如果Point[i]在O[i]的外部，则O[i]不是Point[1]...Point[i]的最小覆盖圆，推出非A；如果Point[i]在O[i]的内部，由于外接圆是由其边界上的点确定的，那么确定O[i]的点来自于Point[1]...Point[i - 1]，而Point[1]...Point[i - 1]已经确定最小覆盖圆O[i-1]，那么O[i] = O[i-1]，则Point[i]不在O[i-1]外部，推出非A。因此A能推出B。

公式1(三点定圆)：

（未完待续....）

 三点定圆算法

       ① 以Point[1]为初始点确定一个圆O(1):O{Point[1]}。

       ②加入Point[i]，确定一个最小覆盖原O(i):O{Point[1]...Point[i]}（i = 2...n），确定方式如下：

                （1）如果Point[i] 在最小覆盖原O(i):O{Point[1]...Point[i-1]} 的内部或边界上，则O(i) = O(i-1)。当 i = n 时结束返回O(n)结束，否则继续加入下一个点（i = i+1 回到②）；

                  (2)  如果Point[i]在最小覆盖原O(i):O{Point[1]...Point[i-1]} 的外部，则Point[i]一定在O(i)的边界上，此时得到三点定圆的一个点，如果三点定圆点数量达到3，通过这三个点算出O（i）返回并结束，否则接着寻找下一个三点定圆点执行③

       ③将Point = {Point[1]，...，Point[n]} 更新为Point =  {Point[i]，Point[1]，...Point[i-1]，Point[i +1]，...，Point[n]}

编程

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#include<opencv2/opencv.hpp>
#define N 100010
using namespace std;
using namespace cv;
struct node { double x, y; }b[N];
node O;
double R;
double sqr(double x) { return x * x; }
double dis(node x, node y)
{
    return sqrt(sqr(x.x - y.x) + sqr(x.y - y.y));
}
bool incircle(node x)
{
    if (dis(O, x) <= R) return true;
    return false;
}
node solve(double a, double b, double c, double d, double e, double f)
{
    double y = (f * a - c * d) / (b * d - e * a);
    double x = (f * b - c * e) / (a * e - b * d);
    return (node{ x, y });
}
int main()
{
    Mat src = Mat(Size(1000, 1000), CV_8UC3,Scalar(0,0,0));
    vector<Point> pts;
    pts.push_back(Point(0, 0));
    pts.push_back(Point(501, 188));
    pts.push_back(Point(456, 207));
    pts.push_back(Point(678, 456));
    pts.push_back(Point(437, 111));
    pts.push_back(Point(193, 209));
    pts.push_back(Point(178, 496));
    pts.push_back(Point(137, 411));
    pts.push_back(Point(193, 309));
    int n = pts.size() - 1;
    for (int i = 1; i < pts.size(); i++) {
        b[i].x = pts[i].x;
        b[i].y = pts[i].y;
        circle(src, Point(b[i].x, b[i].y), 2, Scalar(125,255,0), 2);
        char str[5];
        sprintf_s(str, "%d", i);
        putText(src, str, pts[i], 0, 1, Scalar(0, 255, 0));
    }
    //random_shuffle(b + 1, b + n + 1);
    R = 0;
    for (int i = 1; i <= n; i++) {
        if (!incircle(b[i]))
        {
            O.x = b[i].x; O.y = b[i].y; R = 0;
            for (int j = 1; j < i; j++) {
                if (!incircle(b[j]))
                {
                    O.x = (b[i].x + b[j].x) / 2;
                    O.y = (b[i].y + b[j].y) / 2;
                    R = dis(O, b[i]);
                    for (int k = 1; k < j; k++)
                        if (!incircle(b[k]))
                        {
                            O = solve(
                                b[i].x - b[j].x, b[i].y - b[j].y, (sqr(b[j].x) + sqr(b[j].y) - sqr(b[i].x) - sqr(b[i].y)) / 2,
                                b[i].x - b[k].x, b[i].y - b[k].y, (sqr(b[k].x) + sqr(b[k].y) - sqr(b[i].x) - sqr(b[i].y)) / 2
                            );
                            R = dis(b[i], O);
                        }
                }
            }
        }
  
    }
    printf("%.10lf\n%.10lf %.10lf", R, O.x, O.y);
    circle(src, Point(O.x, O.y), R, Scalar(255, 55, 115), 1);
    imshow("result", src);
    waitKey(0);
    return 0;
}

 参考：

(51条消息) 最小圆覆盖算法_commonc的博客-CSDN博客_最小覆盖圆

三角形外接圆半径公式是什么?_百度知道 (baidu.com)