计算几何模板
几何
点
struct Point {
double x, y;
Point() : x(0), y(0) {}
Point(double _x, double _y) : x(_x), y(_y) {}
Point operator+(Point& p) { return { x + p.x, y + p.y }; }
Point operator-(Point& p) { return { x - p.x, y - p.y }; }
double operator*(Point& p) { return x * p.x + y * p.y; } //点积
double operator^(Point& p) { return x * p.y - y * p.x; } //叉积
};
误差消除
const double eps=1e-9;
int dcmp(double x){
return (fabs(x)<=eps)?0:(x<0?-1:1);
}
直线
struct Line{
//ax+by+c=0
double a,b,c;
//u为直线上一点，v为方向向量
Point u,v;
Line(){}
//两点确定的直线方程
Line(Point p,Point q){
a=p.y-q.y;b=q.x-p.x;c=p.x*q.y-q.x*p.y;
//保证u、v两点逆时针排列
if((p^q)<0) swap(p,q);
u=p;v=q-p;
}
};
点到直线的距离
double disPL(Point u,Line l){
double length;
length = abs(l.a*u.x+l.b*u.y+l.c)/(sqrt(l.a*l.a+l.b*l.b));
return length;
}
点在直线上的投影点
point proPointLine(Point u,Line l){
point v;
double t=(-u.x*l.a-u.y*l.b-l.c)/(l.a*l.a+l.b*l.b);
v.x=u.x+l.a*t;
v.y=u.y+l.b*t;
return v;
}
求两直线的交点
Point itsLineLine(Line l1,Line l2){
Point p;
double k = l1.a*l2.b-l1.b*l2.a;
p.x = -(l1.c*l2.b-l1.b*l2.c)/k;
p.y = -(l1.a*l2.c-l1.c*l2.a)/k;
return p;
}
求三角形面积
//求三角形面积
double triarea(Point u,Point v,Point w){
//叉积方法
return abs((v-u)^(w-u))/2.0;
}
double triarea(double a,double b,double c){
//海伦公式
double p=(a+b+c)/2;
return sqrt(p*(p-a)*(p-b)*(p-c));
}
求多边形的面积
//求多边形的面积
double polygonArea(point *u,int size){
double area=0;
point begin=u[1];
/*
由第一个点起始的顺序叉积，其中，点可以无序，
面积值为边之间连线的封闭部分，叉积能够计算容斥部分
*/
for(int i=3;i<=size;i++) area+=((u[i-1]-begin)^(u[i]-begin))/2;
return area;
}
圆
//圆结构
struct circle{
//圆心
point cc;
//半径
double radius;
};
//求三点所确定的圆
circle concyclic(point u,point v,point w){
circle c;
point o;
double k=2*(v.x-u.x)*(w.y-v.y)-2*(v.y-u.y)*(w.x-v.x);
o.x=(w.y-v.y)*(v.x*v.x+v.y*v.y-u.x*u.x-u.y*u.y)-(v.y-u.y)*(w.x*w.x+w.y*w.y-v.x*v.x-v.y*v.y);
o.y=(v.x-u.x)*(w.x*w.x+w.y*w.y-v.x*v.x-v.y*v.y)-(w.x-v.x)*(v.x*v.x+v.y*v.y-u.x*u.x-u.y*u.y);
o.x/=k;o.y/=k;
c.cc=o;c.radius=disPointPoint(o,u);
return c;
}
//求圆与直线的交点
vector<point> itsStrCir(line l,circle c){
double k=l.u*l.v,a=norm(l.u),b=norm(l.v),r=c.radius;
double d=k*k-b*b*(a*a-r*r);
vector<point>ans;
if(d==0) ans.push_back(l.u+l.v*(-k/(b*b)));
else {
ans.push_back(l.u+l.v*((-k+d)/(b*b)));
ans.push_back(l.u+l.v*((-k-d)/(b*b)));
}
return ans;
}
//求两圆的交点
vector<point> itsCirCir(circle c1,circle c2){
vector<point> ans;
point o1=c1.cc,o2=c2.cc,a=o2-o1,b;
b.x=a.y;b.y=-a.x;
double r1=c1.radius,r2=c2.radius;
double d=disPointPoint(o1,o2);
double S=triarea(r1,r2,d);
double h=2*S/d;
double t=sqrt(r1*r1-h*h);
//目标两点的中垂线与O1O2反向
if(r1*r1+d*d<r2*r2) t=-t;
//不能用半径和来判断，否则会漏掉内切的情况
if(h==0) ans.push_back(o1+a*t/norm(a));
else {
ans.push_back(o1+a*t/norm(a)+b*h/norm(b));
ans.push_back(o1+a*t/norm(a)-b*h/norm(b));
}
return ans;
}
//求一点与圆的切线
vector<line> tlPointCircle(point u,circle c){
vector<line> ans;
circle o;
o.cc=(c.cc+u)/2;
o.radius=disPointPoint(c.cc,u)/2;
vector<point> p=itsCirCir(o,c);
if(p.size()==1) {
point v;
v.x=(u-c.cc).y;v.y=-(u-c.cc).x;
ans.push_back(line(u,u+v));
}
if(p.size()==2){
ans.push_back(line(p[0],u));
ans.push_back(line(p[1],u));
}
return ans;
}
//求两圆的公切线
vector<line> comTangent(circle c1,circle c2){
vector<line> ans,q;
int r1=c1.radius,r2=c2.radius;
int d=disPointPoint(c1.cc,c2.cc);
point u,v,a=c2.cc-c1.cc,t;
if(r1==r2){
u=c1.cc-c2.cc;
v.x=u.y;v.y=-u.x;
ans.push_back(line(c1.cc+v*r1/norm(v),c1.cc+v*r1/norm(v)+u));
ans.push_back(line(c1.cc-v*r1/norm(v),c1.cc-v*r1/norm(v)+u));
}
else {
//内侧切线
if(triarea(r1,r2,d)==0){
t=c1.cc+a*r1/r2;
q=tlPointCircle(t,c1);
while(q.size()) ans.push_back(q.back()),q.pop_back();
}
//外侧切线
t=c1.cc+a*r1/(r1-r2);
q=tlPointCircle(t,c1);
while(q.size()) ans.push_back(q.back()),q.pop_back();
}
return ans;
}
//最小圆覆盖
circle Smallestcir(point *u,int size){
random_shuffle(u+1,u+1+size);
point o=u[1];
double r=0;
for(int i=2;i<=size;i++){
if(disPointPoint(o,u[i])<=r) continue;
o=(u[i]+u[1])/2;
r=disPointPoint(u[i],u[1])/2;
for(int j=2;j<i;j++){
if(disPointPoint(u[j],o)<=r) continue;
o=(u[i]+u[j])/2;
r=disPointPoint(u[i],u[j])/2;
for(int k=1;k<j;k++){
if(disPointPoint(u[k],o)<=r) continue;
circle c=concyclic(u[i],u[j],u[k]);
o=c.cc;r=c.radius;
}
}
}
circle c;
c.cc=o;c.radius=r;
return c;
}
凸包
圈奶牛
Andrew法
手写栈，先求求下凸包，后求上凸包
int stk[N << 1]; // 记得开双倍
vector<Point> andrew(Point* poly,int n) {
sort(poly + 1, poly + n + 1, [](Point u, Point v) {return (u.x == v.x ? u.y > v.y : u.x < v.x);});
int top = 1;
for (int i = 1; i <= n; i++) {
while (top >= 3 && ((poly[stk[top - 1]] - poly[stk[top - 2]]) ^ (poly[i] - poly[stk[top - 1]])) < 0) {
top--;
}
stk[top++] = i;
}
for (int i = n; i >= 1; i--) {
while (top >= 3 && ((poly[stk[top - 1]] - poly[stk[top - 2]]) ^ (poly[i] - poly[stk[top - 1]])) < 0) {
top--;
}
stk[top++] = i;
}
vector<Point> tmp;
for (int i = 1; i < top - 1; i++) {
tmp.push_back(poly[stk[i]]);
}
return tmp;
}