@总结 - 14@ 计算几何相关模板

基础模板。

#include <cmath>
#include <cstdio>
#include <algorithm>
using namespace std;

const double PI = acos(-1.0);
const double INF = 1E9;
const double EPS = 1E-9;
const int MAXN = 100000;

int dcmp(double x) {return (fabs(x) <= EPS ? 0 : (x > 0 ? 1 : -1));}

struct point{
    double x, y;
    point() : x(0), y(0) {}
    point(double _x, double _y) : x(_x), y(_y) {}
    
    friend point operator - (point A, point B) {return point(A.x - B.x, A.y - B.y);}
    friend point operator + (point A, point B) {return point(A.x + B.x, A.y + B.y);}
    friend point operator * (point A, double k) {return point(A.x * k, A.y * k);}
    friend point operator / (point A, double k) {return point(A.x / k, A.y / k);}
    friend double operator * (point A, point B) {return A.x * B.x + A.y * B.y;} //点积 
    friend double operator ^ (point A, point B) {return A.x * B.y - A.y * B.x;} //叉积 
    friend bool operator == (point A, point B) {return dcmp(A.x - B.x) == 0 && dcmp(A.y - B.y) == 0;}
    friend bool operator < (point A, point B) {return (A.x == B.x ? A.y < B.y : A.x < B.x);}
    
    friend double length(point A) {return sqrt(A * A);}
    friend double angle(point P) {
        double a = atan2(P.y, P.x);
        return (a < 0 ? a + 2*PI : a);
    }
    friend double angle(point P, point Q) {
        double a = atan2(P ^ Q, P * Q);
        a = (a < 0 ? a + 2*PI : a);
        return (a > PI ? 2*PI - a : a);
    }
    friend double area(point pnt[], int N) {
        double sum = 0;
        for(int i=0;i<N;i++)
            sum += (pnt[i] ^ pnt[(i + 1) % N]);
        return fabs(sum) / 2;
    }
    friend double dist(point A, point B) {
        return length(A - B);
    }
    
    friend point normal(point A) {return point(A.y, -A.x);} //法向量 
    friend point unit(point A) {return A / length(A);}
    friend point rotate(point A, double t) {
        return point(cos(t)*A.x - sin(t)*A.y, cos(t)*A.y + sin(t)*A.x);
    }
    friend point bisector(point A, point B) { //角平分线 
        return A * length(B) + B * length(A);
    }
	
	friend double slope(point a, point b) {
		if( a.x == b.x )
			return a.y < b.y ? INF : -INF;
		else return (a.y - b.y) / (a.x - b.x);
	}
    friend void read(point &p) {scanf("%lf%lf", &p.x, &p.y);}
};

struct line{
    point a, ab; //直线 a + ab*t
    double ang;
	line() : a(), ab(), ang() {}
	line(point _a, point _ab) : a(_a), ab(_ab), ang(atan2(_ab.y, _ab.x)) {}
	point b() {return a + ab;}
    line(double _a, double _b, double _c) {
        ab = point(-_b, _a);
        a = (_a == 0 ? point(0, -_c/_b) : point(-_c/_a, 0));
    }//直线标准式 Ax + By + C = 0 的转化 
	
	friend bool on_left(line L, point P) {
		return dcmp(L.ab ^ (P - L.a)) > 0;
	}
	friend bool operator < (line a, line b) {return a.ang < b.ang;}
    friend bool is_on_line(line L, point P) {
        return dcmp(((P - L.a)^(P - L.b()))) == 0;
    }
    friend bool is_on_ray(line L, point P) {
        return is_on_line(L, P) && (dcmp((P - L.a)*L.ab) >= 0 );
    }
    friend bool is_on_segment(line L, point P) {
        return is_on_ray(L, P) && (dcmp((P - L.b())*L.ab) <= 0);
    }
    friend bool is_same_side(line L, point P1, point P2) {
        return (dcmp(((P1 - L.a)^(P1 - L.b()))*((P2 - L.a)^(P2 - L.b()))) > 0);
    }
    friend bool is_intersect_ray(line L1, line L2) {
        return dcmp(((L2.a - L1.a)^L2.ab)/(L1.ab^L2.ab)) > 0 &&
            dcmp(((L1.a - L2.a)^L1.ab)/(L2.ab^L1.ab)) > 0;
    }
    friend bool is_intersect_segment(line L1, line L2) {
        return (!is_same_side(L1, L2.a, L2.b())) && (!is_same_side(L2, L1.a, L1.b()));
    }
    friend double dist_line(line L, point P) {
        return fabs(L.ab ^ (P - L.a)) / length(L.ab);
    }
    friend double dist_ray(line L, point P) {
        if( dcmp(L.ab * (P - L.a)) <= 0 ) return dist(P, L.a);
        else return dist_line(L, P);
    }
    friend double dist_segment(line L, point P) {
        if( dcmp(L.ab * (P - L.b())) >= 0 ) return dist(P, L.b());
        else return dist_ray(L, P);
    }
    friend double dist_ray(line L1, line L2) {
        if( is_intersect_ray(L1, L2) ) return 0;
        return min(dist_ray(L2, L1.a), dist_ray(L1, L2.a));
    }
    friend double dist_segment(line L1, line L2) {
        if( is_intersect_segment(L1, L2) ) return 0;
        return min(min(dist_segment(L2, L1.a), dist_segment(L2, L1.b())),
            min(dist_segment(L1, L2.a), dist_segment(L1, L2.b())));
    }
    friend point intersect_line(line L1, line L2) {
        return L1.a + L1.ab*((L2.a - L1.a)^L2.ab)/(L1.ab^L2.ab);
    }
    friend point projection(line L, point P) {
        return L.a + L.ab*((P - L.a) * L.ab)/(L.ab * L.ab);
    }//投影
    friend point reflection(line L, point P) {
        return projection(L, P)*2 - P;
    }//对称点
    friend void transform(line L, double &A, double &B, double &C) {
        A = L.ab.y, B = -L.ab.x, C = L.a.y*L.ab.x - L.a.x*L.ab.y;
    }//直线标准式 Ax + By + C = 0 的转化 
    friend void read(line &l) {
        point P1, P2; read(P1), read(P2);
        l = line(P1, P2 - P1);
    }
};


