计算几何
二维几何基础
1.输入时计算坐标一般是实数,在写程序时一般用精度较高的double类型,不用float类型;double类型读入时用lf%格式,输出是用f%格式;
2.浮点数运算会产生精度误差,所以需要设置一个eps(偏差值),一般取1e-8;
3.判断一个浮点数是否等于0,不能直接用 == 0来判断,需要用sgn()函数判断是都小于eps;判断;
判断两个浮点数是否相等时,也不能直接用 == 判断是否相等,而是要用dcmp()函数判断是否相等;
fabs(x)就是返回x的绝对值
const double pi = acos(-1, 0); //高精度圆周率
const double eps = 1e-8; //偏差值,有时用1e-10
int sgn(double x){
if(fabs(x) < eps) return 0;//判断x是否等于0
else return x<0? -1:1; //判断x是正数还是负数
}
点和向量
1.点
二维平面中的点坐标(x,y);
struct point{
double x,y;
point(){}
point(double x,double y):x(x),y(y){}//构造函数方便赋值
};
2.两点之间的距离
把两点看成直角三角形的斜边
1.用库函数hypot()计算
hypot函数头文件math.h或cmath
hypot(a,b)的返回值为double类型，相当于sqrt(a x a+b x b)
double Distance(Point A, Point B) { return hypot(A.x - B.x, A.y - B.y);}
2.用sqrt()函数计算
double Distance(Point A,Point B){
return sqrt((A.x - B.x) *(A.x - B.x) + (A.y - B.y) * (A.y - B.y));
}
3.向量运算
在struct point中,对向量运算重载运算符;返回一个结构体;
struct Point {
double x, y;
Point() {}
Point(double x, double y) : x(x), y(y) {}
Point operator+(Point B) { return Point(x + B.x, y + B.y); }
Point operator-(Point B) { return Point(x - B.x, y - B.y); }
Point operator*(double k) { return Point(x * k, y * k); }//放大k倍
Point operator/(double k) { return Point(x / k, y / k); }//缩小k倍
bool operator==(Point B) { return sgn(x - B.x) == 0 && sgn(y - B.y) == 0; }
};
点积和叉积
1.点积
A . B = |A| |B| cos;如果已知向量A和向量B,则不需要角度;
typedef Point Vector;
double Dot(Vector A, Vector B) { return A.x * B.x + A.y * B.y;}
点积的作用
1.判断A与B的夹角是钝角还是锐角;
2.求向量A的长度 / 向量A长度的平方
double vector_length(Vector A) { return sqrt(Dot(A, A));}
double vector_length_square(Vector A) { return Dot(A, A);}
3.求向量A和B的夹角大小
double Angle(Vector A, Vector B) {
return acos(Dot(A, B) / vector_length(A) / vector_length(B));
}
2.叉积
A x B = |A| |B| sin;角度表示向量A旋转到向量B所经过的夹角;
可以用右手定则进行判断正负;向量具有先后顺序,不同的顺序有不同的正负答案
double Cross(Vector A, Vector B) {
return A.x * B.y - A.y * B.x;
}
1.判断向量A,B的方向关系;
A x B > 0,B在A的逆时针; A x B < 0,B在A的顺时针;
A x B = 0,共线,方向相同或相反;
2.计算两向量构成的平行四边形的有向面积
3个点A,B,C,以A为公共点,得到向量B - A,C - A,所构成的面积如下
double Area(Point A, Point B, Point C) {
return Cross(B - A, C - A);
}
以B或C为公共点构成的平行四边形面积相同,正负不同;
3.计算3个点构成的三角形的面积
等于三个点构成的平行四边形的一半
double Area(Point A, Point B, Point C) {
return Cross(B - A, C - A) / 2;
}
4.向量旋转
求向量A逆时针旋转角度为rad后的向量
Vector Rotate(Vector A, double rad) {
return Vector(A.x * cos(rad) - A.y * sin(rad), A.x * sin(rad) + A.y * cos(rad));
}
求向量A顺时针旋转角度为rad后的向量
Point rotate(Point a, double angle){
return Point(a.x * cos(angle) + a.y * sin(angle), -a.x * sin(angle) * cos(angle));
}
特殊情况是旋转90度
逆时针旋转90度 : rotate(A, pi /2),返回 vector(-A.y, A.x);
顺时针旋转90度 : rotate(A, - pi /2),返回 vector(A.y, - A.x);
求单位法向量 ,即逆时针旋转 90度,然后取单位值
Vector Normal(Vector A) {
return Vector(-A.y / vector_length(A), A.x / vector_length(A));
}
5.用叉积检查两个向量是否平行或重合
bool Parallel(Vector A, Vector B) {
return sgn(Cross(A, B)) == 0;
}
6.ToLeftTest
//判断折线bc是不是向ab的逆时针方向（左边）转向凸包构造时将会频繁用到此公式
bool ToLeftTest(Point a, Point b, Point c){
return Cross(b - a, c - b) > 0;
}
点和线
1.一条直线可以用线上面两个点来表示;
- 点向式表示直线 : p = p0 + vt;由一个点p0 和方向向量 v决定;p0是直线上任意一点(x,y); v是方向向量,给定两个点A,B; 那么v = B - A;
如果t无限制,p就是直线;
如果t在[0,1]内,p就是A,B之间的线段;
如果p大于0,p就是射线;
- 斜截式 y = kx + b;
- 普通式 ax + by + c = 0;
//点向式(根据参数t来控制)
struct Line {
Point v, p;
Line(Point v, Point p) : v(v), p(p) {}
Point point(double t){
return v + (p - v) * t;
}
};
struct Line {
Point p1, p2;
Line() {}
//根据端点确定直线
Line(Point p1, Point p2) : p1(p1), p2(p2) {}
//根据一个点和倾斜角angel确定直线，0 <= angel < pi
Line(Point p, double angel) {
p1 = p;
if (sgn(angel - pi / 2) == 0)p2 = (p1 + Point(0, 1));
else p2 = p1 + Point(1, tan(angel));
}
//ax + by + c = 0
Line(double a, double b, double c) {
if (sgn(a) == 0) {
p1 = Point(0, -c / b);
p2 = Point(1, -c / b);
} else if (sgn(b) == 0) {
p1 = Point(-c / a, 0);
p2 = Point(-c / a, 1);
} else {
p1 = Point(0, -c / b);
p2 = Point(1, (-c - a) / b);
}
}
};
线段可以用直线表示,加以区分
typedef Line Segment;
- 点和直线的位置关系
用直线上两点p1 和 p2 与点p构成两个向量,用叉积的正负判断方向,从而得到位置关系;
//点和直线关系:1 在左侧;2 在右侧;0 在直线上。根据视线从p1向p2看的左右
int Point_line_relation(Point p, Line v) {
int c = sgn(Cross(p - v.p1, v.p2 - v.p1));
if (c < 0)return 1; //1：p在v的左边
if (c > 0)return 2; //2：p在v的右边
return 0; //0：p在v上
}
- 判断点是否在线上
先用叉积判断是否共线, 然后用点积看p和v的两个端点产生的角是否为180°;
// 点和线段的关系：0 点p不在线段v上；1 点p在线段v上。
bool on_segment(Point p, Point a, Point b) {
return sgn(Cross(p - a, p - b)) == 0 && sgn(Dot(p - a, p - b)) <= 0;
}
bool Point_on_seg(Point p, Line v) {
return sgn(Cross(p - v.p1, v.p2 - v.p1)) == 0 && sgn(Dot(p - v.p1, p - v.p2)) <= 0;
}
4.点到直线的距离
已知点p和直线v(p1,p2),求p到v的距离;首先用叉积先求p,p1,p2构成的平行四边形的面积,然后用面积除以平行四边形的底边长,也就是线段(p1,p2)的长度,就得到了平行四边形的高,即p到直线的距离
double Dis_point_line(Point p, Line v) {
return fabs(Cross(p - v.p1, v.p2 - v.p1)/ Distance(v.p1, v.p2));
}
- 点在直线上的投影
已知直线上两点p1,p2以及直线外的一点p,求投影点p0;
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-56QCXK5P-1643292628101)(C:\Users\14996\AppData\Roaming\Typora\typora-user-images\image-20220120205516703.png)]
k = |p0 - p1| / |p2 - p1|k就是线段p0p1 和 p2p1长度的比值
(p - p1) . (p2 - p1) = |p- p1| * |p2 - p1| * cos = |p0 - p1| *|p2 - p1|
| p0 - p1 | = ((p - p1) . (p2 - p1)) / | p2 - p1 |
k = |p0 - p1| / |p2 - p1| = ((p - p1) . (p2 - p1)) / | p2 - p1 |*| p2 - p1 |;
p0 = p1 + k * (p2 - p1) = p1 + ((p - p1) . (p2 - p1)) / | p2 - p1 |*| p2 - p1 | * (p2 - p1);
Point Point_line_proj(Point p, Line v) {
double k = Dot(v.p2 - v.p1, p - v.p1) / vector_length_square(v.p2 - v.p1);
return v.p1 + (v.p2 - v.p1) * k;
}
6.点关于直线的对称点
先求一个点p对一条直线v的镜像点,先求点p在直线上的投影q,再求对称点p’
Point Point_line_symmetry(Point p, Line v) {
Point q = Point_line_proj(p, v);
return Point(2 * q.x - p.x, 2 * q.y - p.y);
}
7.点到线段的距离
点p到线段AB的距离等于在以下三段距离中去最小
1.从p出发做AB的垂线,如果交点在AB线段上,这个距离就是最小值
2.p到A的距离 3.p到B的距离
double Dis_point_seg(Point p, Segment v) {
if (sgn(Dot(p - v.p1, v.p2 - v.p1)) < 0 || sgn(Dot(p - v.p2, v.p1 - v.p2)) < 0) //点的投影不在线段上
return min(Distance(p, v.p1), Distance(p, v.p2));
return Dis_point_line(p, v); //点的投影在线段上
}
8.两条直线的位置关系
int Line_relation(Line v1, Line v2) {
if (sgn(Cross(v1.p2 - v1.p1, v2.p2 - v2.p1)) == 0) {
if (Point_line_relation(v1.p1, v2) == 0) return 1;//1 重合
else return 0;//0 平行
}
return 2; //2 相交
}
9.两条直线的交点
//求两直线ab和cd的交点
//调用前要保证两直线不平行或重合
//叉积为零则平行或重合
Point Cross_point(Point a, Point b, Point c, Point d) { //Line1:ab, Line2:cd
double s1 = Cross(b - a, c - a);
double s2 = Cross(b - a, d - a); //叉积有正负
return Point(c.x * s2 - d.x * s1, c.y * s2 - d.y * s1) / (s2 - s1);
}
10.判断两个线段是否相交
利用叉积正负来判断;如果一条线段的两个端点在一条线段的两侧,那么这两个端点与另一线段产生的两个叉积正负相反,也就是说两个叉积相乘为负.
//两线段是否相交：1 相交，0不相交
bool Cross_segment(Point a, Point b, Point c, Point d) {//Line1:ab, Line2:cd
double c1 = Cross(b - a, c - a), c2 = Cross(b - a, d - a);
double d1 = Cross(d - c, a - c), d2 = Cross(d - c, b - c);
return sgn(c1) * sgn(c2) < 0 && sgn(d1) * sgn(d2) < 0;//注意交点是端点的情况不算在内
}
11.求两条线段的交点
先判断是否相交,若相交,问题转化为直线交点问题
多边形
三角形和普通多边形定义
1.外心 : 三边中垂线交点，到三角形三个顶点距离相同
2.内心 : 角平分线的交点，到三角形三边的距离相同
3.垂心 : 三条高线的交点
4.重心 : 三条中线的交点，到三角形三顶点距离的平方和最小的点，三角形内到三边距离之积最大的点
5.普通多边形 通常按照逆时针储存所有顶点
- 多边形
由在同一平面且不再同一直线上的多条线段首位顺次连接且不相交所组成的图形交多边形
7.简单多边形
简单多边形是除相邻边外其它边不相交的多边形
8.凸多边形
过多边形的任意一边做一条直线，如果其他各个顶点都在这条直线的同侧，则把这个多边形叫做凸多边形
任意凸多边形外角和均为360°
任意凸多边形内角和为(n−2)180°
判断点在多边形内部
1.转角法思想的一种实现:
以点p为起点引一条水平线,检查与多边形每条边的相交情况
检查以下3个参数:c = cross(p - j,i - j); u = i.y - p.y; v = j.y - p.y
叉积c用来检查p点在线段ij的左侧还是右侧;u,v用来检查经过p的水平线是否击穿过线段ij;
用num计数
if(c > 0 && u < 0 && v >= 0) num ++;
if(c < 0 && u >= 0 && v < 0) num --;
多边形的形状是由各个顶点的排列顺序决定的
//判断点和任意多边形的关系: 3 点上; 2 边上; 1 内部; 0 外部
int Point_in_polygon(Point pt, Point *p, int n) { //点pt，多边形Point *p
for (int i = 0; i < n; i++) { //点在多边形的顶点上
if (p[i] == pt)return 3;
}
for (int i = 0; i < n; i++) {//点在多边形的边上
Line v = Line(p[i], p[(i + 1) % n]);
if (Point_on_seg(pt, v)) return 2;
}
int num = 0;
for (int i = 0; i < n; i++) {
int j = (i + 1) % n;
int c = sgn(Cross(pt - p[j], p[i] - p[j]));
int u = sgn(p[i].y - pt.y);
int v = sgn(p[j].y - pt.y);
if (c > 0 && u < 0 && v >= 0) num++;
if (c < 0 && u >= 0 && v < 0) num--;
}
return num != 0; //1 内部; 0 外部
}
2.pnpoly算法
假设多边形的坐标存放在一个数组里，首先我们需要取得该数组在横坐标和纵坐标的最大值和最小值，根据这四个点算出一个四边型，首先判断目标坐标点是否在这个四边型之内，如果在这个四边型之外，那可以跳过后面较为复杂的计算，直接返回false。
if (p.x < minX || p.x > maxX || p.y < minY || p.y > maxY) {
return false;
// 这个测试都过不了。。。直接返回false；
}
核心算法:
参数nvert 代表多边形有几个点。浮点数testx, testy代表待测试点的横坐标和纵坐标，*vertx,*verty分别指向储存多边形横纵坐标数组的首地址。
int pnpoly (int nvert, float *vertx, float *verty, float testx, float testy) {
int i, j, c = 0;
for (i = 0, j = nvert-1; i < nvert; j = i++) {
if ( ( (verty[i]>testy) != (verty[j]>testy) ) &&
(testx < (vertx[j]-vertx[i]) * (testy-verty[i]) / (verty[j]-verty[i]) + vertx[i]) )
c = !c;
}
return c;
}
3 .射线法
以该点为起点引一条射线，与多边形的边界相交奇数次，说明在多边形的内部。
int pointin(Point p,Point *a,int n) {
int wn = 0;
for (int i = 0;i < n;i ++) {
if (sgn(on_segment(p, a[i], a[(i+1)%n]))==0) return -1;//判断点是否在多边形的边界上
int k = sgn(Cross(a[(i+1-1)%n+1]-a[i], p-a[i]));
int d1 = sgn(a[i].y-p.y);
int d2 = sgn(a[(i+1-1)%n+1].y-p.y);
if (k>0 && d1<=0 && d2>0) wn++;
if (k<0 && d2<=0 && d1>0) wn--;
}
if (wn) return 1;
else return 0;
}
多边形的面积
以原点为中心划分三角形,然后求多边形的面积
double Polygon_area(Point *p, int n) { //Point *p表示多边形。从原点开始划分三角形
double area = 0;
for (int i = 0; i < n; i++)
area += Cross(p[i], p[(i + 1) % n]);
return area / 2; //面积有正负，不能简单地取绝对值
}
求多边形的重心
hdu 1115
//多边形重心
//将多边形三角剖分，算出每个三角形重心（三角形重心是3点坐标平均值）
//对每个三角形有向面积求加权平均值
Point Polygon_center(Point *p, int n) { //求多边形重心。Point *p表示多边形。
Point ans(0, 0);
if (Polygon_area(p, n) == 0) return ans;
for (int i = 0; i < n; i++)
ans = ans + (p[i] + p[(i + 1) % n]) * Cross(p[i], p[(i + 1) % n]); //面积有正负
return ans / Polygon_area(p, n) / 6.;
}
圆
圆的定义
struct Circle{
Point c;
double r;
Circle(){};
Circle(Point c, double r) : c(c), r(r){};
Circle(double x, double y, double _r) { c = Point(x, y); r = _r;}
};
点和圆的关系根据点到圆心的距离判断
int Point_circle_relation(Point p,Circle C){
double dst = Distance(p, C.c);
if(sgn(dst - C.r) < 0) return 0; //点在园内
if(sgn(dst - C.r) == 0) return 1; //点在圆上
return 2; //点在圆外
}
直线和圆的关系根据圆心到直线的距离判断
int Line_circle_relation(Line v,Circle C){
double dst = Dis_point_line(C.c, v);
if(sgn(dst - C.r) < 0) return 0; //直线与圆相交
if (sgn(dst - C.r) == 0) return 1; //直线与圆相切
return 2; //直线在圆外
}
线段和圆的关系根据圆心到线段的距离判断
int Seg_circle_relation(Segment v,Circle C){
double dst = Dis_point_seg(C.c, v);
if(sgn(dst - C.r) < 0) return 0; //线段在园内
if(sgn(dst - C.r) == 0) return 1; //线段在圆上
return 2; //线段在圆外
}
直线和圆的交点pa,pb是交点,返回值是交点个数
先求圆心c在直线上的投影q.再求距离d,然后根据r和d求出长度k,最后求出两个交点pa = q + n*k,pb = q - n * k;其中n是直线的单位向量
int Line_cross_circle(Line v,Circle C,Point &pa,Point &pb){
if(Line_circle_relation(v,C) == 2) return 0;
Point q = Point_line_proj(C.c, v);
double d = Dis_point_line(C.c, v);
double k = sqrt(C.r * C.r - d * d);
if(sgn(k) == 0){
pa = q;pb = q;return 1;
}
Point n = (v.p2 - v.p1) / vector_length(v.p2 - v.p1);
pa = q + n * k;
pb = q - n * k;
return 2;
}
典例 :hdu 5572
https://blog.csdn.net/gpc_123/article/details/122782827
#include <bits/stdc++.h>
using namespace std;
const double eps = 1e-8;//偏差值
int sgn(double x) {
if (fabs(x) < eps)return 0;
else return x < 0 ? -1 : 1;
}
struct Point {
double x, y;
Point() {}
Point(double x, double y) : x(x), y(y) {}
Point operator+(Point B) { return Point(x + B.x, y + B.y); }
Point operator-(Point B) { return Point(x - B.x, y - B.y); }
Point operator*(double k) { return Point(x * k, y * k); }//放大k倍
Point operator/(double k) { return Point(x / k, y / k); }//缩小k倍
};
typedef Point Vector;
double Dot(Vector A, Vector B) { return A.x * B.x + A.y * B.y;}
double vector_length(Vector A) { return sqrt(Dot(A, A));}
double vector_length_square(Vector A) { return Dot(A, A);}
double Cross(Vector A, Vector B) { return A.x * B.y - A.y * B.x;}
double Distance(Point A, Point B) { return hypot(A.x - B.x, A.y - B.y);}
struct Line {
Point p1, p2;
Line() {}
//根据端点确定直线
Line(Point p1, Point p2) : p1(p1), p2(p2) {}
};
typedef Line Segment;
int Point_line_relation(Point p, Line v) {
int c = sgn(Cross(p - v.p1, v.p2 - v.p1));
if (c < 0)return 1; //1：p在v的左边
if (c > 0)return 2; //2：p在v的右边
return 0; //0：p在v上
}
double Dis_point_line(Point p, Line v) {
return fabs(Cross(p - v.p1, v.p2 - v.p1)/ Distance(v.p1, v.p2));
}
double Dis_point_seg(Point p, Segment v) {
if (sgn(Dot(p - v.p1, v.p2 - v.p1)) < 0 || sgn(Dot(p - v.p2, v.p1 - v.p2)) < 0) //点的投影不在线段上
return min(Distance(p, v.p1), Distance(p, v.p2));
return Dis_point_line(p, v); //点的投影在线段上
}
Point Point_line_proj(Point p, Line v) {
double k = Dot(v.p2 - v.p1, p - v.p1) / vector_length_square(v.p2 - v.p1);
return v.p1 + (v.p2 - v.p1) * k;
}
Point Point_line_symmetry(Point p, Line v) {
Point q = Point_line_proj(p, v);
return Point(2 * q.x - p.x, 2 * q.y - p.y);
}
struct Circle{
Point c;
double r;
Circle(){};
Circle(Point c, double r) : c(c), r(r){};
Circle(double x, double y, double _r) { c = Point(x, y); r = _r;}
};
int Seg_circle_relation(Segment v,Circle C){
double dst = Dis_point_seg(C.c, v);
if(sgn(dst - C.r) < 0) return 0; //线段在园内
if(sgn(dst - C.r) == 0) return 1; //线段在圆上
return 2; //线段在圆外
}
int Line_circle_relation(Line v,Circle C){
double dst = Dis_point_line(C.c, v);
if(sgn(dst - C.r) < 0) return 0; //直线与圆相交
if (sgn(dst - C.r) == 0) return 1; //直线与圆相切
return 2; //直线在圆外
}
int Line_cross_circle(Line v,Circle C,Point &pa,Point &pb){
if(Line_circle_relation(v,C) == 2) return 0;
Point q = Point_line_proj(C.c, v);
double d = Dis_point_line(C.c, v);
double k = sqrt(C.r * C.r - d * d);
if(sgn(k) == 0){
pa = q;pb = q;return 1;
}
Point n = (v.p2 - v.p1) / vector_length(v.p2 - v.p1);
pa = q + n * k;
pb = q - n * k;
return 2;
}
int main(){
int T;scanf("%d",&T);
for (int cas = 1; cas <= T;++ cas){
Circle O; Point A, B, V;
scanf("%lf%lf%lf", &O.c.x, &O.c.y, &O.r);
scanf("%lf%lf%lf%lf", &A.x, &A.y, &V.x, &V.y);
scanf("%lf%lf", &B.x, &B.y);
Line l(A, A + V); //射线
Line t(A, B);
//情况一 :直线和圆不相交,而且直线经过点
if(Point_line_relation(B,l) == 0 && Seg_circle_relation(t,O) >= 1 && sgn(Cross(B - A,V) == 0))
printf("Case #%d : YES\n", cas);
else{
Point pa, pb;
//情况二 :直线和圆相切,不经过点
if(Line_cross_circle(l,O,pa,pb) != 2)
printf("Case #%d : NO\n", cas);
else{
//情况三 :直线和圆相交
Point cut; //直线和圆的碰撞点
if(Distance(pa,A) > Distance(pb,A)) cut = pb;
else cut = pa;
Line mid(cut, O.c); //圆心到碰撞点的直线
Point en = Point_line_symmetry(A, mid); //镜像点
Line light(cut, en); //反射线
if(Distance(light.p2,B) > Distance(light.p1,B))
swap(light.p1, light.p2);
if(sgn(Cross(light.p2 - light.p1,Point(B.x - cut.x,B.y - cut.y))) == 0)
printf("Case #%d : YES\n", cas);
else
printf("Case #%d : NO\n", cas);
}
}
}
return 0;
}
最小圆覆盖
给定n个点的平面坐标,求一个半径最小的圆,把n个点全部包围,部分点在圆上;
常见两种算法 : 几何算法,模拟退火
1 . 几何算法
一个最小圆可以由两个点或者三个点来确定;两个点是圆心就是线段AB的中点,三个点的圆心就是三角形的外心;
增量法 :
-
先加入一个点p1,然后c1的圆心就是p1,半径为0;
-
加入第二个点p2,然后c2的圆心就是线段p1,p2的中心,半径为两点距离的一半;
-
加入第三个点p3;有两种情况 :
p3在c2的内部或者圆上,不影响原来的最小圆;
p3在c2的外部,此时需要更新;因为p3一定需要在c3上,现在的任务就是在p1,p2中找到一个点或者两个点,与p3一起两点定圆或者三点定圆;定圆的过程回到1,然后依次加入p1,p2;
算法复杂度接近O(n);
最小圆覆盖裸题 :hdu3007
输入n个点的坐标,n < 500,求最小圆覆盖
#include <bits/stdc++.h> using namespace std; const double eps = 1e-8;//偏差值 const int maxn = 505;//点的数量 int sgn(double x) { if (fabs(x) < eps)return 0; else return x < 0 ? -1 : 1; } struct Point { double x, y; }; double Distance(Point A, Point B) { return hypot(A.x - B.x, A.y - B.y);} //求三角形abc的外接圆的圆心 Point circle_center(const Point a,const Point b,const Point c){ Point center; double a1 = b.x - a.x, b1 = b.y - a.y, c1 = (a1 * a1 + b1 * b1) / 2; double a2 = c.x - a.x, b2 = c.y - a.y, c2 = (a2 * a2 + b2 * b2) / 2; double d = a1 * b2 - a2 * b1; center.x = a.x + (c1 * b2 - c2 * b1) / d; center.y = a.y + (a1 * c2 - a2 * c1) / d; return center; } //求最小覆盖圆,返回圆心c,半径r void min_cover_circle(Point *p,int n,Point &c,double &r){ random_shuffle(p, p + n); //随机函数,打乱所有点 c = p[0];r = 0; //从第一个点p0开始,圆心为p0,半径为0 for (int i = 1; i < n;i++){ //扩展所有点 if(sgn(Distance(p[i],c) - r) > 0){ //点pi在圆外 c = p[i];r = 0; //重新设置圆心为pi,半径为0 for (int j = 0; j < i;j ++){ //重新检测前面所有点 if(sgn(Distance(p[j],c) - r) > 0){ //两点定圆 c.x = (p[i].x + p[j].x) / 2; c.y = (p[i].y + p[j].y) / 2; r = Distance(p[j], c); for (int k = 0; k < j;k ++){ if(sgn(Distance(p[k],c) - r) > 0){ //两点定圆不行就是三点定圆 c = circle_center(p[i], p[j], p[k]); r = Distance(p[i], c); } } } } } } } int main(){ int n; Point p[maxn]; Point c;double r; while(~scanf("%d", & n) && n){ for (int i = 0; i < n;i ++) scanf("%lf %lf",&p[i].x,&p[i].y); min_cover_circle(p, n, c, r); printf("%.2f %.2f %.2f\n", c.x, c.y, r); } return 0; }2 . 模拟退火算法
用模拟退火求最小圆,不断迭代降温;在每次迭代时,找到能覆盖到所有点的一个圆;在多次迭代中,逐步逼近最后要求的圆心和半径;
//求最小覆盖圆,返回圆心c,半径r void min_cover_circle(Point *p,int n,Point &c,double &r){ double T = 100.0; double delta = 0.98; c = p[0]; int pos; while(T > eps){ pos = 0;r = 0; for (int i = 0; i <= n - 1;i ++){ if(Distance(c,p[i]) > r){ r = Distance(c, p[i]); pos = i; } } c.x += (p[pos].x - c.x) / r * T; c.y += (p[pos].y - c.y) / r * T; T *= delta; } }
199
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
