计算几何全家桶
完整板子
#include<bits/stdc++.h>
using namespace std;
const double eps = 1e-8,Pi = acos(-1.0);
inline int dcmp(double x){return (x<-eps)?-1:(x>eps?1:0);}
inline double Abs(double x){return x * dcmp(x);};
namespace Flat_Space{
struct pt{
double x,y; pt(double x_ = 0,double y_ = 0){x = x_;y = y_;}
inline void read(){scanf("%lf%lf",&x,&y);}
};
inline bool cmpx(const pt &a,const pt &b){return (a.x != b.x) ? a.x < b.x : a.y < b.y;}
typedef pt vec;
inline double Len(const vec &a){return sqrt(a.x*a.x + a.y*a.y);}//模长
inline double angle(const vec &a){return atan2(a.y,a.x);}//象限角
inline vec operator + (const vec &a,const vec &b){return vec(a.x+b.x,a.y+b.y);}//加
inline vec operator - (const vec &a,const vec &b){return vec(a.x-b.x,a.y-b.y);}//减
inline double operator * (const vec &a,const vec &b){return a.x*b.x+a.y*b.y;}//点积
inline double operator ^ (const vec &a,const vec &b){return a.x*b.y-a.y*b.x;}//叉积
inline vec operator * (const vec &a,const double b){return vec(a.x*b,a.y*b);}//数乘
inline vec operator / (const vec &a,const double b){return vec(a.x/b,a.y/b);}//数除
inline vec rotate(const vec &a,const double theta){return vec(a.x*cos(theta)-a.y*sin(theta),a.x*sin(theta)+a.y*cos(theta));}//逆时针转theta
/* [ cos , sin ]
[ -sin , cos ] */
inline vec rotate_90(const vec &a){return vec(-a.y,a.x);}//逆时针90
inline vec rotate_P(const vec &a,const vec &b,const double theta){return rotate(a-b,theta)+b;}//a绕b逆时针转theta
struct Line{//Point,Segment,Ray,Line
pt s,t;
Line(pt s_=pt(0,0),pt t_=pt(0,0)){s = s_;t = t_;}
};
inline double maxx(const Line &L){return max(L.s.x,L.t.x);}
inline double minx(const Line &L){return min(L.s.x,L.t.x);}
inline double maxy(const Line &L){return max(L.s.y,L.t.y);}
inline double miny(const Line &L){return min(L.s.y,L.t.y);}
inline double angle(const Line &L){return angle(L.t-L.s);}
inline bool operator < (const Line &a,const Line &b){
double a1 = angle(a),a2 = angle(b);
if(dcmp(a2-a1) != 0) return dcmp(a2-a1) > 0;
return dcmp((b.t-a.s)^(a.t-b.s)) > 0;
}
inline bool operator == (const pt &a,const pt &b){return (!dcmp(a.x-b.x)) && (!dcmp(a.y-b.y));}
inline bool judge_PL(const pt &a,const Line &b){return !dcmp((a-b.s)^(b.t-b.s));}
inline bool judge_PS(const pt &a,const Line &b){return ((!dcmp((a-b.s)^(b.t-b.s)) && (dcmp((a-b.s)*(a-b.t)) <= 0)));}
inline double dis_PP(const pt &a,const pt &b){return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));}
inline double dis_PL(const pt &a,const Line &b){return Abs((b.t-b.s) ^ (a-b.s)) / Len(b.t-b.s);}
inline double dis_PS(const pt &a,const Line &b){
pt x = a-b.s,y = a-b.t,z = b.t-b.s;
if(dcmp(x*z) < 0) return Len(x);
if(dcmp(y*z) > 0) return Len(y);
return dis_PL(a,b);
}
inline pt point_PL(const pt &a,const Line &b){
pt x = a-b.s,y = a-b.t,z = b.t-b.s;
double p1 = x * z,p2 = -(y*z);
return b.s + z * (p1/(p1+p2));
}
inline pt point_PS(const pt &a,const Line &b){
pt x = a-b.s,y = a-b.t,z = b.t-b.s;
if(dcmp(x*z) < 0) return x;
if(dcmp(y*z) > 0) return y;
return point_PL(a,b);
}
inline pt cross_LL(const Line &a,const Line &b){
pt x = a.t-a.s,y = b.t-b.s,z = a.s-b.s;
return a.s + x * (y^z)/(x^y);
}
inline bool judge_cross_SL(const Line &a,const Line &b){
if(!dcmp((a.t-a.s) ^ (b.t-b.s))) return false;
pt z = b.t - b.s;
return (dcmp(a.s ^ z) * dcmp(a.t ^ z)) <= 0;
}
inline bool judge_cross_SS(const Line &a,const Line &b){
if(maxx(a) < minx(b) || maxx(b) < minx(a) || maxy(a) < miny(b) || maxy(b) < miny(a)) return false;
pt x = a.t-a.s,y = b.t-b.s;
return (dcmp(a.s ^ y) * dcmp(a.t ^ y) <= 0) && (dcmp(b.s ^ x) * dcmp(b.t ^ x) <= 0);
}
inline pt mirror(pt a,Line b){return point_PL(a,b) * 2.0 - a;}
}
using namespace Flat_Space;
namespace Polygon{//多边形
struct polygon{
vector<pt> pts;//默认逆时针
inline void insert(pt p){pts.push_back(p);}
inline pt& operator [](int x){return pts[x];}
inline void clear(){pts.clear();}
inline int nxt(int x){return x == pts.size()-1 ? 0 : x+1;}
inline double peri(){//perimeter
double ans = 0;
for(int i = 0; i < pts.size(); ++i) ans += dis_PP(pts[i],pts[nxt(i)]);
return ans;
}
inline double area(){
double ans = 0;
for(int i = 1; i < pts.size()-1; ++i) ans += (pts[i]-pts[0]) ^ (pts[nxt(i)]-pts[0]);
return ans/2.0;
}
inline bool Left(pt x,pt l,pt r){return dcmp((l-x)^(r-x)) <= 0;}// xl 在 xr 的方向
inline void rever(){reverse(pts.begin(),pts.end());}
inline int include(pt p){//out:0 in:1 on:2 做p点平行于x轴直线，看是不是与多边形有奇数个交点
int cnt = 0;
for(int i = 0; i < pts.size(); ++i){
pt s = pts[i],t = pts[nxt(i)];
Line L = Line(s,t);
if(judge_PS(p,L)) return 2;
if(dcmp(p.y-miny(L)) >= 0 && dcmp(p.y-maxy(L)) < 0){
double cx = s.x + ((p.y-s.y)/(t.y-s.y) * (t.x-s.x));
if(dcmp(cx-p.x) > 0) ++cnt;
}
}
return cnt & 1;
}
inline void include_bs(pt p){
int n = pts.size();
if(!Left(pts[0],p,pts[1])) return 0;
if(!Left(pts[n-1],p,pts[0])) return 0;
if(judge_PS(p,Line(pts[0],pts[1]))) return 2;
if(judge_PS(p,Line(pts[0],pts[n-1]))) return 2;
int l = 1, r = n-2,ans = 1;
while(l <= r){
int mid = (l + r) / 2;
if(!Left(pts[0],pts[mid],p)) l = mid + 1,ans = mid;
else r = mid - 1;
}
if(!Left(pts[ans],p,pts[nxt(ans)])) return 0;
if(judge_PS(p,Line(pts[ans],pts[nxt(ans)]))) return 2;
return 1;
}
};
inline bool disjoint(polygon A,polygon B){
for(int i = 0; i < A.pts.size(); ++i){
Line x = Line(A[i],A[A.nxt(i)]);
for(int j = 0; j < B.pts.size(); ++j){
Line y = Line(B[j],B[B.nxt(j)]);
if(judge_cross_SS(x,y)) return false;
}
}
if(A.include_bs(B[0])) return false;
if(B.include_bs(A[0])) return false;
return true;
}
}
using namespace Polygon;
namespace Circle{
struct circle{
pt o; double r;
circle(pt o_=pt(0,0),double r_ = 0){o = o_;r = r_;}
inline bool include(pt p){return dcmp(r-dis_PP(o,p)) >= 0;}
inline void print(int x){
printf("%.*lf\n",x,r);
printf("%.*lf %.*lf",x,o.x,x,o.y);
}
circle circum(pt A,pt B,pt C){//外心 外接圆
pt D = (A+B) / 2,E = (B+C) / 2;
o = cross_LL(Line(D,D+(rotate_90(B-A))),Line(E,E+(rotate_90(C-B))));
r = dis_PP(o,A);
}
};
inline circle mincov(const vector<pt> &p){//最小圆覆盖
int n = p.size();
circle c = circle(0,0);
for(int i = 0; i < n; ++i)
if(!c.include(p[i])){
c = circle(p[i],0);
for(int j = 0; j < i; ++j){
if(!c.include(p[j])){
c = circle((p[i]+p[j]) / 2.0,dis_PP(p[i],p[j]) / 2.0);
for(int k = 0; k < j; ++k)
if(!c.include(p[k])) c.circum(p[i],p[j],p[k]);
}
}
}
return c;
}
}
using namespace Circle;
namespace Hull{// 凸包、旋转卡壳、半平面交、闵可夫斯基和
inline void Andrew(polygon &p){
vector<pt> q(p.pts.size() << 1);
sort(p.pts.begin(),p.pts.end(),cmpx);
int top = -1,n = p.pts.size();
q[++top] = p[0];
for(int i = 1; i < n; ++i){
while(top && dcmp((q[top-1]-q[top]) ^ (p[i]-q[top])) >= 0) --top;
q[++top]=p[i];
}
int now = top;
for(int i = n-2; i >= 0; --i){
while(top > now && dcmp((q[top-1]-q[top]) ^ (p[i]-q[top])) >= 0) -- top;
q[++top] = p[i];
}
if(n > 1) --top;
p.clear();
for(int i = 0; i <= top; ++i) p.insert((q[i]));
}
inline double S(const pt &x,const pt &y,const pt &z){return Abs((y-x)^(z-x));}
inline double diam(polygon &p){//diameter
int n = p.pts.size();double ans = 0;
for(int i = 0,j = 1; i < n; ++i){
for(;;j = p.nxt(j)) if(dcmp(S(p[j],p[i],p[p.nxt(i)]))-S(p[p.nxt(j)],p[i],p[p.nxt(i)]) >= 0) break;
ans = max(ans,max(dis_PP(p[j],p[i]),dis_PP(p[j],p[p.nxt(i)])));
}
return ans;
}
inline polygon SI(vector<Line> q){//左半平面交
int n = q.size();
sort(q.begin(),q.end());
vector<Line> li(n+1),L(n+1);
vector<pt> p(n+1);
int len = 0;
for(int i = 0; i < n; ++i){
if(i > 0 && !dcmp(angle(q[i])-angle(q[i-1]))) continue;
L[++len]=q[i];
}
int l = 1,r = 0;
for(int i = 1; i <= len; ++i){
while(l < r && dcmp((L[i].t - p[r]) ^ (L[i].s - p[r])) > 0) --r;
while(l < r && dcmp((L[i].t - p[l+1]) ^ (L[i].s - p[l+1])) > 0) ++l;
li[++r] = L[i];
if(r > l) p[r] = cross_LL(li[r],li[r-1]);
}
while(l < r && dcmp((L[l].t - p[r]) ^ (L[l].s - p[r])) > 0) --r;
while(l < r && dcmp((L[r].t - p[l+1]) ^ (L[r].s - p[l+1])) > 0) ++l;
p[r+1] = cross_LL(li[r],li[l]);++r;
polygon P;
for(int j = l + 1; j <= r; ++j) P.insert(p[j]);
return P;
}
inline polygon merge(polygon A,polygon B){//闵可夫斯基和
int n = A.pts.size(),m = B.pts.size(),tot1 = 0,tot2 = 0;
vector<pt> p(n+1),q(m+1);
for(int i = 1; i < n; ++i) p[++tot1]=A[i]-A[i-1];p[++tot1]=A[0]-A[n-1];
for(int i = 1; i < m; ++i) q[++tot2]=B[i]-B[i-1];q[++tot2]=B[0]-B[m-1];
int i = 1, j = 1, tot = 0;pt last = pt(0,0);
polygon C;
C.pts.resize(n+m+1); C[0] = last = A[0] + B[0];
while(i <= n && j <= m){
if(dcmp(p[i]^q[j]) >= 0) C[++tot] = last + p[i++], last = C[tot];
else C[++tot] = last+q[j++],last = C[tot];
}
while(i <= n) C[++tot] = last + p[i++],last = C[tot];
while(j <= m) C[++tot] = last + q[j++],last = C[tot];
Andrew(C);
return C;
}
}
using namespace Hull;
using namespace Circle;
int main(){
pt x;
}
一、前置芝士
- 高一平面向量芝士
二、精度相关
eps精度相关的常量（个人通常设为1e-8），主要是为处理C++的精度丢失导致的问题Pi即圆周率 \(\pi\) ，可以用acos(-1.0)求出（因为 \(\cos(\pi) = -1\) 啦inline int dcmp(double x)表示 \(x\) 的符号，正 \(1\) 负 \(-1\) 零 \(0\)inline double Abs(double x)和math.h库中的abs()函数作用相同，但是该函数处理了精度问题
code：
const double eps = 1e-8,Pi = acos(-1.0);
inline int dcmp(double x){return (x<-eps)?-1:(x>eps?1:0);}
inline double Abs(double x){return x * dcmp(x);};
三、平面相关
1. 点与向量
① 定义
struct pt/vec点/向量inline bool cmpx(const pt &a,const pt &b)按横坐标排序
struct pt{
double x,y; pt(double x_ = 0,double y_ = 0){x = x_;y = y_;}
inline void read(){scanf("%lf%lf",&x,&y);}
};
inline bool cmpx(const pt &a,const pt &b){return (a.x != b.x) ? a.x < b.x : a.y < b.y;}
typedef pt vec;
② 一元运算
inline double Len(const vec &a)求向量模长inline double angle(const vec &a)求象限角，使用了atan2()函数，可以避免除数为零等麻烦
inline double Len(const vec &a){return sqrt(a.x*a.x + a.y*a.y);}//模长
inline double angle(const vec &a){return atan2(a.y,a.x);}//象限角
③ 二元运算
operator +-*^...重定义了运算符，让代码写起来更加简洁。（为了让后面的内容更加易懂，我们讲一下点积和叉积）- 点积，在几何上通常表示一个向量在另外一个向量上的投影的长度，正负号可以表示两个向量之间的夹角与 \(90°\) 的大小（正小负大）
- 叉积，在几何上通常表示一个向量和另一个向量形成的平行四边形的面积大小，正负号表示两个向量间的位置关系（后向量对于前向量来说 正逆反顺）
inline vec rotate(const vec &a,const double theta)将向量 \(a\) 逆时针旋转 \(\theta\) / 将点 \(a\) 绕原点旋转 \(\theta\)inline vec rotate_90(const vec &a)将向量 \(a\) 逆时针旋转 \(90°\) / 将点 \(a\) 绕原点旋转 \(90°\)inline vec rotate_P(const vec &a,const vec &b,const double theta)将点 \(a\) 绕点 \(b\) 旋转 \(\theta\)
inline vec operator + (const vec &a,const vec &b){return vec(a.x+b.x,a.y+b.y);}//加
inline vec operator - (const vec &a,const vec &b){return vec(a.x-b.x,a.y-b.y);}//减
inline double operator * (const vec &a,const vec &b){return a.x*b.x+a.y*b.y;}//点积
inline double operator ^ (const vec &a,const vec &b){return a.x*b.y-a.y*b.x;}//叉积
inline vec operator * (const vec &a,const double b){return vec(a.x*b,a.y*b);}//数乘
inline vec operator / (const vec &a,const double b){return vec(a.x/b,a.y/b);}//数除
inline vec rotate(const vec &a,const double theta){return vec(a.x*cos(theta)-a.y*sin(theta),a.x*sin(theta)+a.y*cos(theta));}//逆时针转theta
/* [ cos , sin ]
[ -sin , cos ] */
inline vec rotate_90(const vec &a){return vec(-a.y,a.x);}//逆时针90
inline vec rotate_P(const vec &a,const vec &b,const double theta){return rotate(a-b,theta)+b;}//a绕b逆时针转theta
2. 线
① 定义
struct Line线（线段 \(S\)，射线 \(R\)，直线 \(L\)）- \(s,t\) 分别为这条线上的两个点
Line(pt s_=pt(0,0),pt t_=pt(0,0)Line类的初始化
struct Line{//Point,Segment,Ray,Line
pt s,t;
Line(pt s_=pt(0,0),pt t_=pt(0,0)){s = s_;t = t_;}
};
② 一元运算
max(min) x/y横纵坐标的最大/最小值inline double angle(const Line &L)线 \(L\) 的角度
inline double maxx(const Line &L){return max(L.s.x,L.t.x);}
inline double minx(const Line &L){return min(L.s.x,L.t.x);}
inline double maxy(const Line &L){return max(L.s.y,L.t.y);}
inline double miny(const Line &L){return min(L.s.y,L.t.y);}
inline double angle(const Line &L){return angle(L.t-L.s);}
③ 二元运算
inline bool operator < (const Line &a,const Line &b)重载运算符，通过极角进行比较inline bool operator == (const pt &a,const pt &b)判断两条线段是否重合
inline bool operator < (const Line &a,const Line &b){
double a1 = angle(a),a2 = angle(b);
if(dcmp(a2-a1) != 0) return dcmp(a2-a1) > 0;
return dcmp((b.t-a.s)^(a.t-b.s)) > 0;
}
inline bool operator == (const pt &a,const pt &b){return (!dcmp(a.x-b.x)) && (!dcmp(a.y-b.y));}
3. 点线间的复杂运算
① 判定点是否在线上
inline bool judge_PL(const pt &a,const Line &b)求一个点是否在一条直线上- 向量叉积的正负号的几何含义：向量叉积为 \(0\) 说明两个向量方向相同
inline bool judge_PS(const pt &a,const Line &b)求一个点是否在一条线段上- 判断点是否在直线上
- 向量点积的正负号的几何含义：向量点积为负说明两个向量的夹角大于 \(90°\)
inline bool judge_PL(const pt &a,const Line &b){return !dcmp((a-b.s)^(b.t-b.s));}
inline bool judge_PS(const pt &a,const Line &b){return ((!dcmp((a-b.s)^(b.t-b.s)) && (dcmp((a-b.s)*(a-b.t)) <= 0)));}
② 计算距离
inline double dis_PP(const pt &a,const pt &b)求两点之间的距离- 直接勾股定理即可
inline double dis_PL(const pt &a,const Line &b)求点到直线的距离- 我们首先使用叉积求出该点与直线上两点构成的面积，再除以直线上两点的距离即可求出高，也就是该点到直线的距离
inline double dis_PS(const pt &a,const Line &b)求点到线段的距离- 我们考虑先用点积正负号的性质 判断 点与线段的垂足落在线段上/外
- 若落在线段外，即为两点之间的距离
- 若落在线段上，即为点到直线的距离
inline double dis_PP(const pt &a,const pt &b){return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));}
inline double dis_PL(const pt &a,const Line &b){return Abs((b.t-b.s) ^ (a-b.s)) / Len(b.t-b.s);}
inline double dis_PS(const pt &a,const Line &b){
pt x = a-b.s,y = a-b.t,z = b.t-b.s;
if(dcmp(x*z) < 0) return Len(x);
if(dcmp(y*z) > 0) return Len(y);
return dis_PL(a,b);
}
③ 计算垂足
inline pt point_PL(const pt &a,const Line &b)点到直线距离最近的点（垂足）- 我们通过求出点 \(a\) 与直线上两点 \(s,t\) 的连线在直线上的投影 \(p1,p2\)
- 再运用得出的 \(p1,p2\) 的比例求出对应的垂足的坐标
inline pt point_PS(const pt &a,const Line &b)点到线段距离最近的点- 我们考虑还是用点积正负号的性质 判断 点与线段的垂足落在线段上/外
- 若落在线段外，即为线段端点
- 若落在线段上，即为点到直线的垂足
inline pt point_PL(const pt &a,const Line &b){
pt x = a-b.s,y = a-b.t,z = b.t-b.s;
double p1 = x * z,p2 = -(y*z);
return b.s + z * (p1/(p1+p2));
}
inline pt point_PS(const pt &a,const Line &b){
pt x = a-b.s,y = a-b.t,z = b.t-b.s;
if(dcmp(x*z) < 0) return x;
if(dcmp(y*z) > 0) return y;
return point_PL(a,b);
}
④ 计算交点
inline pt cross_LL(const Line &a,const Line &b)求两直线的交点- 通过将面积比转化为线段之比，再使用线段之比通过向量求出交点（后半部分类似垂足）
inline bool judge_cross_SL(const Line &a,const Line &b)求直线与线段是否相交- 首先不能平行
- 线段的两个端点要分别在直线的两侧（这样才能穿过直线，和直线有交）
inline bool judge_cross_SS(const Line &a,const Line &b)求两个线段之间是否相交- 判断两条直线的横/纵坐标覆盖的区间有交
- 判断 \(a\) 的端点在 \(b\) 的两侧并且 \(b\) 的端点在 \(a\) 的两侧
inline pt cross_LL(const Line &a,const Line &b){
pt x = a.t-a.s,y = b.t-b.s,z = a.s-b.s;
return a.s + x * (y^z)/(x^y);
}
inline bool judge_cross_SL(const Line &a,const Line &b){
if(!dcmp((a.t-a.s) ^ (b.t-b.s))) return false;
pt z = b.t - b.s;
return (dcmp(a.s ^ z) * dcmp(a.t ^ z)) <= 0;
}
inline bool judge_cross_SS(const Line &a,const Line &b){
if(maxx(a) < minx(b) || maxx(b) < minx(a) || maxy(a) < miny(b) || maxy(b) < miny(a)) return false;
pt x = a.t-a.s,y = b.t-b.s;
return (dcmp(a.s ^ y) * dcmp(a.t ^ y) <= 0) && (dcmp(b.s ^ x) * dcmp(b.t ^ x) <= 0);
}
⑤ 翻转
inline pt mirror(pt a,Line b)求点 \(a\) 关于直线 \(b\) 的对称点
inline pt mirror(pt a,Line b){return point_PL(a,b) * 2.0 - a;}
四、多边形
1. 数组构建
vector<pt> pts;逆时针排列的点集表示多边形inline void insert(pt p)插入一个点 \(p\)inline pt& operator [](int x)快速访问struct中数组（简化代码）inline void clear()清空多边形inline int nxt(int x)因为多边形是封闭的，构成一个环，所以有该函数计算下一个点的坐标
2. 操作
inline double peri()求多边形周长inline double area()求多边形面积- 就是将多边形分成三角形，求三角形面积
inline bool Left(pt x,pt l,pt r)判断 \(xl\) 是否在 \(xr\) 的逆时针方向inline void rever()将数组中记录多边形的顺序调换(顺时针/逆时针)inline int include(pt p)判断一个点与多边形的位置关系 (外：0、内：1、上：2)- 该函数实现的是从点 \(p\) 向外延伸一条射线，如果和多边形只有一个交点，那么在内部（需要特判在多边形上）。
inline include_bs(pt p)含义同上，实现方式不同- 该函数实现的是判断 \(p\) 是否在每一条线的右边
inline bool disjoint(polygon A,polygon B)判断两个多边形在面积上是否互不重合- 首先判断多边形的边之间是否有交
- 然后再判断两个多边形是否互相包含（判断一个点即可）
struct polygon{
vector<pt> pts;//默认逆时针
inline void insert(pt p){pts.push_back(p);}
inline pt& operator [](int x){return pts[x];}
inline void clear(){pts.clear();}
inline int nxt(int x){return x == pts.size()-1 ? 0 : x+1;}
inline double peri(){//perimeter
double ans = 0;
for(int i = 0; i < pts.size(); ++i) ans += dis_PP(pts[i],pts[nxt(i)]);
return ans;
}
inline double area(){
double ans = 0;
for(int i = 1; i < pts.size()-1; ++i) ans += (pts[i]-pts[0]) ^ (pts[nxt(i)]-pts[0]);
return ans/2.0;
}
inline bool Left(pt x,pt l,pt r){return dcmp((l-x)^(r-x)) <= 0;}// xl 在 xr 的方向
inline void rever(){reverse(pts.begin(),pts.end());}
inline int include(pt p){//out:0 in:1 on:2 做p点平行于x轴直线，看是不是与多边形有奇数个交点
int cnt = 0;
for(int i = 0; i < pts.size(); ++i){
pt s = pts[i],t = pts[nxt(i)];
Line L = Line(s,t);
if(judge_PS(p,L)) return 2;
if(dcmp(p.y-miny(L)) >= 0 && dcmp(p.y-maxy(L)) < 0){
double cx = s.x + ((p.y-s.y)/(t.y-s.y) * (t.x-s.x));
if(dcmp(cx-p.x) > 0) ++cnt;
}
}
return cnt & 1;
}
inline include_bs(pt p){
int n = pts.size();
if(!Left(pts[0],p,pts[1])) return 0;
if(!Left(pts[n-1],p,pts[0])) return 0;
if(judge_PS(p,Line(pts[0],pts[1]))) return 2;
if(judge_PS(p,Line(pts[0],pts[n-1]))) return 2;
int l = 1, r = n-2,ans = 1;
while(l <= r){
int mid = (l + r) / 2;
if(!Left(pts[0],pts[mid],p)) l = mid + 1,ans = mid;
else r = mid - 1;
}
if(!Left(pts[ans],p,pts[nxt(ans)])) return 0;
if(judge_PS(p,Line(pts[ans],pts[nxt(ans)]))) return 2;
return 1;
}
};
inline bool disjoint(polygon A,polygon B){
for(int i = 0; i < A.pts.size(); ++i){
Line x = Line(A[i],A[A.nxt(i)]);
for(int j = 0; j < B.pts.size(); ++j){
Line y = Line(B[j],B[B.nxt(j)]);
if(judge_cross_SS(x,y)) return false;
}
}
if(A.include_bs(B[0])) return false;
if(B.include_bs(A[0])) return false;
return true;
}
}
五、圆
1. 定义
pt o; double r;圆的圆心和半径circle(pt o_=pt(0,0),double r_ = 0){o = o_;r = r_;}圆的初始化inline bool include(pt p){return dcmp(r-dis_PP(o,p)) >= 0;}判断一个点是否在圆内inline void print(int x)输出该结构体的信息，保留 \(x\) 位小数circle circum(pt A,pt B,pt C)三角形外心，外接圆- 取两条中垂线交点为圆心，圆心到任意一个点的距离即为半径
struct circle{
pt o; double r;
circle(pt o_=pt(0,0),double r_ = 0){o = o_;r = r_;}
inline bool include(pt p){return dcmp(r-dis_PP(o,p)) >= 0;}
inline void print(int x){
printf("%.*lf\n",x,r);
printf("%.*lf %.*lf",x,o.x,x,o.y);
}
circle circum(pt A,pt B,pt C){//外心 外接圆
pt D = (A+B) / 2,E = (B+C) / 2;
o = cross_LL(Line(D,D+(rotate_90(B-A))),Line(E,E+(rotate_90(C-B))));
r = dis_PP(o,A);
}
};
3. 最小圆覆盖
- 最小圆覆盖，顾名思义就是用一个最小的圆覆盖所有点。
- 最小覆盖圆，顾名思义就是覆盖所有点的最小的圆。
有一个定理：如果 \(p\) 不在集合 \(S\) 的最小圆覆盖内，那么它一定在集合 \({p}\cup\)S 的最小覆盖圆上
根据这个定理我们就可以快速求出 \(n\) 个点的最小圆覆盖，求法如下：
- 我们先求出前 \(i-1\) 个点的最小圆覆盖 \(C\)
- 如果第 \(i\) 个点在最小圆覆盖内，那么前 \(i\) 个点的最小圆覆盖就是 \(C\)
- 如果不在，那么 \(i\) 一定在最小圆覆盖上，接着我们需要找到另外两个在最小圆覆盖上的点。
- 然后我们每次找不在当前圆中的点 \(j\) 然后将点对 \((i,j)\) 为直径的点当做当前的最小圆
- 接着我们找找不在当前圆中的第三个点 \(k\) 再把点 \(k\) 加进去。
最后记得要 shuffle 哦！
你可能认为它是 \(n^3\) 的，但是实际上它期望是 \(O(n)\) 的。
证明
由于一堆点的最小圆覆盖中最多只有 \(3\) 个点参与最小圆覆盖的构建，所以一个点在调用下一层的概率是 \(\frac 3 n\) 的
我们可以写出计算复杂度的式子：
\[\begin{aligned}
T_1(n) = & O(n) + \sum_{i=1}^n \frac 3 n T_2(i)\\
T_2(n) = & O(n) + \sum_{i=1}^n \frac 3 n T_3(i)\\
T_3(n) = O(n)
\end{aligned}
\]
可以推出 最后的复杂度是 \(O(n)\) 的
inline circle mincov(const vector<pt> &p){//最小圆覆盖（你需要在外面对 p 数组进行shuffle 操作）
int n = p.size();
circle c = circle(0,0);
for(int i = 0; i < n; ++i)
if(!c.include(p[i])){
c = circle(p[i],0);
for(int j = 0; j < i; ++j){
if(!c.include(p[j])){
c = circle((p[i]+p[j]) / 2.0,dis_PP(p[i],p[j]) / 2.0);
for(int k = 0; k < j; ++k)
if(!c.include(p[k])) c.circum(p[i],p[j],p[k]);
}
}
}
return c;
}
六、凸包
1. 凸包
我们考虑大家应该都会斜率优化中的凸包吧，求凸包有几种方法：
- 一个是将凸包分成上下两半分别求
- 另一个是按极角排序之后求
下面代码中就是先求上班部分再求下半部分
inline void Andrew(polygon &p){
vector<pt> q(p.pts.size() << 1);
sort(p.pts.begin(),p.pts.end(),cmpx);
int top = -1,n = p.pts.size();
q[++top] = p[0];
for(int i = 1; i < n; ++i){
while(top && dcmp((q[top-1]-q[top]) ^ (p[i]-q[top])) >= 0) --top;
q[++top]=p[i];
}
int now = top;
for(int i = n-2; i >= 0; --i){
while(top > now && dcmp((q[top-1]-q[top]) ^ (p[i]-q[top])) >= 0) -- top;
q[++top] = p[i];
}
if(n > 1) --top;
p.clear();
for(int i = 0; i <= top; ++i) p.insert((q[i]));
}
2. 旋转卡壳
旋转卡壳是一种求多边形直径/最小矩形覆盖的方法。
具体方式就是我们旋转两条边，并且将多边形夹在中间即可。
inline double S(const pt &x,const pt &y,const pt &z){return Abs((y-x)^(z-x));}
inline double diam(polygon &p){//diameter
int n = p.pts.size();double ans = 0;
for(int i = 0,j = 1; i < n; ++i){
for(;;j = p.nxt(j)) if(dcmp(S(p[j],p[i],p[p.nxt(i)]))-S(p[p.nxt(j)],p[i],p[p.nxt(i)]) >= 0) break;
ans = max(ans,max(dis_PP(p[j],p[i]),dis_PP(p[j],p[p.nxt(i)])));
}
return ans;
}
3. 半平面交
我们首先将所有边按极角排序
然后我们需要维护一个双端队列记录所有的交点，然后我们像维护斜率优化一样维护凸包即可，只不过需要在队列的两端都要 pop
inline polygon SI(vector<Line> q){//左半平面交
int n = q.size();
sort(q.begin(),q.end());
vector<Line> li(n+1),L(n+1);
vector<pt> p(n+1);
int len = 0;
for(int i = 0; i < n; ++i){
if(i > 0 && !dcmp(angle(q[i])-angle(q[i-1]))) continue;
L[++len]=q[i];
}
int l = 1,r = 0;
for(int i = 1; i <= len; ++i){
while(l < r && dcmp((L[i].t - p[r]) ^ (L[i].s - p[r])) > 0) --r;
while(l < r && dcmp((L[i].t - p[l+1]) ^ (L[i].s - p[l+1])) > 0) ++l;
li[++r] = L[i];
if(r > l) p[r] = cross_LL(li[r],li[r-1]);
}
while(l < r && dcmp((L[l].t - p[r]) ^ (L[l].s - p[r])) > 0) --r;
while(l < r && dcmp((L[r].t - p[l+1]) ^ (L[r].s - p[l+1])) > 0) ++l;
p[r+1] = cross_LL(li[r],li[l]);++r;
polygon P;
for(int j = l + 1; j <= r; ++j) P.insert(p[j]);
return P;
}
4. 闵可夫斯基和
闵可夫斯基和就是将一个多边形 \(B\) 挂在多边形 \(A\) 的每一个定点后构成的凸包
具体操作方式就是我们将所有边长代表的向量求出来，按极角排序。
维护两个指针 \(i,j\) 表示我们分别将多边形 \(A,B\) 的前 \(i/j\) 条边加入了新的凸包中，然后我们每次加入极角更小的边（向量）即可。
inline polygon merge(polygon A,polygon B){//闵可夫斯基和
int n = A.pts.size(),m = B.pts.size(),tot1 = 0,tot2 = 0;
vector<pt> p(n+1),q(m+1);
for(int i = 1; i < n; ++i) p[++tot1]=A[i]-A[i-1];p[++tot1]=A[0]-A[n-1];
for(int i = 1; i < m; ++i) q[++tot2]=B[i]-B[i-1];q[++tot2]=B[0]-B[m-1];
int i = 1, j = 1, tot = 0;pt last = pt(0,0);
polygon C;
C.pts.resize(n+m+1); C[0] = last = A[0] + B[0];
while(i <= n && j <= m){
if(dcmp(p[i]^q[j]) >= 0) C[++tot] = last + p[i++], last = C[tot];
else C[++tot] = last+q[j++],last = C[tot];
}
while(i <= n) C[++tot] = last + p[i++],last = C[tot];
while(j <= m) C[++tot] = last + q[j++],last = C[tot];
Andrew(C);
return C;
}
本文来自博客园，作者：ricky_lin，转载请注明原文链接：https://www.cnblogs.com/rickylin/p/17993311/summary_computation-geometry
