二维计算几何
点和向量化为坐标Coord进行运算，使用stl中的complex实现。
复数相乘的几何意义为长度相乘，极角相加。
用直线上的一点p和方向向量v表示一条经过p的直线，直线上的所有点q满足q=p+t*v，其中t是参数；当限制t≥0时，该参数方程表示射线；限制0≤t≤1时，该参数方程表示线段。
此外，如果已知线段端点a1和a2，可以通过Line(a1,a2-a1)来得到对应的参数形式。
Morley定理：三角形每个内角的三等分线相交成等边三角形。
欧拉定理：平面图的点数V、边数E和面数F满足V+F-E=2。
typedef double lf;
typedef complex<lf> Coord;
const lf EPS=1e-9,PI=acos(-1);
#define X real()
#define Y imag()
struct Line
{
Coord p,v;
Line(Coord p=Coord(),Coord v=Coord()):p(p),v(v) {}
Coord point(lf t)const
{
return p+v*t;
}
};
struct Circle
{
Coord c;
lf r;
Circle(Coord c=Coord(),lf r=0):c(c),r(r) {}
Coord point(lf t)const//t为参数，幅角
{
return c+polar(r,t);
}
};
/*
Coord(lf x=0,lf y=0);//构造函数
lf real(Coord a);//a的实部（复平面的横坐标）,也可写作a.real()
lf imag(Coord a);//a的虚部（复平面的纵坐标）,也可写作a.imag()
lf abs(Coord a);//向量a的模长，或是点a到原点的距离
lf norm(Coord a);//abs的平方，比abs快，但是要注意浮点数精度溢出
lf arg(Coord a);//a的幅角，与atan2(a.real(),a.imag())等价
Coord polar(lf r,lf t);//极坐标生成方式，r为幅值，t为幅角
//运算符重载+、-、*、/（以及对应的赋值运算，但是赋值运算不能写在表达式中，详见参考地址）、<<、>>（输出括号形式的坐标）
*/
int sgn(lf d)
{
return (d>EPS)-(d<-EPS);
}
bool operator!=(const Coord &A,const Coord &B)//不等运算符，涉及到浮点数比较要重写
{
return sgn(A.X-B.X)||sgn(A.Y-B.Y);
}
bool operator==(const Coord &A,const Coord &B)
{
return !(A!=B);
}
bool cmpCoord(const Coord &A,const Coord &B)//复数没有小于运算，只能这样定义一个比较函数
{
return sgn(A.X-B.X)?sgn(A.X-B.X)<0:sgn(A.Y-B.Y)<0;
}
bool cmpLine(const Line &A,const Line &B)//按极角排序，求凸包中使用
{
return sgn(arg(A.v)-arg(B.v))<0;
}
lf Dot(const Coord &A,const Coord &B)
{
return A.X*B.X+A.Y*B.Y;
}
lf Cross(const Coord &A,const Coord &B)
{
return A.X*B.Y-B.X*A.Y;
}
lf Angle(const Coord &A,const Coord &B)
{
return acos(Dot(A,B)/abs(A)/abs(B));
}
lf Area2(const Coord &A,const Coord &B,const Coord &C)//三角形ABC有向面积的两倍
{
return Cross(B-A,C-A);
}
Coord Rotate(const Coord &A,lf rad)//向量A逆时针旋转rad弧度
{
return A*polar(1.0,rad);
}
Coord Normal(const Coord &A)//A的法向量，把A逆时针旋转九十度并长度化为1
{
lf L=abs(A);
return Coord(-A.Y/L,A.X/L);
}
bool onLeft(const Coord &P,const Line &L)//p是否在有向直线L左侧，不含线上
{
return sgn(Cross(L.v,P-L.p))>0;
}
lf DistanceToLine(const Coord &P,const Line &L)//点到直线距离（有向）
{
return Cross(L.v,P-L.p)/abs(L.v);
}
lf DistanceToLine(const Coord &P,const Coord &A,const Coord &B)
{
return DistanceToLine(P,Line(A,B-A));
}
lf DistanceToSegment(const Coord &P,const Coord &A,const Coord &B)//点到线段的距离（无向）
{
if(A==B)return abs(P-A);
Coord v1=B-A,v2=P-A,v3=P-B;
if(sgn(Dot(v1,v2))<0)return abs(v2);
if(sgn(Dot(v1,v3))>0)return abs(v3);
return fabs(DistanceToLine(P,Line(A,B-A)));
}
Coord getLineProjection(const Coord &P,const Line &L)//点在直线上的投影
{
return L.point(Dot(L.v,P-L.p)/norm(L.v));
}
Coord getLineProjection(const Coord &P,const Coord &A,const Coord &B)
{
return getLineProjection(P,Line(A,B-A));
}
Coord getSymmetry(const Coord &P,const Coord &O)//P关于O的对称点
{
return O+O-P;
}
Coord getSymmetry(const Coord &P,const Line &L)//P关于L的对称点
{
return getSymmetry(P,getLineProjection(P,L));
}
Coord getLineIntersection(const Line &L1,const Line &L2)//直线交点,须确保两直线相交
{
return L1.point(Cross(L2.v,L1.p-L2.p)/Cross(L1.v,L2.v));
}
Coord getLineIntersection(const Coord &A1,const Coord &A2,const Coord &B1,const Coord &B2)
{
return getLineIntersection(Line(A1,A2-A1),Line(B1,B2-B1));
}
bool SegmentProperIntersection(const Coord &A1,const Coord &A2,const Coord &B1,const Coord &B2)//线段相交判定，交点不在一条线段的端点
{
lf C1=Cross(A2-A1,B1-A1),C2=Cross(A2-A1,B2-A1),
C3=Cross(B2-B1,A1-B1),C4=Cross(B2-B1,A2-B1);
return sgn(C1)*sgn(C2)<0&&sgn(C3)*sgn(C4)<0;
}
bool onSegment(const Coord &P,const Coord &A1,const Coord &A2)//判断点是否在线段上，不包含端点
{
return sgn(Dot(A1-P,A2-P))<0&&!sgn(Cross(A1-P,A2-P));
}
lf PolygonArea(const vector<Coord> &p)//计算多边形的有向面积，凸多边形即为面积
{
lf s=0;
for(int i=2; i<p.size(); ++i)
s+=Area2(p[0],p[i-1],p[i]);
return s/2;
}
int inPolygon(const Coord &p,const vector<Coord> &poly)//点在多边形内的判定，转角法，正值为内部，0为外部，-1在边界上
{
int ans=0;
for(int i=0,k,d1,d2,n=poly.size(); i!=n; ++i)
{
if(onSegment(p,poly[i],poly[(i+1)%n]))return -1;//在边界上
k=sgn(Cross(poly[(i+1)%n]-poly[i],p-poly[i]));
d1=sgn(poly[i].Y-p.Y);
d2=sgn(poly[(i+1)%n].Y-p.Y);
if(k>0&&d1<=0&&d2>0)++ans;
if(k<0&&d2<=0&&d1>0)--ans;
}
return ans;
}
vector<Coord> ConvexHull(vector<Coord> p,int collineation=1)//获得凸包，不希望凸包的边上有输入点第二个参数传0
{
vector<Coord> ans;
sort(p.begin(),p.end(),cmpCoord);//先比横坐标再比纵坐标
for(int i=0; i<p.size(); ++i)//求出下凸包
{
while(ans.size()>1&&sgn(Area2(ans[ans.size()-2],ans[ans.size()-1],p[i]))<collineation)
ans.pop_back();
ans.push_back(p[i]);
}
for(int i=p.size()-2,k=ans.size(); i>=0; --i)//求出上凸包
{
while(ans.size()>k&&-sgn(Area2(ans[ans.size()-1],ans[ans.size()-2],p[i]))<collineation)
ans.pop_back();
ans.push_back(p[i]);
}
if(p.size()>1)ans.pop_back();
return ans;
}
vector<Coord> cutPolygon(const vector<Coord> &poly,const Coord &A,const Coord &B)//用有向直线A->B切割多边形poly， 返回“左侧”。 如果退化，可能会返回一个单点或者线段，复杂度O(n^2)
{
vector<Coord> newpoly;
for(int i=0,n=poly.size(); i!=n; ++i)
{
Coord C=poly[i],D=poly[(i+1)%n];
if(sgn(Cross(B-A,C-A))>=0)newpoly.push_back(C);
if(!sgn(Cross(B-A, C-D)))
{
Coord ip=getLineIntersection(Line(A,B-A),Line(C,D-C));
if(onSegment(ip,C,D))newpoly.push_back(ip);
}
}
return newpoly;
}
vector<Coord> getHalfPlaneIntersection(vector<Line> L)//半平面交
{
sort(L.begin(),L.end(),cmpLine);//按极角排序
vector<Coord> p(L.size(),Coord()); //p[i]为q[i]和q[i+1]的交点
int first=0,last=0;//双端队列的第一个元素和最后一个元素
vector<Line> q(L.size(),Line());//双端队列
q[0]=L[0]; //队列初始化为只有一个半平面L[0]
for(int i=0,n=L.size(); i!=n; ++i)
{
while(first<last&&!onLeft(p[last-1],L[i]))
--last;
while(first<last&&!onLeft(p[first],L[i]))
++first;
q[++last]=L[i];
if(!sgn(Cross(q[last].v,q[last-1].v)))
{
--last;
if(onLeft(L[i].p,q[last]))
q[last]=L[i];
}
if(first<last)
p[last-1]=getLineIntersection(q[last-1], q[last]);
}
while(first<last&&!onLeft(p[last-1],q[first]))
--last;//删除无用平面
if(last-first<=1)return vector<Coord>();//空集
p[last]=getLineIntersection(q[last],q[first]);
return vector<Coord>(p.begin()+first,p.begin()+last+1);//从deque复制到输出中
}
int getLineCircleIntersection(const Line &L,const Circle &C,vector<Coord> &sol)
{
lf a=L.v.X,
b=L.p.X-C.c.X,
c=L.v.Y,
d=L.p.Y-C.c.Y,
e=a*a+c*c,
f=2*(a*b+c*d),
g=b*b+d*d-C.r*C.r,
delta=f*f-4*e*g;
if(sgn(delta)<0)return 0;
if(!sgn(delta))
return sol.push_back(L.point(-f/(2*e))),1;
sol.push_back(L.point((-f-sqrt(delta))/(2*e)));
sol.push_back(L.point((-f+sqrt(delta))/(2*e)));
return 2;
}
int getCircleIntersection(const Circle &C1,const Circle &C2,vector<Coord> &sol)
{
lf d=abs(C1.c-C2.c);
if(!sgn(d))
return sgn(C1.r-C2.r)?0:-1;//重合返回-1
if(sgn(C1.r+C2.r-d)<0||sgn(fabs(C1.r-C2.r)-d)>0)//外离或内含
return 0;
lf a=arg(C2.c-C1.c),
da=acos((C1.r*C1.r+d*d-C2.r*C2.r)/(2*C1.r*d));
Coord p1=C1.point(a-da),p2=C1.point(a+da);
sol.push_back(p1);
if(p1==p2)return 1;//相切
return sol.push_back(p2),2;
}
Line getTangent(const Coord &C,const Coord &P)//圆心C，圆上一点P处切线
{
return Line(P,Normal(C-P));
}
int getTangents(const Coord &p,const Circle &C,vector<Coord> &sol)//点到圆的切点，返回个数
{
Coord u=p-C.c;
lf d=abs(u);
if(d<C.r)return 0;//点在圆内
if(!sgn(d-C.r))//点在圆上
return sol.push_back(p),1;
lf base=arg(u),ang=acos(C.r/d);
sol.push_back(C.point(base+ang));
sol.push_back(C.point(base-ang));
return 2;
}
int getTangents(Circle A,Circle &B,vector<Coord> &a,vector<Coord> &b)//公共切线的切点
{
int cnt=0;
if(A.r<B.r)
swap(A,B),swap(a,b);//有时需标记交换
lf d=abs(A.c-B.c),
rdiff=A.r-B.r,
rsum=A.r+B.r;
if(sgn(d-rdiff)<0)return 0;//内含
lf base=arg(B.c-A.c);
if(!sgn(d)&&!sgn(rdiff))return -1;//重合，无穷多条切线
if(!sgn(d-rdiff))//内切，外公切线
{
a.push_back(A.point(base));
b.push_back(B.point(base));
return 1;
}
//有外公切线的情形
lf ang=acos(rdiff/d);
a.push_back(A.point(base+ang));
b.push_back(B.point(base+ang));
a.push_back(A.point(base-ang));
b.push_back(B.point(base-ang));
cnt+=2;
if(!sgn(d-rsum))
{
a.push_back(A.point(base));
b.push_back(B.point(base+PI));
++cnt;
}
else if(sgn(d-rsum)>0)
{
lf ang_in=acos(rsum/d);
a.push_back(A.point(base+ang_in));
b.push_back(B.point(base+ang_in+PI));
a.push_back(A.point(base-ang_in));
b.push_back(B.point(base-ang_in+PI));
cnt+=2;
}
return cnt;
}
lf AreaCircleWithTriangle(const Circle &C,Coord A,Coord B)//C和三角形OAB的相交面积，如果三角形顶点不在O上则把圆和三角形同时平移，直到有一个顶点在O上
{
int sg=sgn(Cross(A,B));
if(!sg||A==C.c||B==C.c)return 0;
lf OA=abs(A-C.c),OB=abs(B-C.c),angle=Angle(A,B),
d=DistanceToLine(Coord(),A,B);
if(sgn(OA-C.r)<=0&&sgn(OB-C.r)<=0)
return Cross(A,B)/2;
if(sgn(OA-C.r)>=0&&sgn(OB-C.r)>=0&&sgn(d-C.r)>=0)
return sg*C.r*C.r*angle/2;
if(sgn(OA-C.r)>=0&&sgn(OB-C.r)>=0&&sgn(d-C.r)<0)
{
Coord prj=getLineProjection(Coord(),A,B);
if(!onSegment(prj,A,B))return sg*C.r*C.r*angle/2;
vector<Coord> p;
Line L=Line(A,B-A);
getLineCircleIntersection(L,C,p);
lf s1=C.r*C.r*angle/2,
s2=C.r*C.r*Angle(p[0],p[1])/2;
s2-=fabs(Cross(p[0],p[1])/2);
s1=s1-s2;
return sg*s1;
}
if(sgn(OB-C.r)<0)swap(A,B);
Line L=Line(A,B-A);
vector<Coord> inter;
getLineCircleIntersection(L,C,inter);
Coord inter_point=inter[!onSegment(inter[0],A,B)];
lf s=fabs(Cross(inter_point,A)/2);
s+=C.r*C.r*Angle(inter_point,B)/2;
return s*sg;
}
lf AreaCircleWithPolygon(const Circle &C,const vector<Coord> &p)
{
lf ans=0;
for(int i=0; i<p.size(); ++i)
ans+=AreaCircleWithTriangle(C,p[i],p[(i+1)%p.size()]);
return fabs(ans);
}
Coord getGravityCenter(const vector<Coord> &p)//多边形重心
{
Coord a(0,0);
lf am=0,mj;
for(int i=0; i<p.size(); ++i)
{
mj=Cross(p[i],p[(i+1)%p.size()]);
a+=mj*(p[i]+p[(i+1)%p.size()]);
am+=mj;
}
return a/am/3.0;
}
/*
struct Area//扫描线求各边平行于坐标轴的矩形面积交并
{
struct Seg
{
ll l,r,h;
int f;//矩形上边缘或下边缘
bool operator<(const Seg &s)const
{
return h<s.h;
}
};
vector<Seg> s;
Ranker rk;
void add(ll l,ll b,ll r,ll t)//(l,b)~(r,t)
{
rk.push_back(l),rk.push_back(r);
s.push_back({l,r,b,1}),s.push_back({l,r,t,0});
}
ll ask(int k)//返回重叠层数至少为k的面积
{
ll ans=0,sum=0;
sort(s.begin(),s.end());
rk.init();
vector<int> vis(rk.size(),0);
for(int i=0; i+1<s.size(); ++i)
{
for(int j=rk.ask(s[i].l),e=rk.ask(s[i].r); j<e; ++j)//此处可线段树维护
{
if(s[i].f&&++vis[j]==k)sum+=rk[j+1]-rk[j];
if(!s[i].f&&vis[j]--==k)sum-=rk[j+1]-rk[j];
}
ans+=sum*(s[i+1].h-s[i].h);
}
return ans;
}
};
*/
1092
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
