计算几何模板--更新中
2017-02-03 23:22
363 查看
weeping的私人模板,需要的话就拿去用吧(虽然可能会有小错,如果有请提醒我)
这是现在日常用的模板
2017.10.11 update:
把白书上的半平面交模板换成别的了,原来的那个好像有问题,一直wa。
把求三角形外心的模板改了精度更高的模板。
2017.10.08 update:
原来的模板太啰嗦,今天全部重新简化了,但同时保留以前模板。
2017.09.04 update:
1.发现射线法判断点在多边形内部的代码有误,已修正
2017.08.06 update:
1.卷包裹法求凸包
2.三角形和圆的相关模板
3.O(logn)判断点在凸包内
现在模板:
#include <iostream> #include <cstdio> #include <cmath> #include <algorithm> using namespace std; const double PI = acos(-1.0); const double eps = 1e-10; /****************常用函数***************/ //判断ta与tb的大小关系 int sgn( double ta, double tb) { if(fabs(ta-tb)<eps)return 0; if(ta<tb) return -1; return 1; } //点 class Point { public: double x, y; Point(){} Point( double tx, double ty){ x = tx, y = ty;} bool operator < (const Point &_se) const { return x<_se.x || (x==_se.x && y<_se.y); } friend Point operator + (const Point &_st,const Point &_se) { return Point(_st.x + _se.x, _st.y + _se.y); } friend Point operator - (const Point &_st,const Point &_se) { return Point(_st.x - _se.x, _st.y - _se.y); } //点位置相同(double类型) bool operator == (const Point &_off)const { return sgn(x, _off.x) == 0 && sgn(y, _off.y) == 0; } }; /****************常用函数***************/ //点乘 double dot(const Point &po,const Point &ps,const Point &pe) { return (ps.x - po.x) * (pe.x - po.x) + (ps.y - po.y) * (pe.y - po.y); } //叉乘 double xmult(const Point &po,const Point &ps,const Point &pe) { return (ps.x - po.x) * (pe.y - po.y) - (pe.x - po.x) * (ps.y - po.y); } //两点间距离的平方 double getdis2(const Point &st,const Point &se) { return (st.x - se.x) * (st.x - se.x) + (st.y - se.y) * (st.y - se.y); } //两点间距离 double getdis(const Point &st,const Point &se) { return sqrt((st.x - se.x) * (st.x - se.x) + (st.y - se.y) * (st.y - se.y)); } //两点表示的向量 class Line { public: Point s, e;//两点表示,起点[s],终点[e] double a, b, c;//一般式,ax+by+c=0 double angle;//向量的角度,[-pi,pi] Line(){} Line( Point ts, Point te):s(ts),e(te){}//get_angle();} Line(double _a,double _b,double _c):a(_a),b(_b),c(_c){} //排序用 bool operator < (const Line &ta)const { return angle<ta.angle; } //向量与向量的叉乘 friend double operator / ( const Line &_st, const Line &_se) { return (_st.e.x - _st.s.x) * (_se.e.y - _se.s.y) - (_st.e.y - _st.s.y) * (_se.e.x - _se.s.x); } //向量间的点乘 friend double operator *( const Line &_st, const Line &_se) { return (_st.e.x - _st.s.x) * (_se.e.x - _se.s.x) - (_st.e.y - _st.s.y) * (_se.e.y - _se.s.y); } //从两点表示转换为一般表示 //a=y2-y1,b=x1-x2,c=x2*y1-x1*y2 bool pton() { a = e.y - s.y; b = s.x - e.x; c = e.x * s.y - e.y * s.x; return true; } //半平面交用 //点在向量左边(右边的小于号改成大于号即可,在对应直线上则加上=号) friend bool operator < (const Point &_Off, const Line &_Ori) { return (_Ori.e.y - _Ori.s.y) * (_Off.x - _Ori.s.x) < (_Off.y - _Ori.s.y) * (_Ori.e.x - _Ori.s.x); } //求直线或向量的角度 double get_angle( bool isVector = true) { angle = atan2( e.y - s.y, e.x - s.x); if(!isVector && angle < 0) angle += PI; return angle; } //点在线段或直线上 1:点在直线上 2点在s,e所在矩形内 bool has(const Point &_Off, bool isSegment = false) const { bool ff = sgn( xmult( s, e, _Off), 0) == 0; if( !isSegment) return ff; return ff && sgn(_Off.x - min(s.x, e.x), 0) >= 0 && sgn(_Off.x - max(s.x, e.x), 0) <= 0 && sgn(_Off.y - min(s.y, e.y), 0) >= 0 && sgn(_Off.y - max(s.y, e.y), 0) <= 0; } //点到直线/线段的距离 double dis(const Point &_Off, bool isSegment = false) { ///化为一般式 pton(); //到直线垂足的距离 double td = (a * _Off.x + b * _Off.y + c) / sqrt(a * a + b * b); //如果是线段判断垂足 if(isSegment) { double xp = (b * b * _Off.x - a * b * _Off.y - a * c) / ( a * a + b * b); double yp = (-a * b * _Off.x + a * a * _Off.y - b * c) / (a * a + b * b); double xb = max(s.x, e.x); double yb = max(s.y, e.y); double xs = s.x + e.x - xb; double ys = s.y + e.y - yb; if(xp > xb + eps || xp < xs - eps || yp > yb + eps || yp < ys - eps) td = min( getdis(_Off,s), getdis(_Off,e)); } return fabs(td); } //关于直线对称的点 Point mirror(const Point &_Off) { ///注意先转为一般式 Point ret; double d = a * a + b * b; ret.x = (b * b * _Off.x - a * a * _Off.x - 2 * a * b * _Off.y - 2 * a * c) / d; ret.y = (a * a * _Off.y - b * b * _Off.y - 2 * a * b * _Off.x - 2 * b * c) / d; return ret; } //计算两点的中垂线 static Line ppline(const Point &_a,const Point &_b) { Line ret; ret.s.x = (_a.x + _b.x) / 2; ret.s.y = (_a.y + _b.y) / 2; //一般式 ret.a = _b.x - _a.x; ret.b = _b.y - _a.y; ret.c = (_a.y - _b.y) * ret.s.y + (_a.x - _b.x) * ret.s.x; //两点式 if(fabs(ret.a) > eps) { ret.e.y = 0.0; ret.e.x = - ret.c / ret.a; if(ret.e == ret. s) { ret.e.y = 1e10; ret.e.x = - (ret.c - ret.b * ret.e.y) / ret.a; } } else { ret.e.x = 0.0; ret.e.y = - ret.c / ret.b; if(ret.e == ret. s) { ret.e.x = 1e10; ret.e.y = - (ret.c - ret.a * ret.e.x) / ret.b; } } return ret; } //------------直线和直线(向量)------------- //向量向左边平移t的距离 Line& moveLine( double t) { Point of; of = Point( -( e.y - s.y), e.x - s.x); double dis = sqrt( of.x * of.x + of.y * of.y); of.x= of.x * t / dis, of.y = of.y * t / dis; s = s + of, e = e + of; return *this; } //直线重合 static bool equal(const Line &_st,const Line &_se) { return _st.has( _se.e) && _se.has( _st.s); } //直线平行 static bool parallel(const Line &_st,const Line &_se) { return sgn( _st / _se, 0) == 0; } //两直线(线段)交点 //返回-1代表平行,0代表重合,1代表相交 static bool crossLPt(const Line &_st,const Line &_se, Point &ret) { if(parallel(_st,_se)) { if(Line::equal(_st,_se)) return 0; return -1; } ret = _st.s; double t = ( Line(_st.s,_se.s) / _se) / ( _st / _se); ret.x += (_st.e.x - _st.s.x) * t; ret.y += (_st.e.y - _st.s.y) * t; return 1; } //------------线段和直线(向量)---------- //直线和线段相交 //参数:直线[_st],线段[_se] friend bool crossSL( Line &_st, Line &_se) { return sgn( xmult( _st.s, _se.s, _st.e) * xmult( _st.s, _st.e, _se.e), 0) >= 0; } //判断线段是否相交(注意添加eps) static bool isCrossSS( const Line &_st, const Line &_se) { //1.快速排斥试验判断以两条线段为对角线的两个矩形是否相交 //2.跨立试验(等于0时端点重合) return max(_st.s.x, _st.e.x) >= min(_se.s.x, _se.e.x) && max(_se.s.x, _se.e.x) >= min(_st.s.x, _st.e.x) && max(_st.s.y, _st.e.y) >= min(_se.s.y, _se.e.y) && max(_se.s.y, _se.e.y) >= min(_st.s.y, _st.e.y) && sgn( xmult( _se.s, _st.s, _se.e) * xmult( _se.s, _se.e, _st.s), 0) >= 0 && sgn( xmult( _st.s, _se.s, _st.e) * xmult( _st.s, _st.e, _se.s), 0) >= 0; } }; //寻找凸包的graham 扫描法所需的排序函数 Point gsort; bool gcmp( const Point &ta, const Point &tb)/// 选取与最后一条确定边夹角最小的点,即余弦值最大者 { double tmp = xmult( gsort, ta, tb); if( fabs( tmp) < eps) return getdis( gsort, ta) < getdis( gsort, tb); else if( tmp > 0) return 1; return 0; } class Polygon { public: const static int maxpn = 5e4+7; Point pt[maxpn];//点(顺时针或逆时针) Line dq[maxpn]; //求半平面交打开注释 int n;//点的个数 //求多边形面积,多边形内点必须顺时针或逆时针 double area() { double ans = 0.0; for(int i = 0; i < n; i ++) { int nt = (i + 1) % n; ans += pt[i].x * pt[nt].y - pt[nt].x * pt[i].y; } return fabs( ans / 2.0); } //求多边形重心,多边形内点必须顺时针或逆时针 Point gravity() { Point ans; ans.x = ans.y = 0.0; double area = 0.0; for(int i = 0; i < n; i ++) { int nt = (i + 1) % n; double tp = pt[i].x * pt[nt].y - pt[nt].x * pt[i].y; area += tp; ans.x += tp * (pt[i].x + pt[nt].x); ans.y += tp * (pt[i].y + pt[nt].y); } ans.x /= 3 * area; ans.y /= 3 * area; return ans; } //判断点是否在任意多边形内[射线法],O(n) bool ahas( Point &_Off) { int ret = 0; double infv = 1e20;//坐标系最大范围 Line l = Line( _Off, Point( -infv ,_Off.y)); for(int i = 0; i < n; i ++) { Line ln = Line( pt[i], pt[(i + 1) % n]); if(fabs(ln.s.y - ln.e.y) > eps) { Point tp = (ln.s.y > ln.e.y)? ln.s: ln.e; if( ( fabs( tp.y - _Off.y) < eps && tp.x < _Off.x + eps) || Line::isCrossSS( ln, l)) ret++; } else if( Line::isCrossSS( ln, l)) ret++; } return ret&1; } //判断任意点是否在凸包内,O(logn) bool bhas( Point & p) { if( n < 3) return false; if( xmult( pt[0], p, pt[1]) > eps) return false; if( xmult( pt[0], p, pt[n-1]) < -eps) return false; int l = 2,r = n-1; int line = -1; while( l <= r) { int mid = ( l + r) >> 1; if( xmult( pt[0], p, pt[mid]) >= 0) line = mid,r = mid - 1; else l = mid + 1; } return xmult( pt[line-1], p, pt[line]) <= eps; } //凸多边形被直线分割 Polygon split( Line &_Off) { //注意确保多边形能被分割 Polygon ret; Point spt[2]; double tp = 0.0, np; bool flag = true; int i, pn = 0, spn = 0; for(i = 0; i < n; i ++) { if(flag) pt[pn ++] = pt[i]; else ret.pt[ret.n ++] = pt[i]; np = xmult( _Off.s, _Off.e, pt[(i + 1) % n]); if(tp * np < -eps) { flag = !flag; Line::crossLPt( _Off, Line(pt[i], pt[(i + 1) % n]), spt[spn++]); } tp = (fabs(np) > eps)?np: tp; } ret.pt[ret.n ++] = spt[0]; ret.pt[ret.n ++] = spt[1]; n = pn; return ret; } /** 卷包裹法求点集凸包,_p为输入点集,_n为点的数量 **/ void ConvexClosure( Point _p[], int _n) { sort( _p, _p + _n); n = 0; for(int i = 0; i < _n; i++) { while( n > 1 && sgn( xmult( pt[n-2], pt[n-1], _p[i]), 0) <= 0) n--; pt[n++] = _p[i]; } int _key = n; for(int i = _n - 2; i >= 0; i--) { while( n > _key && sgn( xmult( pt[n-2], pt[n-1], _p[i]), 0) <= 0) n--; pt[n++] = _p[i]; } if(n>1) n--;//除去重复的点,该点已是凸包凸包起点 } /****** 寻找凸包的graham 扫描法********************/ /****** _p为输入的点集,_n为点的数量****************/ void graham( Point _p[], int _n) { int cur=0; for(int i = 1; i < _n; i++) if( sgn( _p[cur].y, _p[i].y) > 0 || ( sgn( _p[cur].y, _p[i].y) == 0 && sgn( _p[cur].x, _p[i].x) > 0) ) cur = i; swap( _p[cur], _p[0]); n = 0, gsort = pt[n++] = _p[0]; if( _n <= 1) return; sort( _p + 1, _p+_n ,gcmp); pt[n++] = _p[1]; for(int i = 2; i < _n; i++) { while(n>1 && sgn( xmult( pt[n-2], pt[n-1], _p[i]), 0) <= 0)// 当凸包退化成直线时需特别注意n n--; pt[n++] = _p[i]; } } //凸包旋转卡壳(注意点必须顺时针或逆时针排列) //返回值凸包直径的平方(最远两点距离的平方) pair<Point,Point> rotating_calipers() { int i = 1 % n; double ret = 0.0; pt = pt[0]; pair<Point,Point>ans=make_pair(pt[0],pt[0]); for(int j = 0; j < n; j ++) { while( fabs( xmult( pt[i+1], pt[j], pt[j + 1])) > fabs( xmult( pt[i], pt[j], pt[j + 1])) + eps) i = (i + 1) % n; //pt[i]和pt[j],pt[i + 1]和pt[j + 1]可能是对踵点 if(ret < getdis2(pt[i],pt[j])) ret = getdis2(pt[i],pt[j]), ans = make_pair(pt[i],pt[j]); if(ret < getdis2(pt[i+1],pt[j+1])) ret = getdis(pt[i+1],pt[j+1]), ans = make_pair(pt[i+1],pt[j+1]); } return ans; } //凸包旋转卡壳(注意点必须逆时针排列) //返回值两凸包的最短距离 double rotating_calipers( Polygon &_Off) { int i = 0; double ret = 1e10;//inf pt = pt[0]; _Off.pt[_Off.n] = _Off.pt[0]; //注意凸包必须逆时针排列且pt[0]是左下角点的位置 while( _Off.pt[i + 1].y > _Off.pt[i].y) i = (i + 1) % _Off.n; for(int j = 0; j < n; j ++) { double tp; //逆时针时为 >,顺时针则相反 while((tp = xmult(_Off.pt[i + 1],pt[j], pt[j + 1]) - xmult(_Off.pt[i], pt[j], pt[j + 1])) > eps) i = (i + 1) % _Off.n; //(pt[i],pt[i+1])和(_Off.pt[j],_Off.pt[j + 1])可能是最近线段 ret = min(ret, Line(pt[j], pt[j + 1]).dis(_Off.pt[i], true)); ret = min(ret, Line(_Off.pt[i], _Off.pt[i + 1]).dis(pt[j + 1], true)); if(tp > -eps)//如果不考虑TLE问题最好不要加这个判断 { ret = min(ret, Line(pt[j], pt[j + 1]).dis(_Off.pt[i + 1], true)); ret = min(ret, Line(_Off.pt[i], _Off.pt[i + 1]).dis(pt[j], true)); } } return ret; } //-----------半平面交------------- //复杂度:O(nlog2(n)) //获取半平面交的多边形(多边形的核) //参数:向量集合[l],向量数量[ln];(半平面方向在向量左边) //函数运行后如果n[即返回多边形的点数量]为0则不存在半平面交的多边形(不存在区域或区域面积无穷大) int judege( Line &_lx, Line &_ly, Line &_lz) { Point tmp; Line::crossLPt(_lx,_ly,tmp); return sgn(xmult(_lz.s,tmp,_lz.e),0); } int halfPanelCross(Line L[], int ln) { int i, tn, bot, top; for(int i = 0; i < ln; i++) L[i].get_angle(); sort(L, L + ln); //平面在向量左边的筛选 for(i = tn = 1; i < ln; i ++) if(fabs(L[i].angle - L[i - 1].angle) > eps) L[tn ++] = L[i]; ln = tn, n = 0, bot = 0, top = 1; dq[0] = L[0], dq[1] = L[1]; for(i = 2; i < ln; i ++) { while(bot < top && judege(dq[top],dq[top-1],L[i]) > 0) top --; while(bot < top && judege(dq[bot],dq[bot+1],L[i]) > 0) bot ++; dq[++ top] = L[i]; } while(bot < top && judege(dq[top],dq[top-1],dq[bot]) > 0) top --; while(bot < top && judege(dq[bot],dq[bot+1],dq[top]) > 0) bot ++; //若半平面交退化为点或线 // if(top <= bot + 1) // return 0; dq[++top] = dq[bot]; for(i = bot; i < top; i ++) Line::crossLPt(dq[i],dq[i + 1],pt[n++]); return n; } }; class Circle { public: Point c;//圆心 double r;//半径 double db, de;//圆弧度数起点, 圆弧度数终点(逆时针0-360) //-------圆--------- //判断圆在多边形内 bool inside( Polygon &_Off) { if(_Off.ahas(c) == false) return false; for(int i = 0; i < _Off.n; i ++) { Line l = Line(_Off.pt[i], _Off.pt[(i + 1) % _Off.n]); if(l.dis(c, true) < r - eps) return false; } return true; } //判断多边形在圆内(线段和折线类似) bool has( Polygon &_Off) { for(int i = 0; i < _Off.n; i ++) if( getdis2(_Off.pt[i],c) > r * r - eps) return false; return true; } //-------圆弧------- //圆被其他圆截得的圆弧,参数:圆[_Off] Circle operator-(Circle &_Off) const { //注意圆必须相交,圆心不能重合 double d2 = getdis2(c,_Off.c); double d = getdis(c,_Off.c); double ans = acos((d2 + r * r - _Off.r * _Off.r) / (2 * d * r)); Point py = _Off.c - c; double oans = atan2(py.y, py.x); Circle res; res.c = c; res.r = r; res.db = oans + ans; res.de = oans - ans + 2 * PI; return res; } //圆被其他圆截得的圆弧,参数:圆[_Off] Circle operator+(Circle &_Off) const { //注意圆必须相交,圆心不能重合 double d2 = getdis2(c,_Off.c); double d = getdis(c,_Off.c); double ans = acos((d2 + r * r - _Off.r * _Off.r) / (2 * d * r)); Point py = _Off.c - c; double oans = atan2(py.y, py.x); Circle res; res.c = c; res.r = r; res.db = oans - ans; res.de = oans + ans; return res; } //过圆外一点的两条切线 //参数:点[_Off](必须在圆外),返回:两条切线(切线的s点为_Off,e点为切点) pair<Line, Line> tangent( Point &_Off) { double d = getdis(c,_Off); //计算角度偏移的方式 double angp = acos(r / d), ango = atan2(_Off.y - c.y, _Off.x - c.x); Point pl = Point(c.x + r * cos(ango + angp), c.y + r * sin(ango + angp)), pr = Point(c.x + r * cos(ango - angp), c.y + r * sin(ango - angp)); return make_pair(Line(_Off, pl), Line(_Off, pr)); } //计算直线和圆的两个交点 //参数:直线[_Off](两点式),返回两个交点,注意直线必须和圆有两个交点 pair<Point, Point> cross(Line _Off) { _Off.pton(); //到直线垂足的距离 double td = fabs(_Off.a * c.x + _Off.b * c.y + _Off.c) / sqrt(_Off.a * _Off.a + _Off.b * _Off.b); //计算垂足坐标 double xp = (_Off.b * _Off.b * c.x - _Off.a * _Off.b * c.y - _Off.a * _Off.c) / ( _Off.a * _Off.a + _Off.b * _Off.b); double yp = (- _Off.a * _Off.b * c.x + _Off.a * _Off.a * c.y - _Off.b * _Off.c) / (_Off.a * _Off.a + _Off.b * _Off.b); double ango = atan2(yp - c.y, xp - c.x); double angp = acos(td / r); return make_pair(Point(c.x + r * cos(ango + angp), c.y + r * sin(ango + angp)), Point(c.x + r * cos(ango - angp), c.y + r * sin(ango - angp))); } }; class triangle { public: Point a, b, c;//顶点 triangle(){} triangle(Point a, Point b, Point c): a(a), b(b), c(c){} //计算三角形面积 double area() { return fabs( xmult(a, b, c)) / 2.0; } //计算三角形外心 //返回:外接圆圆心 Point circumcenter() { double pa = a.x * a.x + a.y * a.y; double pb = b.x * b.x + b.y * b.y; double pc = c.x * c.x + c.y * c.y; double ta = pa * ( b.y - c.y) - pb * ( a.y - c.y) + pc * ( a.y - b.y); double tb = -pa * ( b.x - c.x) + pb * ( a.x - c.x) - pc * ( a.x - b.x); double tc = a.x * ( b.y - c.y) - b.x * ( a.y - c.y) + c.x * ( a.y - b.y); return Point( ta / 2.0 / tc, tb / 2.0 / tc); } //计算三角形内心 //返回:内接圆圆心 Point incenter() { Line u, v; double m, n; u.s = a; m = atan2(b.y - a.y, b.x - a.x); n = atan2(c.y - a.y, c.x - a.x); u.e.x = u.s.x + cos((m + n) / 2); u.e.y = u.s.y + sin((m + n) / 2); v.s = b; m = atan2(a.y - b.y, a.x - b.x); n = atan2(c.y - b.y, c.x - b.x); v.e.x = v.s.x + cos((m + n) / 2); v.e.y = v.s.y + sin((m + n) / 2); Point ret; Line::crossLPt(u,v,ret); return ret; } //计算三角形垂心 //返回:高的交点 Point perpencenter() { Line u,v; u.s = c; u.e.x = u.s.x - a.y + b.y; u.e.y = u.s.y + a.x - b.x; v.s = b; v.e.x = v.s.x - a.y + c.y; v.e.y = v.s.y + a.x - c.x; Point ret; Line::crossLPt(u,v,ret); return ret; } //计算三角形重心 //返回:重心 //到三角形三顶点距离的平方和最小的点 //三角形内到三边距离之积最大的点 Point barycenter() { Line u,v; u.s.x = (a.x + b.x) / 2; u.s.y = (a.y + b.y) / 2; u.e = c; v.s.x = (a.x + c.x) / 2; v.s.y = (a.y + c.y) / 2; v.e = b; Point ret; Line::crossLPt(u,v,ret); return ret; } //计算三角形费马点 //返回:到三角形三顶点距离之和最小的点 Point fermentPoint() { Point u, v; double step = fabs(a.x) + fabs(a.y) + fabs(b.x) + fabs(b.y) + fabs(c.x) + fabs(c.y); int i, j, k; u.x = (a.x + b.x + c.x) / 3; u.y = (a.y + b.y + c.y) / 3; while (step > eps) { for (k = 0; k < 10; step /= 2, k ++) { for (i = -1; i <= 1; i ++) { for (j =- 1; j <= 1; j ++) { v.x = u.x + step * i; v.y = u.y + step * j; if (getdis(u,a) + getdis(u,b) + getdis(u,c) > getdis(v,a) + getdis(v,b) + getdis(v,c)) u = v; } } } } return u; } }; int main(void) { return 0; }
以前模板:
#include <iostream> #include <cstdio> #include <cmath> #include <algorithm> using namespace std; const double PI = acos(-1.0); const double eps = 1e-10; //点 class Point { public: double x, y; Point(){} Point(double x, double y):x(x),y(y){} bool operator < (const Point &_se) const { return x<_se.x || (x==_se.x && y<_se.y); } /*******判断ta与tb的大小关系*******/ static int sgn(double ta,double tb) { if(fabs(ta-tb)<eps)return 0; if(ta<tb) return -1; return 1; } static double xmult(const Point &po, const Point &ps, const Point &pe) { return (ps.x - po.x) * (pe.y - po.y) - (pe.x - po.x) * (ps.y - po.y); } friend Point operator + (const Point &_st,const Point &_se) { return Point(_st.x + _se.x, _st.y + _se.y); } friend Point operator - (const Point &_st,const Point &_se) { return Point(_st.x - _se.x, _st.y - _se.y); } //点位置相同(double类型) bool operator == (const Point &_off) const { return Point::sgn(x, _off.x) == 0 && Point::sgn(y, _off.y) == 0; } //点位置不同(double类型) bool operator != (const Point &_Off) const { return ((*this) == _Off) == false; } //两点间距离的平方 static double dis2(const Point &_st,const Point &_se) { return (_st.x - _se.x) * (_st.x - _se.x) + (_st.y - _se.y) * (_st.y - _se.y); } //两点间距离 static double dis(const Point &_st, const Point &_se) { return sqrt((_st.x - _se.x) * (_st.x - _se.x) + (_st.y - _se.y) * (_st.y - _se.y)); } }; //两点表示的向量 class Line { public: Point s, e;//两点表示,起点[s],终点[e] double a, b, c;//一般式,ax+by+c=0 double angle;//向量的角度,[-pi,pi] Line(){} Line(const Point &s, const Point &e):s(s),e(e){get_angle(1);} Line(double _a,double _b,double _c):a(_a),b(_b),c(_c){} //向量与点的叉乘,参数:点[_Off] //[点相对向量位置判断] double operator /(const Point &_Off) const { return (_Off.y - s.y) * (e.x - s.x) - (_Off.x - s.x) * (e.y - s.y); } //向量与向量的叉乘,参数:向量[_Off] friend double operator /(const Line &_st,const Line &_se) { return (_st.e.x - _st.s.x) * (_se.e.y - _se.s.y) - (_st.e.y - _st.s.y) * (_se.e.x - _se.s.x); } friend double operator *(const Line &_st,const Line &_se) { return (_st.e.x - _st.s.x) * (_se.e.x - _se.s.x) - (_st.e.y - _st.s.y) * (_se.e.y - _se.s.y); } //从两点表示转换为一般表示 //a=y2-y1,b=x1-x2,c=x2*y1-x1*y2 bool pton() { a = e.y - s.y; b = s.x - e.x; c = e.x * s.y - e.y * s.x; return true; } //求直线或向量的角度 double get_angle(bool isVector) { angle=atan2(e.y-s.y,e.x-s.x); if(!isVector && angle<0) angle+=PI; return angle; } // bool operator < (const Line &ta)const { return angle<ta.angle; } //-----------点和直线(向量)----------- //点在向量左边(右边的小于号改成大于号即可,在对应直线上则加上=号) //参数:点[_Off],向量[_Ori] friend bool operator<(const Point &_Off, const Line &_Ori) { return (_Ori.e.y - _Ori.s.y) * (_Off.x - _Ori.s.x) < (_Off.y - _Ori.s.y) * (_Ori.e.x - _Ori.s.x); } //点在直线上,参数:点[_Off] bool lhas(const Point &_Off) const { return Point::sgn((*this) / _Off, 0) == 0; } //点在线段上,参数:点[_Off] bool shas(const Point &_Off) const { return lhas(_Off) && Point::sgn(_Off.x - min(s.x, e.x), 0) > 0 && Point::sgn(_Off.x - max(s.x, e.x), 0) < 0 && Point::sgn(_Off.y - min(s.y, e.y), 0) > 0 && Point::sgn(_Off.y - max(s.y, e.y), 0) < 0; } //点到直线/线段的距离 //参数: 点[_Off], 是否是线段[isSegment](默认为直线) double dis(const Point &_Off, bool isSegment = false) { ///化为一般式 pton(); //到直线垂足的距离 double td = (a * _Off.x + b * _Off.y + c) / sqrt(a * a + b * b); //如果是线段判断垂足 if(isSegment) { double xp = (b * b * _Off.x - a * b * _Off.y - a * c) / ( a * a + b * b); double yp = (-a * b * _Off.x + a * a * _Off.y - b * c) / (a * a + b * b); double xb = max(s.x, e.x); double yb = max(s.y, e.y); double xs = s.x + e.x - xb; double ys = s.y + e.y - yb; if(xp > xb + eps || xp < xs - eps || yp > yb + eps || yp < ys - eps) td = min(Point::dis(_Off,s), Point::dis(_Off,e)); } return fabs(td); } //关于直线对称的点 Point mirror(const Point &_Off) const { ///注意先转为一般式 Point ret; double d = a * a + b * b; ret.x = (b * b * _Off.x - a * a * _Off.x - 2 * a * b * _Off.y - 2 * a * c) / d; ret.y = (a * a * _Off.y - b * b * _Off.y - 2 * a * b * _Off.x - 2 * b * c) / d; return ret; } //计算两点的中垂线 static Line ppline(const Point &_a, const Point &_b) { Line ret; ret.s.x = (_a.x + _b.x) / 2; ret.s.y = (_a.y + _b.y) / 2; //一般式 ret.a = _b.x - _a.x; ret.b = _b.y - _a.y; ret.c = (_a.y - _b.y) * ret.s.y + (_a.x - _b.x) * ret.s.x; //两点式 if(std::fabs(ret.a) > eps) { ret.e.y = 0.0; ret.e.x = - ret.c / ret.a; if(ret.e == ret. s) { ret.e.y = 1e10; ret.e.x = - (ret.c - ret.b * ret.e.y) / ret.a; } } else { ret.e.x = 0.0; ret.e.y = - ret.c / ret.b; if(ret.e == ret. s) { ret.e.x = 1e10; ret.e.y = - (ret.c - ret.a * ret.e.x) / ret.b; } } return ret; } //------------直线和直线(向量)------------- //直线向左边平移t的距离 Line& moveLine(double t) { Point of; of=Point(-(e.y-s.y),e.x-s.x); double dis=sqrt(of.x*of.x+of.y*of.y); of.x=of.x*t/dis,of.y=of.y*t/dis; s=s+of,e=e+of; return *this; } //直线重合,参数:直线向量[_st],[_se] static bool equal(const Line &_st, const Line &_se) { return _st.lhas(_se.e) && _se.lhas(_se.s); } //直线平行,参数:直线向量[_st],[_se] static bool parallel(const Line &_st,const Line &_se) { return Point::sgn(_st / _se, 0) == 0; } //两直线(线段)交点,参数:直线向量[_st],[_se],交点 //返回-1代表平行,0代表重合,1代表相交 static bool crossLPt(const Line &_st,const Line &_se,Point &ret) { if(Line::parallel(_st,_se)) { if(Line::equal(_st,_se)) return 0; return -1; } ret = _st.s; double t = (Line(_st.s,_se.s)/_se)/(_st/_se); ret.x += (_st.e.x - _st.s.x) * t; ret.y += (_st.e.y - _st.s.y) * t; return 1; } //------------线段和直线(向量)---------- //线段和直线交 //参数:直线[_st],线段[_se] friend bool crossSL(const Line &_st,const Line &_se) { return Point::sgn((_st / _se.s) * (_st / _se.e) ,0) <= 0; } //------------线段和线段(向量)---------- //判断线段是否相交(注意添加eps),参数:线段[_st],线段[_se] static bool isCrossSS(const Line &_st,const Line &_se) { //1.快速排斥试验判断以两条线段为对角线的两个矩形是否相交 //2.跨立试验(等于0时端点重合) return max(_st.s.x, _st.e.x) >= min(_se.s.x, _se.e.x) && max(_se.s.x, _se.e.x) >= min(_st.s.x, _st.e.x) && max(_st.s.y, _st.e.y) >= min(_se.s.y, _se.e.y) && max(_se.s.y, _se.e.y) >= min(_st.s.y, _st.e.y) && Point::sgn((_st / Line(_st.s, _se.s)) * (_st / Line(_st.s, _se.e)), 0) <= 0 && Point::sgn((_se / Line(_se.s, _st.s)) * (_se / Line(_se.s, _st.e)), 0) <= 0; } }; class Polygon { public: const static int maxpn = 5e4+7; Point pt[maxpn];//点(顺时针或逆时针) int n;//点的个数 //求多边形面积,多边形内点必须顺时针或逆时针 double area() const { double ans = 0.0; for(int i = 0; i < n; i ++) { int nt = (i + 1) % n; ans += pt[i].x * pt[nt].y - pt[nt].x * pt[i].y; } return fabs(ans / 2.0); } //求多边形重心,多边形内点必须顺时针或逆时针 Point gravity() const { Point ans; ans.x = ans.y = 0.0; double area = 0.0; for(int i = 0; i < n; i ++) { int nt = (i + 1) % n; double tp = pt[i].x * pt[nt].y - pt[nt].x * pt[i].y; area += tp; ans.x += tp * (pt[i].x + pt[nt].x); ans.y += tp * (pt[i].y + pt[nt].y); } ans.x /= 3 * area; ans.y /= 3 * area; return ans; } //判断点在凸多边形内,参数:点[_Off] bool chas(const Point &_Off) const { double tp = 0, np; for(int i = 0; i < n; i ++) { np = Line(pt[i], pt[(i + 1) % n]) / _Off; if(tp * np < -eps) return false; tp = (fabs(np) > eps)?np: tp; } return true; } //判断点是否在任意多边形内[射线法],O(n) bool ahas(const Point &_Off) const { int ret = 0; double infv = 1e20;//坐标系最大范围 Line l = Line(_Off, Point( -infv ,_Off.y)); for(int i = 0; i < n; i ++) { Line ln = Line(pt[i], pt[(i + 1) % n]); if(fabs(ln.s.y - ln.e.y) > eps) { Point tp = (ln.s.y > ln.e.y)? ln.s: ln.e; if((fabs(tp.y - _Off.y) < eps && tp.x < _Off.x + eps)||Line::isCrossSS(ln,l)) ret ++; } else if(Line::isCrossSS(ln,l)) ret ++; } return ret&1; } //判断任意点是否在凸包内,O(lgn) bool bhas(const Point & p) { if(n<3) return false; if(Point::xmult(pt[0],p,pt[1])>eps) return false; if(Point::xmult(pt[0],p,pt[n-1])<-eps) return false; int l=2,r=n-1; int line=-1; while(l<=r) { int mid=(l+r)>>1; if(Point::xmult(pt[0],p,pt[mid])>=0) line=mid,r=mid-1; else l=mid+1; } return Point::xmult(pt[line-1],p,pt[line])<=eps; } //凸多边形被直线分割,参数:直线[_Off] Polygon split(Line _Off) { //注意确保多边形能被分割 Polygon ret; Point spt[2]; double tp = 0.0, np; bool flag = true; int i, pn = 0, spn = 0; for(i = 0; i < n; i ++) { if(flag) pt[pn ++] = pt[i]; else ret.pt[ret.n ++] = pt[i]; np = _Off / pt[(i + 1) % n]; if(tp * np < -eps) { flag = !flag; Line::crossLPt(_Off,Line(pt[i], pt[(i + 1) % n]),spt[spn++]); } tp = (fabs(np) > eps)?np: tp; } ret.pt[ret.n ++] = spt[0]; ret.pt[ret.n ++] = spt[1]; n = pn; return ret; } /** 卷包裹法求点集凸包,_p为输入点集,_n为点的数量 **/ void ConvexClosure(Point _p[],int _n) { sort(_p,_p+_n); n=0; for(int i=0;i<_n;i++) { while(n>1&&Point::sgn(Line(pt[n-2],pt[n-1])/Line(pt[n-2],_p[i]),0)<=0) n--; pt[n++]=_p[i]; } int _key=n; for(int i=_n-2;i>=0;i--) { while(n>_key&&Point::sgn(Line(pt[n-2],pt[n-1])/Line(pt[n-2],_p[i]),0)<=0) n--; pt[n++]=_p[i]; } if(n>1) n--;//除去重复的点,该点已是凸包凸包起点 } // /****** 寻找凸包的graham 扫描法********************/ // /****** _p为输入的点集,_n为点的数量****************/ // /**使用时需把gmp函数放在Polygon类上面L,ine类下面,并且看情况修改pt[0]**/ // bool gcmp(const Point &ta,const Point &tb)/// 选取与最后一条确定边夹角最小的点,即余弦值最大者 // { // double tmp=Line(pt[0],ta)/Line(pt[0],tb); // if(Point::sgn(tmp,0)==0) // return Point::dis(pt[0],ta)<Point::dis(pt[0],tb); // else if(tmp>0) // return 1; // return 0; // } // void graham(Point _p[],int _n) // { // int cur=0; // for(int i=1;i<_n;i++) // if(Point::sgn(_p[cur].y,_p[i].y)>0 || (Point::sgn(_p[cur].y,_p[i].y)==0 && Point::sgn(_p[cur].x,_p[i].x)>0)) // cur=i; // swap(_p[cur],_p[0]); // n=0,pt[n++]=_p[0]; // if(_n==1) return; // sort(_p+1,_p+_n,Polygon::gcmp); // pt[n++]=_p[1],pt[n++]=_p[2]; // for(int i=3;i<_n;i++) // { // while(n>1 && Point::sgn(Line(pt[n-2],pt[n-1])/Line(pt[n-2],_p[i]),0)<0)// 当凸包退化成直线时需特别注意n // n--; // pt[n++]=_p[i]; // } // } //凸包旋转卡壳(注意点必须顺时针或逆时针排列) //返回值凸包直径的平方(最远两点距离的平方) double rotating_calipers() { int i = 1; double ret = 0.0; pt = pt[0]; for(int j = 0; j < n; j ++) { while(fabs(Point::xmult(pt[i+1],pt[j], pt[j + 1])) > fabs(Point::xmult(pt[i],pt[j], pt[j + 1])) + eps) i = (i + 1) % n; //pt[i]和pt[j],pt[i + 1]和pt[j + 1]可能是对踵点 ret = max(ret, max(Point::dis(pt[i],pt[j]), Point::dis(pt[i + 1],pt[j + 1]))); } return ret; } //凸包旋转卡壳(注意点必须逆时针排列) //返回值两凸包的最短距离 double rotating_calipers(Polygon &_Off) { int i = 0; double ret = 1e10;//inf pt = pt[0]; _Off.pt[_Off.n] = _Off.pt[0]; //注意凸包必须逆时针排列且pt[0]是左下角点的位置 while(_Off.pt[i + 1].y > _Off.pt[i].y) i = (i + 1) % _Off.n; for(int j = 0; j < n; j ++) { double tp; //逆时针时为 >,顺时针则相反 while((tp = Point::xmult(_Off.pt[i + 1],pt[j], pt[j + 1]) - Point::xmult(_Off.pt[i], pt[j], pt[j + 1])) > eps) i = (i + 1) % _Off.n; //(pt[i],pt[i+1])和(_Off.pt[j],_Off.pt[j + 1])可能是最近线段 ret = min(ret, Line(pt[j], pt[j + 1]).dis(_Off.pt[i], true)); ret = min(ret, Line(_Off.pt[i], _Off.pt[i + 1]).dis(pt[j + 1], true)); if(tp > -eps)//如果不考虑TLE问题最好不要加这个判断 { ret = min(ret, Line(pt[j], pt[j + 1]).dis(_Off.pt[i + 1], true)); ret = min(ret, Line(_Off.pt[i], _Off.pt[i + 1]).dis(pt[j], true)); } } return ret; } //-----------半平面交------------- //复杂度:O(nlog2(n)) //#include <algorithm> //获取半平面交的多边形(多边形的核) //参数:向量集合[l],向量数量[ln];(半平面方向在向量左边) //函数运行后如果n[即返回多边形的点数量]为0则不存在半平面交的多边形(不存在区域或区域面积无穷大) int halfPanelCross(Line *L, int ln) { sort(L,L + ln); //极角排序 int first,last; //双端队列的队首和队尾的下标 Point *p = new Point[ln]; //p[i]为q[i]和q[i+1]的交点 Line *q = new Line[ln]; //双端队列 q[first=last=0] = L[0]; //双端队列初始只有一个半平面L[0] for(int i = 1; i < ln; i++) { while(first < last && !(p[last-1] < L[i])) last--; while(first < last && !(p[first] < L[i])) first++; q[++last] = L[i]; if(fabs(q[last]/q[last-1]) < eps) { if(L[i].s < q[--last]) q[last] = L[i]; //两向量平行,取内侧的一个 } if(first < last) Line::crossLPt(q[last-1], q[last], p[last-1]); } while(first < last && !(p[last-1] < q[first])) last--; //删除无用平面 if(last - first <= 1) return 0; //空集 Line::crossLPt(q[last], q[first], p[last]); //计算首尾两个半平面的交点 int m=0; for(int i = 1; i <= last; i++) pt[m++] = p[i]; //复制到pt中 return m; } }; class Circle { public: Point c;//圆心 double r;//半径 double db, de;//圆弧度数起点, 圆弧度数终点(逆时针0-360) //-------圆--------- //判断圆在多边形内 bool inside(const Polygon &_Off) const { if(_Off.ahas(c) == false) return false; for(int i = 0; i < _Off.n; i ++) { Line l = Line(_Off.pt[i], _Off.pt[(i + 1) % _Off.n]); if(l.dis(c, true) < r - eps) return false; } return true; } //判断多边形在圆内(线段和折线类似) bool has(const Polygon &_Off) const { for(int i = 0; i < _Off.n; i ++) if(Point::dis2(_Off.pt[i],c) > r * r - eps) return false; return true; } //-------圆弧------- //圆被其他圆截得的圆弧,参数:圆[_Off] Circle operator-(Circle &_Off) const { //注意圆必须相交,圆心不能重合 double d2 = Point::dis2(c,_Off.c); double d = Point::dis(c,_Off.c); double ans = acos((d2 + r * r - _Off.r * _Off.r) / (2 * d * r)); Point py = _Off.c - c; double oans = atan2(py.y, py.x); Circle res; res.c = c; res.r = r; res.db = oans + ans; res.de = oans - ans + 2 * PI; return res; } //圆被其他圆截得的圆弧,参数:圆[_Off] Circle operator+(Circle &_Off) const { //注意圆必须相交,圆心不能重合 double d2 = Point::dis2(c,_Off.c); double d = Point::dis(c,_Off.c); double ans = acos((d2 + r * r - _Off.r * _Off.r) / (2 * d * r)); Point py = _Off.c - c; double oans = atan2(py.y, py.x); Circle res; res.c = c; res.r = r; res.db = oans - ans; res.de = oans + ans; return res; } //过圆外一点的两条切线 //参数:点[_Off](必须在圆外),返回:两条切线(切线的s点为_Off,e点为切点) std::pair<Line, Line> tangent(const Point &_Off) const { double d = Point::dis(c,_Off); //计算角度偏移的方式 double angp = std::acos(r / d), ango = std::atan2(_Off.y - c.y, _Off.x - c.x); Point pl = Point(c.x + r * std::cos(ango + angp), c.y + r * std::sin(ango + angp)), pr = Point(c.x + r * std::cos(ango - angp), c.y + r * std::sin(ango - angp)); return std::make_pair(Line(_Off, pl), Line(_Off, pr)); } //计算直线和圆的两个交点 //参数:直线[_Off](两点式),返回两个交点,注意直线必须和圆有两个交点 std::pair<Point, Point> cross(Line _Off) const { _Off.pton(); //到直线垂足的距离 double td = fabs(_Off.a * c.x + _Off.b * c.y + _Off.c) / sqrt(_Off.a * _Off.a + _Off.b * _Off.b); //计算垂足坐标 double xp = (_Off.b * _Off.b * c.x - _Off.a * _Off.b * c.y - _Off.a * _Off.c) / ( _Off.a * _Off.a + _Off.b * _Off.b); double yp = (- _Off.a * _Off.b * c.x + _Off.a * _Off.a * c.y - _Off.b * _Off.c) / (_Off.a * _Off.a + _Off.b * _Off.b); double ango = std::atan2(yp - c.y, xp - c.x); double angp = std::acos(td / r); return std::make_pair(Point(c.x + r * std::cos(ango + angp), c.y + r * std::sin(ango + angp)), Point(c.x + r * std::cos(ango - angp), c.y + r * std::sin(ango - angp))); } }; class triangle { public: Point a, b, c;//顶点 triangle(){} triangle(Point a, Point b, Point c): a(a), b(b), c(c){} //计算三角形面积 double area() { return fabs(Point::xmult(a, b, c)) / 2.0; } //计算三角形外心 //返回:外接圆圆心 Point circumcenter() { Line u,v; u.s.x = (a.x + b.x) / 2; u.s.y = (a.y + b.y) / 2; u.e.x = u.s.x - a.y + b.y; u.e.y = u.s.y + a.x - b.x; v.s.x = (a.x + c.x) / 2; v.s.y = (a.y + c.y) / 2; v.e.x = v.s.x - a.y + c.y; v.e.y = v.s.y + a.x - c.x; Point ret; Line::crossLPt(u,v,ret); return ret; } //计算三角形内心 //返回:内接圆圆心 Point incenter() { Line u, v; double m, n; u.s = a; m = atan2(b.y - a.y, b.x - a.x); n = atan2(c.y - a.y, c.x - a.x); u.e.x = u.s.x + cos((m + n) / 2); u.e.y = u.s.y + sin((m + n) / 2); v.s = b; m = atan2(a.y - b.y, a.x - b.x); n = atan2(c.y - b.y, c.x - b.x); v.e.x = v.s.x + cos((m + n) / 2); v.e.y = v.s.y + sin((m + n) / 2); Point ret; Line::crossLPt(u,v,ret); return ret; } //计算三角形垂心 //返回:高的交点 Point perpencenter() { Line u,v; u.s = c; u.e.x = u.s.x - a.y + b.y; u.e.y = u.s.y + a.x - b.x; v.s = b; v.e.x = v.s.x - a.y + c.y; v.e.y = v.s.y + a.x - c.x; Point ret; Line::crossLPt(u,v,ret); return ret; } //计算三角形重心 //返回:重心 //到三角形三顶点距离的平方和最小的点 //三角形内到三边距离之积最大的点 Point barycenter() { Line u,v; u.s.x = (a.x + b.x) / 2; u.s.y = (a.y + b.y) / 2; u.e = c; v.s.x = (a.x + c.x) / 2; v.s.y = (a.y + c.y) / 2; v.e = b; Point ret; Line::crossLPt(u,v,ret); return ret; } //计算三角形费马点 //返回:到三角形三顶点距离之和最小的点 Point fermentPoint() { Point u, v; double step = fabs(a.x) + fabs(a.y) + fabs(b.x) + fabs(b.y) + fabs(c.x) + fabs(c.y); int i, j, k; u.x = (a.x + b.x + c.x) / 3; u.y = (a.y + b.y + c.y) / 3; while (step > eps) { for (k = 0; k < 10; step /= 2, k ++) { for (i = -1; i <= 1; i ++) { for (j =- 1; j <= 1; j ++) { v.x = u.x + step * i; v.y = u.y + step * j; if (Point::dis(u,a) + Point::dis(u,b) + Point::dis(u,c) > Point::dis(v,a) + Point::dis(v,b) + Point::dis(v,c)) u = v; } } } } return u; } }; int main(void) { return 0; }
相关文章推荐
- 计算几何模板 更新中
- 计算几何模板——不断更新
- 计算几何基础模板(不间断更新)
- 计算几何模板 (更新中)
- 【计算几何】【计算几何模板】【持续更新】
- [置顶] 【计算几何各种小模板总结贴】[不定期更新]
- 计算几何基础模板 以后还会更新
- 【模板整合】【及时更新】【天坑】计算几何模板
- uva 10652 凸包 + 更新版计算几何模板
- 三维计算几何模板[不定期更新]
- UVa 11178计算几何 模板题
- 计算几何初步模板
- 计算几何模板1 点部分
- 计算几何:点积模板
- 计算几何:规范相交模板
- 二维下计算几何程序头 (模板)
- 求凸包的周长(计算几何模板)
- 计算几何的简单模板
- ACM计算几何模板
- 计算几何模板——平面最近点对(附C++源代码) - [转载经典]