G56 最小圆覆盖 随机增量法【计算几何】

#include <iostream>
#include <cstring>
#include <algorithm>
#include <cmath>
#define x first
#define y second
using namespace std;

const int N=50010;
const double PI=acos(-1);
int n;
struct Point{double x,y;} p[N];
struct Circle{Point p; double r;} C;

Point operator+(Point a, Point b){
  return {a.x+b.x, a.y+b.y};
}
Point operator-(Point a, Point b){
  return {a.x-b.x, a.y-b.y};
}
Point operator*(Point a, double t){
  return {a.x*t, a.y*t};
}
Point operator/(Point a, double t){
  return {a.x/t, a.y/t};
}
double operator*(Point a, Point b){
  return a.x*b.y-a.y*b.x;
}
Point rotate(Point a, double b){
  return {a.x*cos(b)-a.y*sin(b),a.x*sin(b)+a.y*cos(b)};
}
double dis(Point a, Point b){
  return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
Point cross(Point a,Point u,Point b,Point v){ 
  double t=(a-b)*v/(v*u);
  return a+u*t; //直线交点
}
pair<Point, Point> midperp(Point a, Point b){
  return {(a+b)/2, rotate(b-a,PI/2)}; //中垂线
}
Circle cover(Point a, Point b){
  return {(a+b)/2, dis(a,b)/2}; //覆盖两点的圆
}
Circle cover(Point a, Point b, Point c){
  auto u=midperp(a, b), v=midperp(a, c);
  auto p=cross(u.x, u.y, v.x, v.y);
  return {p, dis(p, a)}; //覆盖三点的圆
}
void increment(){ //增量法
  C={p[1], 0};
  for(int i=2; i<=n; i++){
    if(C.r<dis(C.p, p[i])){
      C={p[i], 0}; //一点圆
      for(int j=1; j<i; j++){
        if(C.r<dis(C.p, p[j])){
          C=cover(p[i], p[j]); //两点圆
          for(int k=1; k<j; k++)
            if(C.r<dis(C.p, p[k]))
              C=cover(p[i],p[j],p[k]); //三点圆
        }
      }
    }
  }  
}
int main(){
  scanf("%d", &n);
  for(int i=1;i<=n;i++)scanf("%lf%lf",&p[i].x,&p[i].y);
  double a, b;
  scanf("%lf%lf", &a, &b);
  for(int i=1; i<=n; i++){
    p[i]=rotate(p[i], -a/180*PI);
    p[i].x/=b; //椭圆压缩成圆
  }
  random_shuffle(p+1, p+n+1); //随机化
  increment(); //增量法
  printf("%.3lf\n", C.r);
}

G56 最小圆覆盖 随机增量法【计算几何】

问题描述

给定平面上n个点，求一个半径最小的圆，使得所有点都在圆内或圆上。

算法思路

随机增量法。基本思想是：

1. 随机打乱点的顺序，避免特定输入顺序下的最坏情况。
2. 从空集开始，逐个加入点，维护当前的最小覆盖圆。
3. 当加入新点时，若该点在当前圆内则圆不变；否则，新点必在新圆边界，以该点为起点，逐步扩展圆，确保所有点被覆盖。

步骤

1. 初始化当前圆为空圆（半径0）。
2. 随机打乱点的顺序。
3. 遍历每个点p：
   • 若p在当前圆内，继续下一点。
   • 否则，以p为新圆起点（半径0），遍历之前所有点q：
     • 若q在新圆内，继续。
     • 否则，以p、q为直径构建新圆。
     • 遍历之前所有点r：
       • 若r在当前新圆内，继续。
       • 否则，以p、q、r三点确定外接圆为新圆。
   • 更新当前圆为新圆。
4. 最终圆即为最小覆盖圆。

#include <iostream>
#include <vector>
#include <cmath>
#include <algorithm>
#include <random>

using namespace std;

struct Point {
    double x, y;
    Point(double x = 0, double y = 0) : x(x), y(y) {}
};

double dist2(const Point& a, const Point& b) {
    double dx = a.x - b.x;
    double dy = a.y - b.y;
    return dx*dx + dy*dy;
}

double dist(const Point& a, const Point& b) {
    return sqrt(dist2(a, b));
}

struct Circle {
    Point center;
    double radius;
    Circle() : center(Point(0, 0)), radius(0) {}
    Circle(Point c, double r) : center(c), radius(r) {}
    Circle(Point a, Point b, Point c) {
        double A = b.x - a.x, B = b.y - a.y;
        double C = c.x - a.x, D = c.y - a.y;
        double E = A*(a.x + b.x) + B*(a.y + b.y);
        double F = C*(a.x + c.x) + D*(a.y + c.y);
        double G = 2*(A*(c.y - b.y) - B*(c.x - b.x));
        if (fabs(G) < 1e-8) { // 三点共线
            if (dist2(a, b) >= dist2(a, c) && dist2(a, b) >= dist2(b, c)) {
                center = Point((a.x + b.x)/2, (a.y + b.y)/2);
                radius = dist(a, b)/2;
            } else if (dist2(a, c) >= dist2(a, b) && dist2(a, c) >= dist2(b, c)) {
                center = Point((a.x + c.x)/2, (a.y + c.y)/2);
                radius = dist(a, c)/2;
            } else {
                center = Point((b.x + c.x)/2, (b.y + c.y)/2);
                radius = dist(b, c)/2;
            }
        } else {
            center.x = (D*E - B*F)/G;
            center.y = (A*F - C*E)/G;
            radius = dist(center, a);
        }
    }
    Circle(Point a, Point b) {
        center.x = (a.x + b.x)/2;
        center.y = (a.y + b.y)/2;
        radius = dist(a, b)/2;
    }
    Circle(Point a) : center(a), radius(0) {}
    bool contains(const Point& p) {
        return dist2(p, center) <= radius*radius + 1e-8;
    }
};

Circle minimal_circle_cover(vector<Point> points) {
    random_shuffle(points.begin(), points.end());
    Circle C;
    int n = points.size();
    for (int i = 0; i < n; ++i) {
        if (!C.contains(points[i])) {
            C = Circle(points[i]);
            for (int j = 0; j < i; ++j) {
                if (!C.contains(points[j])) {
                    C = Circle(points[i], points[j]);
                    for (int k = 0; k < j; ++k) {
                        if (!C.contains(points[k])) {
                            C = Circle(points[i], points[j], points[k]);
                        }
                    }
                }
            }
        }
    }
    return C;
}

int main() {
    int n;
    cin >> n;
    vector<Point> points(n);
    for (int i = 0; i < n; ++i) {
        cin >> points[i].x >> points[i].y;
    }
    Circle c = minimal_circle_cover(points);
    printf("Center: (%.2f, %.2f), Radius: %.2f\n", c.center.x, c.center.y, c.radius);
    return 0;
}