void convex(point *A, int n, point *B, int &m) {
    static point t[MAXN + 5], s[MAXN + 5];
    for(int i=0;i<n;i++) t[i] = A[i]; sort(t, t + n);
    
    int cnt = 0, tp = 0;
    for(int i=0;i<n;i++) {
    	if( i && A[i] == A[i-1] ) continue; // 注意判重点 
        while( tp >= 2 && slope(s[tp - 1], s[tp]) > slope(s[tp], t[i]) )
            tp--;
        s[++tp] = t[i];
    }
    for(int i=1;i<=tp;i++) B[cnt++] = s[i];
    tp = 0;
    for(int i=0;i<n;i++) {
    	if( i && A[i] == A[i-1] ) continue; // 注意判重点
        while( tp >= 2 && slope(s[tp - 1], s[tp]) < slope(s[tp], t[i]) )
            tp--;
        s[++tp] = t[i];
    }
    for(int i=tp-1;i>=2;i--) B[cnt++] = s[i];
    m = cnt;
}
/*
求凸包注意两个常常需要讨论的点:
(1)N 个点共线。
(2)凸包上的边是否允许多点共线(代码中给的是允许共线的版本)。

update in 2020/06/03:写凸包可以不用求斜率,直接判叉积。
void convex(point *A, int n, point *B, int &m) {
    static point t[N + 5]; int k = n;
	for(int i=0;i<k;i++) t[i] = A[i]; sort(t, t + k);
    
    int o = 1; m = 0;
    for(int i=0;i<k;i++) {
        while( m > o && ((B[m - 1] - B[m - 2]) ^ (t[i] - B[m - 1])) >= 0 )
            m--;
        B[m++] = t[i];
    }
    
    o = m;
    for(int i=k-2;i>=0;i--) {
        while( m > o && ((B[m - 1] - B[m - 2]) ^ (t[i] - B[m - 1])) >= 0 )
            m--;
        B[m++] = t[i];
	}
	B[--m] = (point){0, 0};
}
*/

line l[MAXN + 5], ql[MAXN + 5]; point qp[MAXN + 5]; int s, t;
void insert(line x) {
	while( s < t && !on_left(x, qp[t-1]) ) t--;
	while( s < t && !on_left(x, qp[s]) ) s++;
//注意这里必须先弹队尾再弹队首 
	if( dcmp(x.ab ^ ql[t].ab) == 0 ) {
		if( on_left(ql[t], x.a) ) ql[t] = x;
	}
	else ql[++t] = x;
	if( s < t ) qp[t-1] = intersect_line(ql[t-1], ql[t]);
}
bool half_plane(line *A, int n) {
	for(int i=0;i<n;i++) l[i] = A[i];
	l[n++] = line(point(-INF, -INF), point(1, 0));
	l[n++] = line(point(-INF, INF), point(0, -1));
	l[n++] = line(point(INF, INF), point(-1, 0));
	l[n++] = line(point(INF, -INF), point(0, 1));
	
	sort(l, l + n); ql[s = t = 1] = l[0];
	for(int i=1;i<n;i++) insert(l[i]);
	while( s < t && !on_left(ql[s], qp[t-1]) ) t--; //模拟插入队首直线 
	
	return t - s + 1 >= 3;
}//只判断是否存在半平面交,其他用途同理
/*
用半平面交时注意一个细节:直线上的点是否在半平面上。
主要是最后半平面交可能退化成单点,如果按上面代码来写就是无解的
*/
posted @ 2020-01-01 18:29  Tiw_Air_OAO  阅读(66)  评论(0编辑  收藏
把博客园图标替换成自己的图标
把博客园图标替换成自己的图标end