2013-team5/geo2d
从 Trac 迁移的文章
这是从旧校内 Wiki 迁移的文章，可能存在一些样式问题，您可以向 memset0 反馈。
原文章内容如下：
{{{
#include <iostream>
#include <cstdio>
#include <cmath>
#include <cstring>
#include <string>
#include <vector>
#include <algorithm>
using namespace std;
const int MAXN = 1000;
const double eps = 1e-8, PI = atan2(0, -1);
const double INF=1e10;
inline double sqr(double x){ return x * x; }
inline bool zero(double x){
return (x > 0 ? x : -x) < eps;
}
inline int sgn(double x){
return (x > eps ? 1 : (x + eps < 0 ? -1 : 0));
}
struct point{
double x, y;
point(double x, double y):x(x), y(y) {}
point() {}
bool operator == (const point & a) const{
return sgn(x - a.x) == 0 && sgn(y - a.y) == 0;
}
bool operator != (const point & a) const{
return sgn(x - a.x) != 0 || sgn(y - a.y) != 0;
}
bool operator < (const point & a) const{
return sgn(x - a.x) < 0 || sgn(x - a.x) == 0 && sgn(y - a.y) < 0;
}
point operator + (const point & a) const{
return point(x + a.x, y + a.y);
}
point operator - (const point & a) const{
return point(x - a.x, y - a.y);
}
point operator * (const double & a) const{
return point(x * a, y * a);
}
point operator / (const double & a) const{
return point(x / a, y / a);
}
double operator * (const point & a) const{
return x * a.y - y * a.x;
}//xmult
double operator ^ (const point & a) const{
return x * a.x + y * a.y;
}//dmult
double sqrlen() const{ return sqr(x) + sqr(y); }
double length() const{ return sqrt(sqrlen()); }
point trunc(double a) const{
return (*this) * (a / length());
}
point rotate(double ang) const{
point p(sin(ang), cos(ang));
return point((*this) * p, (*this) ^ p);
}
point rotate(const point & a) const{
point p(-a.y, a.x); p = p.trunc(1.0);
return point((*this) * p, (*this) ^ p);
}
};
struct line{
point s,t;
double ang;
line() {}
line(point s,point t):s(s), t(t){ ang=atan2(t.y-s.y, t.x-s.x); }
bool operator < (const line &l) const { return ang<l.ang; }
};
typedef vector <point> VP;
/***************** Part one: Dot&Line *****************/
inline bool DotOnSeg(point &p, point &l1, point &l2){
return zero((l2-l1)*(p-l1)) && sgn((l1-p)^(l2-p))<0;
} /* 判断点是否在线段上(不含端点);若包含线段端点,需将<0改成<=0*/ /* tested */
inline double DotToSeg(point &p, point &l1, point &l2){
point L=l2-l1;
if (sgn((p-l1)^L)<=0) return (p-l1).length();/* 1 */
if (sgn((p-l2)^L)> 0) return (p-l2).length();/* 2 */
return fabs(L*(p-l1)/L.length());
} /* 返回点到线段距离;去掉2为点到射线距离;去掉1,2为点到直线距离 */ /* tested */
inline double DotToLine(point &p, point &l1, point &l2){
point L=l2-l1; return fabs(L*(p-l1)/L.length());
} /* 返回点到直线距离 */ /* tested */
inline point DotProLine(point &p, point &l1, point &l2){
point L=l2-l1;
return l1+L*((p-l1)^L)/L.sqrlen();
} /* 返回点在直线上的投影点 */ /* tested */
inline int SideLD(point l1, point l2, point p){
return sgn((l2-l1)*(p-l1));
} /* 判断点在有向直线<l1l2>的哪测,1为左侧,0在线上,-1在右侧 */
inline bool Parallal(point &u1, point &u2, point &v1, point &v2){
return zero((u2-u1)*(v2-v1));
} /* 判断两直线,射线,线段是否平行 */ /* tested */
inline bool SegAcross(point &u1, point &u2, point &v1, point &v2){
return sgn(((u2-u1)*(v1-u1))*((u2-u1)*(v2-u1)))<0
&& sgn(((v2-v1)*(u1-v1))*((v2-v1)*(u2-v1)))<0;
} /* 判断两线段是否规范相交(不包括端点) */
inline point Intersection(point &u1, point &u2, point &v1, point &v2){
double t=((u2-u1)*(v1-u1))/((v2-v1)*(u2-u1));
return v1+(v2-v1)*t;
} /* 使用前需保证两直线,射线或者线段相交,返回其交点 */ /* tested */
/***************** Part two: Polygon *****************/
inline bool DotInsidePoly(point &p, point *v, int n){
int cnt=0;
point t=p; t.y=INF;
for (int i=0; i<n-1 ;i++){
if (SegAcross(p,t,v[i],v[i+1])){ cnt++; continue; }
if (v[i].x<v[i+1].x) cnt+=DotOnSeg(v[i],p,t);
else cnt+=DotOnSeg(v[i+1],p,t);
}
return cnt&1;
} /* 返回点是否严格在多边形内(不包含边界);采用垂直射线法 */
inline double PolyArea(point *v, int n){
double area=0;
for (int i=1; i<n-1; i++) area+=(v[i]-v[0])*(v[i+1]-v[0]);
return fabs(area/2);
} /* 返回多边形面积 */
int ConvexHull (point *p, int n, VP &v){
sort(p,p+n);
point q[n];
int r=0;
for (int i=0; i<n; i++){
while (r>1 && SideLD(q[r-2],q[r-1],p[i])!=1) r--;
q[r++]=p[i];
} /* 从左下向右扫 */
int l=r;
for (int i=n-2; i>=0; i--){
while (r>l && SideLD(q[r-2],q[r-1],p[i])!=1) r--;
q[r++]=p[i];
} /* 从右上向左扫 */
if (n!=1) r--; /* 两次扫都会让左下角的点入队,所以要删去一个 */
for (int i=0; i<r; i++) v.push_back(q[i]);
return r;
} /* 求二维凸包,p为二维平面上的点,凸包顶点存在v内,逆时针方向;SideLD!=1则无三点共线,<0则保留共线点 */ /* tested */
/***************** Part three: Circle *****************/
inline int LineCircleIntersection(point &c, double r, point &l1, point &l2, VP &p){
point t=DotProLine(c,l1,l2);
double r2=r*r, d2=sqr(DotToLine(c,l1,l2));
if (sgn(r2-d2)<0) return 0;
if (zero(r2-d2)){ p.push_back(t); return 1; }
point v=(l1-t)*sqrt((r2-d2)/(l1-t).sqrlen());
p.push_back(t+v);
p.push_back(t-v);
return 2;
} /* 求圆与直线交点,返回值为交点个数,交点存在VP p内 */ /* tested */
inline int CirclesIntersection(point &c1, double r1, point &c2, double r2, VP &p){
double rr=sqr(r1+r2), d2=(c2-c1).sqrlen(), d=sqrt(d2);
if (zero(d2) && zero(r1-r2)) return -1; /* 重合 */
if (sgn(fabs(r1-r2)-d)>0) return 0; /* 内含 */
if (sgn(rr-d2)<0) return 0; /* 相离 */
if (sgn(rr-d2)==0){ p.push_back(c1+(c2-c1).trunc(r1)); return 1; } /* 相切 */
double ang=acos((d2+r1*r1-r2*r2)/(2*d*r1));
point p1=c1+((c2-c1)*r1).rotate(ang)/d;
point p2=c1+((c2-c1)*r1).rotate(-ang)/d;
p.push_back(p1);
p.push_back(p2);
return 2;
} /* 求圆与圆的交点,返回值为交点个数,交点存在VP p内 */ /* tested */
double InscribedCircle(point &p1, point &p2, point &p3, point &c){
double A=(p2-p3).length(), B=(p3-p1).length(), C=(p1-p2).length();
c=(p1*A+p2*B+p3*C)/(A+B+C);
return DotToSeg(c,p1,p2);
} /* 求三角形内切圆,圆心为c,返回值为半径(需满足三点不共线) */ /* tested */
double CircumscribedCircle(point &p1, point &p2, point &p3, point &c){
point B=p2-p1, C=p3-p1;
double S=(B*C)*2.0;
c.x=(C.y*(B.sqrlen())-B.y*(C.sqrlen()))/S+p1.x;
c.y=(B.x*(C.sqrlen())-C.x*(B.sqrlen()))/S+p1.y;
return (c-p1).length();
} /* 求三角形外切圆,圆心为c,返回值为半径(需满足三点不共线) */ /* tested */
int TangentsP(point &p, point &c, double &r, VP &v){
double d=(c-p).length();
if (sgn(d-r)<0) return 0;
if (sgn(d-r)==0){ v.push_back((c-p).rotate(PI/2)); return 1; }
double ang=asin(r/d);
v.push_back((c-p).rotate(ang));
v.push_back((c-p).rotate(-ang));
return 2;
} /* 求所有过p点的圆c的切线,v中为这些切线的方向向量,返回值为切线数 */ /* tested */
/***************** Part four: Plane *****************/
int HalfPlaneIntersection(line *L, int n, VP p){
sort(L,L+n);
point v[n];
line q[n];
q[0]=L[0];
int l=0,r=0;
for (int i=1; i<n; i++){
while (l<r && SideLD(L[i].s,L[i].t,v[r-1])!=1) r--;
while (l<r && SideLD(L[i].s,L[i].t,v[l])!=1) l++;
q[++r]=L[i];
if (Parallal(q[r-1].s,q[r-1].t,q[r].s,q[r].t)){
r--;
if (SideLD(q[r].s,q[r].t,L[i].s)==1) q[r]=L[i];
}
if (l<r) v[r-1]=Intersection(q[r-1].s,q[r-1].t,q[r].s,q[r].t);
}
while (l<r && SideLD(q[l].s,q[l].t,v[r-1])!=1) r--;
if (r-l<=1) return 0;
v[r]=Intersection(q[r].s,q[r].t,q[l].s,q[l].t);
for (int i=l; i<=r; i++) p.push_back(v[i]);
return p.size();
} /* 求半平面交,p存围成相交区域的凸包上的点;若相交区域无界则需添加一个大矩形代替平面 */ /* tested */
}}}#include <iostream>
#include <cstdio>
#include <cmath>
#include <cstring>
#include <string>
#include <vector>
#include <algorithm>
using namespace std;
const int MAXN = 1000;
const double eps = 1e-8, PI = atan2(0, -1);
const double INF=1e10;
inline double sqr(double x){ return x * x; }
inline bool zero(double x){
return (x > 0 ? x : -x) < eps;
}
inline int sgn(double x){
return (x > eps ? 1 : (x + eps < 0 ? -1 : 0));
}
struct point{
double x, y;
point(double x, double y):x(x), y(y) {}
point() {}
bool operator == (const point & a) const{
return sgn(x - a.x) == 0 && sgn(y - a.y) == 0;
}
bool operator != (const point & a) const{
return sgn(x - a.x) != 0 || sgn(y - a.y) != 0;
}
bool operator < (const point & a) const{
return sgn(x - a.x) < 0 || sgn(x - a.x) == 0 && sgn(y - a.y) < 0;
}
point operator + (const point & a) const{
return point(x + a.x, y + a.y);
}
point operator - (const point & a) const{
return point(x - a.x, y - a.y);
}
point operator * (const double & a) const{
return point(x * a, y * a);
}
point operator / (const double & a) const{
return point(x / a, y / a);
}
double operator * (const point & a) const{
return x * a.y - y * a.x;
}//xmult
double operator ^ (const point & a) const{
return x * a.x + y * a.y;
}//dmult
double sqrlen() const{ return sqr(x) + sqr(y); }
double length() const{ return sqrt(sqrlen()); }
point trunc(double a) const{
return (*this) * (a / length());
}
point rotate(double ang) const{
point p(sin(ang), cos(ang));
return point((*this) * p, (*this) ^ p);
}
point rotate(const point & a) const{
point p(-a.y, a.x); p = p.trunc(1.0);
return point((*this) * p, (*this) ^ p);
}
};
struct line{
point s,t;
double ang;
line() {}
line(point s,point t):s(s), t(t){ ang=atan2(t.y-s.y, t.x-s.x); }
bool operator < (const line &l) const { return ang<l.ang; }
};
typedef vector <point> VP;
/***************** Part one: Dot&Line *****************/
inline bool DotOnSeg(point &p, point &l1, point &l2){
return zero((l2-l1)*(p-l1)) && sgn((l1-p)^(l2-p))<0;
} /* 判断点是否在线段上(不含端点);若包含线段端点,需将<0改成<=0*/ /* tested */
inline double DotToSeg(point &p, point &l1, point &l2){
point L=l2-l1;
if (sgn((p-l1)^L)<=0) return (p-l1).length();/* 1 */
if (sgn((p-l2)^L)> 0) return (p-l2).length();/* 2 */
return fabs(L*(p-l1)/L.length());
} /* 返回点到线段距离;去掉2为点到射线距离;去掉1,2为点到直线距离 */ /* tested */
inline double DotToLine(point &p, point &l1, point &l2){
point L=l2-l1; return fabs(L*(p-l1)/L.length());
} /* 返回点到直线距离 */ /* tested */
inline point DotProLine(point &p, point &l1, point &l2){
point L=l2-l1;
return l1+L*((p-l1)^L)/L.sqrlen();
} /* 返回点在直线上的投影点 */ /* tested */
inline int SideLD(point l1, point l2, point p){
return sgn((l2-l1)*(p-l1));
} /* 判断点在有向直线<l1l2>的哪测,1为左侧,0在线上,-1在右侧 */
inline bool Parallal(point &u1, point &u2, point &v1, point &v2){
return zero((u2-u1)*(v2-v1));
} /* 判断两直线,射线,线段是否平行 */ /* tested */
inline bool SegAcross(point &u1, point &u2, point &v1, point &v2){
return sgn(((u2-u1)*(v1-u1))*((u2-u1)*(v2-u1)))<0
&& sgn(((v2-v1)*(u1-v1))*((v2-v1)*(u2-v1)))<0;
} /* 判断两线段是否规范相交(不包括端点) */
inline point Intersection(point &u1, point &u2, point &v1, point &v2){
double t=((u2-u1)*(v1-u1))/((v2-v1)*(u2-u1));
return v1+(v2-v1)*t;
} /* 使用前需保证两直线,射线或者线段相交,返回其交点 */ /* tested */
/***************** Part two: Polygon *****************/
inline bool DotInsidePoly(point &p, point *v, int n){
int cnt=0;
point t=p; t.y=INF;
for (int i=0; i<n-1 ;i++){
if (SegAcross(p,t,v[i],v[i+1])){ cnt++; continue; }
if (v[i].x<v[i+1].x) cnt+=DotOnSeg(v[i],p,t);
else cnt+=DotOnSeg(v[i+1],p,t);
}
return cnt&1;
} /* 返回点是否严格在多边形内(不包含边界);采用垂直射线法 */
inline double PolyArea(point *v, int n){
double area=0;
for (int i=1; i<n-1; i++) area+=(v[i]-v[0])*(v[i+1]-v[0]);
return fabs(area/2);
} /* 返回多边形面积 */
int ConvexHull (point *p, int n, VP &v){
sort(p,p+n);
point q[n];
int r=0;
for (int i=0; i<n; i++){
while (r>1 && SideLD(q[r-2],q[r-1],p[i])!=1) r--;
q[r++]=p[i];
} /* 从左下向右扫 */
int l=r;
for (int i=n-2; i>=0; i--){
while (r>l && SideLD(q[r-2],q[r-1],p[i])!=1) r--;
q[r++]=p[i];
} /* 从右上向左扫 */
if (n!=1) r--; /* 两次扫都会让左下角的点入队,所以要删去一个 */
for (int i=0; i<r; i++) v.push_back(q[i]);
return r;
} /* 求二维凸包,p为二维平面上的点,凸包顶点存在v内,逆时针方向;SideLD!=1则无三点共线,<0则保留共线点 */ /* tested */
/***************** Part three: Circle *****************/
inline int LineCircleIntersection(point &c, double r, point &l1, point &l2, VP &p){
point t=DotProLine(c,l1,l2);
double r2=r*r, d2=sqr(DotToLine(c,l1,l2));
if (sgn(r2-d2)<0) return 0;
if (zero(r2-d2)){ p.push_back(t); return 1; }
point v=(l1-t)*sqrt((r2-d2)/(l1-t).sqrlen());
p.push_back(t+v);
p.push_back(t-v);
return 2;
} /* 求圆与直线交点,返回值为交点个数,交点存在VP p内 */ /* tested */
inline int CirclesIntersection(point &c1, double r1, point &c2, double r2, VP &p){
double rr=sqr(r1+r2), d2=(c2-c1).sqrlen(), d=sqrt(d2);
if (zero(d2) && zero(r1-r2)) return -1; /* 重合 */
if (sgn(fabs(r1-r2)-d)>0) return 0; /* 内含 */
if (sgn(rr-d2)<0) return 0; /* 相离 */
if (sgn(rr-d2)==0){ p.push_back(c1+(c2-c1).trunc(r1)); return 1; } /* 相切 */
double ang=acos((d2+r1*r1-r2*r2)/(2*d*r1));
point p1=c1+((c2-c1)*r1).rotate(ang)/d;
point p2=c1+((c2-c1)*r1).rotate(-ang)/d;
p.push_back(p1);
p.push_back(p2);
return 2;
} /* 求圆与圆的交点,返回值为交点个数,交点存在VP p内 */ /* tested */
double InscribedCircle(point &p1, point &p2, point &p3, point &c){
double A=(p2-p3).length(), B=(p3-p1).length(), C=(p1-p2).length();
c=(p1*A+p2*B+p3*C)/(A+B+C);
return DotToSeg(c,p1,p2);
} /* 求三角形内切圆,圆心为c,返回值为半径(需满足三点不共线) */ /* tested */
double CircumscribedCircle(point &p1, point &p2, point &p3, point &c){
point B=p2-p1, C=p3-p1;
double S=(B*C)*2.0;
c.x=(C.y*(B.sqrlen())-B.y*(C.sqrlen()))/S+p1.x;
c.y=(B.x*(C.sqrlen())-C.x*(B.sqrlen()))/S+p1.y;
return (c-p1).length();
} /* 求三角形外切圆,圆心为c,返回值为半径(需满足三点不共线) */ /* tested */
int TangentsP(point &p, point &c, double &r, VP &v){
double d=(c-p).length();
if (sgn(d-r)<0) return 0;
if (sgn(d-r)==0){ v.push_back((c-p).rotate(PI/2)); return 1; }
double ang=asin(r/d);
v.push_back((c-p).rotate(ang));
v.push_back((c-p).rotate(-ang));
return 2;
} /* 求所有过p点的圆c的切线,v中为这些切线的方向向量,返回值为切线数 */ /* tested */
/***************** Part four: Plane *****************/
int HalfPlaneIntersection(line *L, int n, VP p){
sort(L,L+n);
point v[n];
line q[n];
q[0]=L[0];
int l=0,r=0;
for (int i=1; i<n; i++){
while (l<r && SideLD(L[i].s,L[i].t,v[r-1])!=1) r--;
while (l<r && SideLD(L[i].s,L[i].t,v[l])!=1) l++;
q[++r]=L[i];
if (Parallal(q[r-1].s,q[r-1].t,q[r].s,q[r].t)){
r--;
if (SideLD(q[r].s,q[r].t,L[i].s)==1) q[r]=L[i];
}
if (l<r) v[r-1]=Intersection(q[r-1].s,q[r-1].t,q[r].s,q[r].t);
}
while (l<r && SideLD(q[l].s,q[l].t,v[r-1])!=1) r--;
if (r-l<=1) return 0;
v[r]=Intersection(q[r].s,q[r].t,q[l].s,q[l].t);
for (int i=l; i<=r; i++) p.push_back(v[i]);
return p.size();
} /* 求半平面交,p存围成相交区域的凸包上的点;若相交区域无界则需添加一个大矩形代替平面 */ /* tested */