计算几何【模板】
TLE_Automat · · 个人记录
此博客是为了存放模板开的blog。
#include <bits/stdc++.h>
using namespace std;
#define print_v(a) printf("%.2lf %.2lf\n",a.x,a.y)
#define print_l(x) putchar('\n'),print_v(x.a),print_v(x.b)
#define SZ(x) ((int)((x).size()))
typedef long long ll;
typedef long double ldb;
typedef pair<int,int> pii;
const int MAXN=2e5+10;
const double eps=1e-9;
const double pi=acos(-1);
inline bool dcmp(double x,double y=0){ return fabs(x-y)<=eps; }
struct Vec // 向量
{
double x,y;
Vec(double xx=0,double yy=0){ x=xx; y=yy; }
// 向量四则运算
inline Vec operator + (const Vec &p) const { return Vec(x+p.x,y+p.y); }
inline Vec operator - (const Vec &p) const { return Vec(x-p.x,y-p.y); }
inline Vec operator * (const double &d) const { return Vec(x*d,y*d); }
inline Vec operator / (const double &d) const { return Vec(x/d,y/d); }
};
typedef Vec Pt; // 把点看成从原点出发，终点为该点的向量
inline double dot(const Vec &a,const Vec &b){ return a.x*b.x+a.y*b.y; } // 向量点乘
inline double cross(const Vec &a,const Vec &b){ return a.x*b.y-a.y*b.x; } // 向量叉乘
inline double len(const Vec &a){ return sqrt(dot(a,a)); } // 向量模长
inline double len_sq(const Vec &a){ return dot(a,a); } // 向量模长的平方
inline double dis(const Pt &a,const Pt &b){ return len(b-a); } // 两点距离
// 将向量 a 逆时针旋转弧度 delta
inline Vec rolate(const Vec &a,const double &rad)
{
double t=atan2(a.y,a.x)+rad,l=len(a);
return Vec(l*cos(t),l*sin(t));
}
struct Line // 直线
{
Pt a,b; //用两点表示一条直线
double ang; //直线与 x 轴的夹角（极角）
// 注意这里的直线是有向直线 , 规定表示方向的向量是 b-a
// 因为我们要分出这条直线的左侧和右侧两个部分 , 所以直线必须有向
Line(){}
Line(Pt aa,Pt bb){ a=aa; b=bb; ang=atan2((b-a).y,(b-a).x); }
inline bool include(const Pt &p){ return dcmp(cross(a-p,b-p)); } // 判断一点是否被直线包含
};
struct Seg // 基本跟直线一模一样 , 有时用的时候方便 , 所以就写了
{
Pt a,b; // 两个端点
Seg(){}
Seg(Pt aa,Pt bb){ a=aa; b=bb; }
inline bool include(const Pt &p){ return dcmp(cross(a-p,b-p)) && dot(a-p,b-p)<=0; }
};
/*
* 判断两直线的位置关系
* 相交 : 1
* 平行 : 0
* 重合 : -1
*/
inline int relation(Line &a,Line &b)
{
if(a.include(b.a) && a.include(b.b)) return -1;
else if(dcmp(cross(a.b-a.a,b.b-b.a))) return 0;
return 1;
}
// 求两直线交点，注意求之前一定要先判定两直线有无交点
inline Pt jiao_dian(Line &a,Line &b)
{
assert(relation(a,b)==1);
double s1=cross(b.a-a.a,b.b-a.a),s2=cross(b.b-a.b,b.a-a.b);
return a.a+(a.b-a.a)*s1/(s1+s2);
}
struct Poly // 多边形
{
vector<Pt> pts; // 顶点按逆时针排列
inline bool include(const Pt &p) // 判断点是否在多边形内 , 射线法
{
int cnt=0;
for(int i=0;i<SZ(pts);i++)
{
Pt a=pts[i],b=pts[(i+1)%SZ(pts)];
if(Seg(a,b).include(p)) return true;
double d1=a.y-p.y,d2=b.y-p.y,crs=cross(a-p,b-p);
if((crs<=0 && d1>0 && d2<0) || (crs>=0 && d1<0 &&d2>0)) cnt++;
}
return cnt&1;
}
inline double area() // 多边形面积 , 三角剖分法
{
double res=0;
for(int i=0;i<SZ(pts);i++)
{
Pt a=pts[i],b=pts[(i+1)%SZ(pts)];
res+=cross(a,b);
}
return res/2;
}
};
// Andrew 算法求凸包时需要的排序
inline bool cmp1(const Pt &a,const Pt &b)
{
if(a.x==b.x) return a.y>b.y;
return a.x<b.x;
}
/*
** Andrew 算法求凸包
** 主要思想 : 用单调栈维护斜率的单调性 , 求出凸包
*/
Pt h[MAXN];
int cnt_dian;
inline void tubao(Pt *p,int n)
{
sort(p+1,p+n+1,cmp1);
int tp=0;
static int st[MAXN];
static bool used[MAXN];
st[++tp]=1;
for(int i=2;i<=n;i++) // 求下凸壳
{
while(tp>=2 && cross(p[st[tp]]-p[st[tp-1]],p[i]-p[st[tp]])<=0)
used[st[tp--]]=false;
used[i]=true;
st[++tp]=i;
}
int tmp=tp;
for(int i=n-1;i>=1;i--) // 求上凸壳
if(!used[i])
{
while(tp>tmp && cross(p[st[tp]]-p[st[tp-1]],p[i]-p[st[tp]])<=0)
used[st[tp--]]=false;
used[i]=true;
st[++tp]=i;
}
for(int i=1;i<=tp;i++)
h[i]=p[st[i]];
cnt_dian=tp-1;
//#define OUT_TUBAO
#ifdef OUT_TUBAO
double ans=0;
for(int i=2;i<=tp;i++)
ans+=len(h[i]-h[i-1]);
printf("%.2lf\n",ans);
#endif
}
// 判断点 p 是否在直线 l 右侧
inline bool onr(const Line &l,const Pt &p)
{
Pt a=l.a,b=l.b; // 再次强调 , 直线是有向直线 , 方向向量为 b-a
return cross(b-p,a-p)>eps;
}
// 求半平面交对直线的极角排序
inline bool cmp2(const Line &x,const Line &y)
{
if(dcmp(x.ang,y.ang)) return cross(y.a-x.a,y.b-x.a)>eps; // 如果极角相等 , 因为是取左侧半平面交 , 所以把左侧的直线排在前面
return x.ang<y.ang;
}
// 增量法 求半平面交
// 求向量左侧半平面的交集
inline void half_section(Line *line,int n)
{
sort(line+1,line+n+1,cmp2);
int head=1,tail=0;
static int q[MAXN];
static Pt dian[MAXN];
for(int i=1;i<=n;i++)
{
if(dcmp(line[i].ang,line[i-1].ang)) continue;
// 注意顺序 , 一定要先弹队尾 , 再弹队首 , 具体原因画一下三条线的图就知道了
while(head<tail && onr(line[i],dian[tail])) tail--; // 最后一个交点在当前加入直线右侧 , 弹出队尾
while(head<tail && onr(line[i],dian[head+1])) head++; // 第一个交点在当前加入直线右侧 , 弹出队首
q[++tail]=i;
if(head<tail)
dian[tail]=jiao_dian(line[q[tail]],line[q[tail-1]]);
}
while(head<tail && onr(line[q[head]],dian[tail])) tail--; // 用队首弹出队尾
dian[tail+1]=jiao_dian(line[q[head]],line[q[tail]]); tail++; // 求队首与队尾的交点
//#define OUT_AREA
#ifdef OUT_AREA
Poly poly;
for(int i=head+1;i<=tail;i++)
poly.pts.push_back(dian[i]);
printf("%.3lf\n",poly.area());
#endif
}
// 旋转卡壳算法 求凸包直径
// 传入的 p[] 必须是已经逆时针排好序之后的凸包的顶点
inline void get_Mx(Pt *p,int n)
{
double Mx=0;
int cur=2;
p[n+1]=p[1];
for(int i=1;i<=n;i++)
{
while(cross(p[i]-p[cur],p[i+1]-p[cur])<cross(p[i]-p[cur+1],p[i+1]-p[cur+1]))
{
cur++;
if(cur>n) cur=1;
}
//Mx=max(Mx,len(p[i]-p[cur]));
Mx=max(Mx,len_sq(p[i]-p[cur]));
}
//#define OUT_MAX
#ifdef OUT_MAX
// printf("%.3lf\n",Mx);
printf("%lld\n",(ll)Mx);
#endif
}
int main()
{
return 0;
}