一,凸包

二维凸包:可用极角排序后每次倍增一个点看是否在栈顶两个点内侧,在内侧就加入栈,否则弹出一个。优化是在省去了计算极角的计算量,采用xy排序,从x最小倍增一遍求上边界,再从x最大反过来求下边界。二维凸包构造想法比较简单。

凸多边形有一些有趣的性质

设边数为n

内角和 =(n-2)×180°

外角和 = 360°

对角线的条数=C(n,2)-n=n(n-3)/2

欧拉公式 凸多边形有n个点,m条边,r个面,则 n - m + r = 2

1.POJ1113 注意精度,这里有个公式 城堡围墙长度最小值 = 城堡顶点坐标构成的散点集的凸包总边长 + 半径为L的圆周长

#include<iostream>
#include <algorithm>
#include<math.h>
using namespace std;
#define eps 1e-8
const double PI = acos(-1.0);
struct point {double x, y;point() {}point(int a, int b) {x = a, y = b;}
};int sgn(double a) {if (fabs(a) < eps)return 0;if (a > eps)return 1;return -1;
}
double dis(point a, point b) {return sqrt((b.x - a.x) * (b.x - a.x) + (b.y - a.y) * (b.y - a.y));
}double cross(point op, point sp, point ep) {return (sp.x - op.x) * (ep.y - op.y) - (ep.x - op.x) * (sp.y - op.y);
}
int cmp(point a, point b) {return (a.y == b.y) ? (a.x < b.x) : a.y < b.y;
}
double graham(point pnt[], int n, point res[]) {int i, len, top = 1;sort(pnt, pnt + n, cmp);res[0] = pnt[0];res[1] = pnt[1];for (i = 2; i < n; i++) {while (top && sgn(cross(pnt[i], res[top], res[top - 1])) >= 0)top--;res[++top] = pnt[i];}len = top;for (i = n - 2; i >= 0; i--) {while (top != len && sgn(cross(pnt[i], res[top], res[top - 1])) >= 0)top--;res[++top] = pnt[i];}res[top] = res[0];double ans = 0;for (int p = 0; p < top; p++) {ans += dis(res[p], res[p + 1]);}return ans;
}int main() {
//    freopen("data3.txt", "r", stdin);point pnt[50003], res[50003];int i, n;double L;while (scanf("%d%lf", &n, &L) != EOF) {if (n == 2) {for (i = 0; i < 2; i++)scanf("%lf%lf", &pnt[i].x, &pnt[i].y);printf("%.0lf\n", dis(pnt[0], pnt[1]) + L * 2 * PI);continue;}for (i = 0; i < n; i++)scanf("%lf%lf", &pnt[i].x, &pnt[i].y);printf("%.0lf\n", graham(pnt, n, res) + L * 2 * PI);}return 0;
}

2.HDU 1872 枚举状态,使用位运算

#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
#include <iostream>
using namespace std;
#define MAXN 20
int N;
#define eps 1e-8
const double PI = acos(-1.0);struct point {double x, y;point() {}point(double a, double b) {x = a, y = b;}
};
struct Node {point pos;int val, len;
} tree[MAXN];
int sgn(double a) {if (fabs(a) < eps)return 0;if (a > eps)return 1;return -1;
}
double dis(point a, point b) {return sqrt((b.x - a.x) * (b.x - a.x) + (b.y - a.y) * (b.y - a.y));
}double cross(point op, point sp, point ep) {return (sp.x - op.x) * (ep.y - op.y) - (ep.x - op.x) * (sp.y - op.y);
}
int cmp(point a, point b) {if (a.y < b.y)return 1;if (a.y == b.y)if (a.x < b.x)return 1;return 0;
}
double graham(point pnt[], int n) {point res[MAXN];int i, len, top = 1;sort(pnt, pnt + n, cmp);res[0] = pnt[0];res[1] = pnt[1];for (i = 2; i < n; i++) {while (top && sgn(cross(pnt[i], res[top], res[top - 1])) >= 0)top--;res[++top] = pnt[i];}len = top;for (i = n - 2; i >= 0; i--) {while (top != len && sgn(cross(pnt[i], res[top], res[top - 1])) >= 0)top--;res[++top] = pnt[i];}res[top] = res[0];double ans = 0;for (int p = 0; p < top; p++) {ans += dis(res[p], res[p + 1]);}return ans;
}
struct Ans {int bit;double ex_len;Ans() {bit = 0, ex_len = 0;}
};
int main() {
//    freopen("data3.txt", "r", stdin);
    point pnt[MAXN];int tot;int min_val, min_num;int cas = 1;while (scanf("%d", &N) && N) {if (cas > 1)puts("");printf("Forest %d\n", cas++);Ans ans;min_val = min_num = 0x3fffffff;for (int i = 0; i < N; i++) {scanf("%lf%lf%d%d", &tree[i].pos.x, &tree[i].pos.y, &tree[i].val,&tree[i].len);}int cut_len, cut_val;for (int bit = 0; bit < (1 << N); bit++) {cut_len = cut_val = tot = 0;for (int j = 0; j < N; j++) {if (bit & (1 << j)) {cut_len += tree[j].len;cut_val += tree[j].val;} else {pnt[tot].x = tree[j].pos.x, pnt[tot++].y = tree[j].pos.y;}}if (cut_val > min_val)continue;double need_len = graham(pnt, tot);if (need_len <= cut_len) {if (min_val > cut_val|| (min_val == cut_val && N - tot < min_num)) {min_val = cut_val, min_num = N - tot;ans.bit = bit;ans.ex_len = cut_len - need_len;}}}printf("Cut these trees:");for (int i = 0; i < N; i++)if (ans.bit & 1 << i)printf(" %d", i + 1);puts("");printf("Extra wood: %.2lf\n", ans.ex_len);}return 0;
}

三维凸包 倍增思想

#include<iostream>
#include<cmath>
#include<cstring>
#include<algorithm>
using namespace std;const int MAXN = 505;
const double eps = 1e-8;int zero(double x) {if (fabs(x) < eps)return 0;return (x > eps) ? 1 : -1;
}
struct Point {double x, y, z;Point() {}Point(double xx, double yy, double zz) :x(xx), y(yy), z(zz) {}Point operator -(const Point p1) {  //两向量之差return Point(x - p1.x, y - p1.y, z - p1.z);}Point operator *(Point p) {   //叉乘return Point(y * p.z - z * p.y, z * p.x - x * p.z, x * p.y - y * p.x);}double operator ^(Point p) { //点乘return (x * p.x + y * p.y + z * p.z);}
};struct CH3D {struct face {int a, b, c; //凸包一个面上三个点的编号bool ok;  //表示该面是否属于最终凸包中的面
    };int n;   //顶点数Point P[MAXN];  //初始顶点int triangle_num;  //凸包表面的三角形数face F[8 * MAXN];int g[MAXN][MAXN]; //凸包表面的三角形double vlen(Point a) {return sqrt(a.x * a.x + a.y * a.y + a.z * a.z);}Point cross(const Point &a, const Point &b, const Point &c) {return Point((b.y - a.y) * (c.z - a.z) - (b.z - a.z) * (c.y - a.y),-((b.x - a.x) * (c.z - a.z) - (b.z - a.z) * (c.x - a.x)),(b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x));}//三角形面积*2double area(Point a, Point b, Point c) {return vlen((b - a) * (c - a));}//四面体a到其他三点有向距离有向体积*6double volume(Point a, Point b, Point c, Point d) {return (b - a) * (c - a) ^ (d - a);}//点在面同向返回正数int dblcmp(Point &p, face &f) {double v = volume(P[f.a], P[f.b], P[f.c], p);if (fabs(v) < eps)return 0;return (v > eps) ? 1 : -1;}void deal(int p, int a, int b) {int f = g[a][b];face add;if (F[f].ok) {if (dblcmp(P[p], F[f]) > 0)dfs(p, f);else {add.a = b;add.b = a;add.c = p;add.ok = 1;g[p][b] = g[a][p] = g[b][a] = triangle_num;F[triangle_num++] = add;}}}void dfs(int p, int now) {F[now].ok = 0;deal(p, F[now].b, F[now].a);deal(p, F[now].c, F[now].b);deal(p, F[now].a, F[now].c);}//调整使前四个点不共面bool check_4point() {int i, flag = 1;for (i = 1; i < n; i++) {if (vlen(P[0] - P[i]) > eps) {swap(P[1], P[i]);flag = 0;break;}}if (flag)return 0;flag = 1;for (i = 2; i < n; i++) {if (vlen((P[0] - P[1]) * (P[1] - P[i])) > eps) {swap(P[2], P[i]);flag = 0;break;}}if (flag)return 0;flag = 1;for (i = 3; i < n; i++) {if (fabs((P[0] - P[1]) * (P[1] - P[2]) ^ (P[0] - P[i])) > eps) {swap(P[3], P[i]);return 1;}}return 0;}//构建三维凸包void build() {int i, j, tmp;face add;triangle_num = 0;if (n < 4)return;if (!check_4point()) //调整前4个点,找不到不共面4个点则无法构建凸包; 若以保证,则可去掉return;for (i = 0; i < 4; i++) { //根据前4个点产生4个三角面add.a = (i + 1) % 4;add.b = (i + 2) % 4;add.c = (i + 3) % 4;add.ok = 1;if (dblcmp(P[i], add) > 0)swap(add.b, add.c);g[add.a][add.b] = g[add.b][add.c] = g[add.c][add.a] = triangle_num;F[triangle_num++] = add;}for (i = 4; i < n; i++) {for (j = 0; j < triangle_num; j++) {if (F[j].ok && dblcmp(P[i], F[j]) > 0) {dfs(i, j);break;}}}tmp = triangle_num;for (i = triangle_num = 0; i < tmp; i++)if (F[i].ok) {F[triangle_num++] = F[i];}}//三维凸包表面积double area() {double res = 0.0;if (n == 3) {Point p = cross(P[0], P[1], P[2]);res = vlen(p) / 2.0;return res;}for (int i = 0; i < triangle_num; i++)res += fabs(area(P[F[i].a], P[F[i].b], P[F[i].c]));return res / 2.0;}//三维凸包体积double volume() {double res = 0.0;Point tmp(0, 0, 0);for (int i = 0; i < triangle_num; i++)res += volume(tmp, P[F[i].a], P[F[i].b], P[F[i].c]);return fabs(res / 6.0);}//判断两个三角形共面bool same(int s, int t) {Point a = P[F[s].a];Point b = P[F[s].b];Point c = P[F[s].c];return !zero(volume(a, b, c, P[F[t].a]))&& !zero(volume(a, b, c, P[F[t].b]))&& !zero(volume(a, b, c, P[F[t].c]));}//三维凸包表面多边形个数int polygon() {int i, j, res, flag;for (i = res = 0; i < triangle_num; i++) {flag = 1;for (j = 0; j < i; j++)if (same(i, j)) {flag = 0;break;}res += flag;}return res;}
};CH3D hull;int main() {
//    freopen("data2.txt", "r", stdin);while (~scanf("%d", &hull.n)) {for (int i = 0; i < hull.n; i++)scanf("%lf%lf%lf", &hull.P[i].x, &hull.P[i].y, &hull.P[i].z);hull.build();printf("%d\n", hull.polygon());
//        res = hull.area();
//        printf("%.3lf\n", res);
    }return 0;
}

二,半平面

对所有线极角排序去重,之后维护一个直线双向栈,倍增一条线时,检查双向栈两端栈顶两条线的交点是否在新线内侧,内侧就加入栈,否则弹出一条线,最后双向栈中的各条线组成图形的核,可以在这个区域里看到原图形的每一个位置。用这个求内切的最大的圆,就是把各边按二分半径收缩后看是否还存在核的最大半径即解。

#include <iostream>
#include <algorithm>
#include <cmath>
using namespace std;const double eps = 1e-8;
const int maxn = 105;int deq[maxn], tail, head, order[maxn], ln ,pn;struct Point {double x, y;
} p[maxn];struct Line {Point a, b;double angle;
} l[maxn];int sgn(double a){if(fabs(a)<eps) return 0;if(a>eps) return 1;return -1;
}int cross(Point op, Point sp, Point ep)
{return (sp.x-op.x) * (ep.y-op.y) - (ep.x-op.x) * (sp.y-op.y);
}bool cmp(int u, int v) {int d = sgn(l[u].angle-l[v].angle);if (!d) return sgn(cross(l[u].a, l[v].a, l[v].b)) < 0; //极角相同时,只保留最靠里面的那条,将更靠半平面里面的放在前面去重时删除return d < 0;
}void getIntersect(Line l1, Line l2, Point& p) { //求两线交点double dot1,dot2;dot1 = cross(l2.a, l1.b, l1.a);dot2 = cross(l1.b, l2.b, l1.a);p.x = (l2.a.x * dot2 + l2.b.x * dot1) / (dot2 + dot1);p.y = (l2.a.y * dot2 + l2.b.y * dot1) / (dot2 + dot1);
}bool judge(Line l0, Line l1, Line l2) {Point p;getIntersect(l1, l2, p);return sgn(cross(p, l0.a, l0.b)) > 0; //大于小于符号与上面cmp()中注释处相反
}bool halfPlaneIntersection() {int i, j;for (i = 0; i < ln; i++) order[i] = i;sort(order, order+ln, cmp);for (i = 1, j = 0; i < ln; i++)if (sgn(l[order[i]].angle-l[order[j]].angle) > 0)order[++j] = order[i];ln = j + 1;deq[0] = order[0];deq[1] = order[1];head = 0;tail = 1;for (i = 2; i < ln; i++) {while (head < tail && judge(l[order[i]], l[deq[tail-1]], l[deq[tail]])) tail--;while (head < tail && judge(l[order[i]], l[deq[head+1]], l[deq[head]])) head++;deq[++tail] = order[i];}while (head < tail && judge(l[deq[head]], l[deq[tail-1]], l[deq[tail]])) tail--;while (head < tail && judge(l[deq[tail]], l[deq[head+1]], l[deq[head]])) head++;if (head + 1 >= tail) return 0;//求出核多边形点p[]deq[++tail] = deq[head];for(pn=0,i=head;i<tail;i++,pn++)getIntersect(l[deq[i]],l[deq[i+1]],p[pn]);return 1;
}
int main(){//freopen("data.txt","r",stdin);int N,i;int cnt=1;while(scanf("%d",&N)&&N){if(cnt>1) printf("\n");printf("Floor #%d\n",cnt++);for(i=0;i<N;i++)scanf("%lf%lf",&p[i].x,&p[i].y);for(i=1;i<N;i++){l[i-1].a=p[i-1],l[i-1].b=p[i];l[i-1].angle = atan2(p[i].y-p[i-1].y, p[i].x-p[i-1].x); //atan2(dy, dx)返回弧度取值范围为-PI到PI
        }l[N-1].a=p[N-1],l[N-1].b=p[0];l[N-1].angle = atan2(p[0].y-p[N-1].y, p[0].x-p[N-1].x);ln=N;if(halfPlaneIntersection()) printf("Surveillance is possible.\n");else printf("Surveillance is impossible.\n");}
return 0;
}

 

三,三维几何 线段

与二维相似综合使用点乘叉乘计算,三维的叉乘比较特殊,两个向量的叉乘为垂直这两个向量的向量,根据这个性质适合求法线,四面体的高等等。和二维一样,点乘判垂直,叉乘判平行,点乘另一个应用是求一个向量在另一个单位向量上的投影长度。

HDU 4741

//三维几何函数库
#include <math.h>
#define eps 1e-8
#define zero(x) (((x)>0?(x):-(x))<eps)
struct point3{double x,y,z;};
struct line3{point3 a,b;};
struct plane3{point3 a,b,c;};//计算叉积 U x V
point3 xmult(point3 u,point3 v){point3 ret;ret.x=u.y*v.z-v.y*u.z;ret.y=u.z*v.x-u.x*v.z;ret.z=u.x*v.y-u.y*v.x;return ret;
}//计算点积 U . V
double dmult(point3 u,point3 v){return u.x*v.x+u.y*v.y+u.z*v.z;
}//矢量差 U - V
point3 subt(point3 u,point3 v){point3 ret;ret.x=u.x-v.x;ret.y=u.y-v.y;ret.z=u.z-v.z;return ret;
}
//矢量和 U + V
point3 plusv(point3 u, point3 v) {point3 ret;ret.x = u.x + v.x;ret.y = u.y + v.y;ret.z = u.z + v.z;return ret;
}
//取平面法向量
point3 pvec(plane3 s){return xmult(subt(s.a,s.b),subt(s.b,s.c));
}
point3 pvec(point3 s1,point3 s2,point3 s3){return xmult(subt(s1,s2),subt(s2,s3));
}//两点距离,单参数取向量大小
double distance(point3 p1,point3 p2){return sqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y)+(p1.z-p2.z)*(p1.z-p2.z));
}//向量大小
double vlen(point3 p){return sqrt(p.x*p.x+p.y*p.y+p.z*p.z);
}//判三点共线
int dots_inline(point3 p1,point3 p2,point3 p3){return vlen(xmult(subt(p1,p2),subt(p2,p3)))<eps;
}//判四点共面
int dots_onplane(point3 a,point3 b,point3 c,point3 d){return zero(dmult(pvec(a,b,c),subt(d,a)));
}//判点是否在线段上,包括端点和共线
int dot_online_in(point3 p,line3 l){return zero(vlen(xmult(subt(p,l.a),subt(p,l.b))))&&(l.a.x-p.x)*(l.b.x-p.x)<eps&&(l.a.y-p.y)*(l.b.y-p.y)<eps&&(l.a.z-p.z)*(l.b.z-p.z)<eps;
}
int dot_online_in(point3 p,point3 l1,point3 l2){return zero(vlen(xmult(subt(p,l1),subt(p,l2))))&&(l1.x-p.x)*(l2.x-p.x)<eps&&(l1.y-p.y)*(l2.y-p.y)<eps&&(l1.z-p.z)*(l2.z-p.z)<eps;
}//判点是否在线段上,不包括端点
int dot_online_ex(point3 p,line3 l){return dot_online_in(p,l)&&(!zero(p.x-l.a.x)||!zero(p.y-l.a.y)||!zero(p.z-l.a.z))&&(!zero(p.x-l.b.x)||!zero(p.y-l.b.y)||!zero(p.z-l.b.z));
}
int dot_online_ex(point3 p,point3 l1,point3 l2){return dot_online_in(p,l1,l2)&&(!zero(p.x-l1.x)||!zero(p.y-l1.y)||!zero(p.z-l1.z))&&(!zero(p.x-l2.x)||!zero(p.y-l2.y)||!zero(p.z-l2.z));
}//判点是否在空间三角形上,包括边界,三点共线无意义
int dot_inplane_in(point3 p,plane3 s){return zero(vlen(xmult(subt(s.a,s.b),subt(s.a,s.c)))-vlen(xmult(subt(p,s.a),subt(p,s.b)))-vlen(xmult(subt(p,s.b),subt(p,s.c)))-vlen(xmult(subt(p,s.c),subt(p,s.a))));
}
int dot_inplane_in(point3 p,point3 s1,point3 s2,point3 s3){return zero(vlen(xmult(subt(s1,s2),subt(s1,s3)))-vlen(xmult(subt(p,s1),subt(p,s2)))-vlen(xmult(subt(p,s2),subt(p,s3)))-vlen(xmult(subt(p,s3),subt(p,s1))));
}//判点是否在空间三角形上,不包括边界,三点共线无意义
int dot_inplane_ex(point3 p,plane3 s){return dot_inplane_in(p,s)&&vlen(xmult(subt(p,s.a),subt(p,s.b)))>eps&&vlen(xmult(subt(p,s.b),subt(p,s.c)))>eps&&vlen(xmult(subt(p,s.c),subt(p,s.a)))>eps;
}
int dot_inplane_ex(point3 p,point3 s1,point3 s2,point3 s3){return dot_inplane_in(p,s1,s2,s3)&&vlen(xmult(subt(p,s1),subt(p,s2)))>eps&&vlen(xmult(subt(p,s2),subt(p,s3)))>eps&&vlen(xmult(subt(p,s3),subt(p,s1)))>eps;
}//判两点在线段同侧,点在线段上返回0,不共面无意义
int same_side(point3 p1,point3 p2,line3 l){return dmult(xmult(subt(l.a,l.b),subt(p1,l.b)),xmult(subt(l.a,l.b),subt(p2,l.b)))>eps;
}
int same_side(point3 p1,point3 p2,point3 l1,point3 l2){return dmult(xmult(subt(l1,l2),subt(p1,l2)),xmult(subt(l1,l2),subt(p2,l2)))>eps;
}//判两点在线段异侧,点在线段上返回0,不共面无意义
int opposite_side(point3 p1,point3 p2,line3 l){return dmult(xmult(subt(l.a,l.b),subt(p1,l.b)),xmult(subt(l.a,l.b),subt(p2,l.b)))<-eps;
}
int opposite_side(point3 p1,point3 p2,point3 l1,point3 l2){return dmult(xmult(subt(l1,l2),subt(p1,l2)),xmult(subt(l1,l2),subt(p2,l2)))<-eps;
}//判两点在平面同侧,点在平面上返回0
int same_side(point3 p1,point3 p2,plane3 s){return dmult(pvec(s),subt(p1,s.a))*dmult(pvec(s),subt(p2,s.a))>eps;
}
int same_side(point3 p1,point3 p2,point3 s1,point3 s2,point3 s3){return dmult(pvec(s1,s2,s3),subt(p1,s1))*dmult(pvec(s1,s2,s3),subt(p2,s1))>eps;
}//判两点在平面异侧,点在平面上返回0
int opposite_side(point3 p1,point3 p2,plane3 s){return dmult(pvec(s),subt(p1,s.a))*dmult(pvec(s),subt(p2,s.a))<-eps;
}
int opposite_side(point3 p1,point3 p2,point3 s1,point3 s2,point3 s3){return dmult(pvec(s1,s2,s3),subt(p1,s1))*dmult(pvec(s1,s2,s3),subt(p2,s1))<-eps;
}//判两直线平行
int parallel(line3 u,line3 v){return vlen(xmult(subt(u.a,u.b),subt(v.a,v.b)))<eps;
}
int parallel(point3 u1,point3 u2,point3 v1,point3 v2){return vlen(xmult(subt(u1,u2),subt(v1,v2)))<eps;
}//判两平面平行
int parallel(plane3 u,plane3 v){return vlen(xmult(pvec(u),pvec(v)))<eps;
}
int parallel(point3 u1,point3 u2,point3 u3,point3 v1,point3 v2,point3 v3){return vlen(xmult(pvec(u1,u2,u3),pvec(v1,v2,v3)))<eps;
}//判直线与平面平行
int parallel(line3 l,plane3 s){return zero(dmult(subt(l.a,l.b),pvec(s)));
}
int parallel(point3 l1,point3 l2,point3 s1,point3 s2,point3 s3){return zero(dmult(subt(l1,l2),pvec(s1,s2,s3)));
}//判两直线垂直
int perpendicular(line3 u,line3 v){return zero(dmult(subt(u.a,u.b),subt(v.a,v.b)));
}
int perpendicular(point3 u1,point3 u2,point3 v1,point3 v2){return zero(dmult(subt(u1,u2),subt(v1,v2)));
}//判两平面垂直
int perpendicular(plane3 u,plane3 v){return zero(dmult(pvec(u),pvec(v)));
}
int perpendicular(point3 u1,point3 u2,point3 u3,point3 v1,point3 v2,point3 v3){return zero(dmult(pvec(u1,u2,u3),pvec(v1,v2,v3)));
}//判直线与平面平行
int perpendicular(line3 l,plane3 s){return vlen(xmult(subt(l.a,l.b),pvec(s)))<eps;
}
int perpendicular(point3 l1,point3 l2,point3 s1,point3 s2,point3 s3){return vlen(xmult(subt(l1,l2),pvec(s1,s2,s3)))<eps;
}//判两线段相交,包括端点和部分重合
int intersect_in(line3 u,line3 v){if (!dots_onplane(u.a,u.b,v.a,v.b))return 0;if (!dots_inline(u.a,u.b,v.a)||!dots_inline(u.a,u.b,v.b))return !same_side(u.a,u.b,v)&&!same_side(v.a,v.b,u);return dot_online_in(u.a,v)||dot_online_in(u.b,v)||dot_online_in(v.a,u)||dot_online_in(v.b,u);
}
int intersect_in(point3 u1,point3 u2,point3 v1,point3 v2){if (!dots_onplane(u1,u2,v1,v2))return 0;if (!dots_inline(u1,u2,v1)||!dots_inline(u1,u2,v2))return !same_side(u1,u2,v1,v2)&&!same_side(v1,v2,u1,u2);return dot_online_in(u1,v1,v2)||dot_online_in(u2,v1,v2)||dot_online_in(v1,u1,u2)||dot_online_in(v2,u1,u2);
}//判两线段相交,不包括端点和部分重合
int intersect_ex(line3 u,line3 v){return dots_onplane(u.a,u.b,v.a,v.b)&&opposite_side(u.a,u.b,v)&&opposite_side(v.a,v.b,u);
}
int intersect_ex(point3 u1,point3 u2,point3 v1,point3 v2){return dots_onplane(u1,u2,v1,v2)&&opposite_side(u1,u2,v1,v2)&&opposite_side(v1,v2,u1,u2);
}//判线段与空间三角形相交,包括交于边界和(部分)包含
int intersect_in(line3 l,plane3 s){return !same_side(l.a,l.b,s)&&!same_side(s.a,s.b,l.a,l.b,s.c)&&!same_side(s.b,s.c,l.a,l.b,s.a)&&!same_side(s.c,s.a,l.a,l.b,s.b);
}
int intersect_in(point3 l1,point3 l2,point3 s1,point3 s2,point3 s3){return !same_side(l1,l2,s1,s2,s3)&&!same_side(s1,s2,l1,l2,s3)&&!same_side(s2,s3,l1,l2,s1)&&!same_side(s3,s1,l1,l2,s2);
}//判线段与空间三角形相交,不包括交于边界和(部分)包含
int intersect_ex(line3 l,plane3 s){return opposite_side(l.a,l.b,s)&&opposite_side(s.a,s.b,l.a,l.b,s.c)&&opposite_side(s.b,s.c,l.a,l.b,s.a)&&opposite_side(s.c,s.a,l.a,l.b,s.b);
}
int intersect_ex(point3 l1,point3 l2,point3 s1,point3 s2,point3 s3){return opposite_side(l1,l2,s1,s2,s3)&&opposite_side(s1,s2,l1,l2,s3)&&opposite_side(s2,s3,l1,l2,s1)&&opposite_side(s3,s1,l1,l2,s2);
}//计算两直线交点,注意事先判断直线是否共面和平行!
//线段交点请另外判线段相交(同时还是要判断是否平行!)
point3 intersection(line3 u,line3 v){point3 ret=u.a;double t=((u.a.x-v.a.x)*(v.a.y-v.b.y)-(u.a.y-v.a.y)*(v.a.x-v.b.x))/((u.a.x-u.b.x)*(v.a.y-v.b.y)-(u.a.y-u.b.y)*(v.a.x-v.b.x));ret.x+=(u.b.x-u.a.x)*t;ret.y+=(u.b.y-u.a.y)*t;ret.z+=(u.b.z-u.a.z)*t;return ret;
}
point3 intersection(point3 u1,point3 u2,point3 v1,point3 v2){point3 ret=u1;double t=((u1.x-v1.x)*(v1.y-v2.y)-(u1.y-v1.y)*(v1.x-v2.x))/((u1.x-u2.x)*(v1.y-v2.y)-(u1.y-u2.y)*(v1.x-v2.x));ret.x+=(u2.x-u1.x)*t;ret.y+=(u2.y-u1.y)*t;ret.z+=(u2.z-u1.z)*t;return ret;
}//计算直线与平面交点,注意事先判断是否平行,并保证三点不共线!
//线段和空间三角形交点请另外判断
point3 intersection(line3 l,plane3 s){point3 ret=pvec(s);double t=(ret.x*(s.a.x-l.a.x)+ret.y*(s.a.y-l.a.y)+ret.z*(s.a.z-l.a.z))/(ret.x*(l.b.x-l.a.x)+ret.y*(l.b.y-l.a.y)+ret.z*(l.b.z-l.a.z));ret.x=l.a.x+(l.b.x-l.a.x)*t;ret.y=l.a.y+(l.b.y-l.a.y)*t;ret.z=l.a.z+(l.b.z-l.a.z)*t;return ret;
}
point3 intersection(point3 l1,point3 l2,point3 s1,point3 s2,point3 s3){point3 ret=pvec(s1,s2,s3);double t=(ret.x*(s1.x-l1.x)+ret.y*(s1.y-l1.y)+ret.z*(s1.z-l1.z))/(ret.x*(l2.x-l1.x)+ret.y*(l2.y-l1.y)+ret.z*(l2.z-l1.z));ret.x=l1.x+(l2.x-l1.x)*t;ret.y=l1.y+(l2.y-l1.y)*t;ret.z=l1.z+(l2.z-l1.z)*t;return ret;
}//计算两平面交线,注意事先判断是否平行,并保证三点不共线!
line3 intersection(plane3 u,plane3 v){line3 ret;ret.a=parallel(v.a,v.b,u.a,u.b,u.c)?intersection(v.b,v.c,u.a,u.b,u.c):intersection(v.a,v.b,u.a,u.b,u.c);ret.b=parallel(v.c,v.a,u.a,u.b,u.c)?intersection(v.b,v.c,u.a,u.b,u.c):intersection(v.c,v.a,u.a,u.b,u.c);return ret;
}
line3 intersection(point3 u1,point3 u2,point3 u3,point3 v1,point3 v2,point3 v3){line3 ret;ret.a=parallel(v1,v2,u1,u2,u3)?intersection(v2,v3,u1,u2,u3):intersection(v1,v2,u1,u2,u3);ret.b=parallel(v3,v1,u1,u2,u3)?intersection(v2,v3,u1,u2,u3):intersection(v3,v1,u1,u2,u3);return ret;
}//点到直线距离
double ptoline(point3 p,line3 l){return vlen(xmult(subt(p,l.a),subt(l.b,l.a)))/distance(l.a,l.b);
}
double ptoline(point3 p,point3 l1,point3 l2){return vlen(xmult(subt(p,l1),subt(l2,l1)))/distance(l1,l2);
}//点到平面距离
double ptoplane(point3 p,plane3 s){return fabs(dmult(pvec(s),subt(p,s.a)))/vlen(pvec(s));
}
double ptoplane(point3 p,point3 s1,point3 s2,point3 s3){return fabs(dmult(pvec(s1,s2,s3),subt(p,s1)))/vlen(pvec(s1,s2,s3));
}//直线到直线距离
double linetoline(line3 u,line3 v){point3 n=xmult(subt(u.a,u.b),subt(v.a,v.b));return fabs(dmult(subt(u.a,v.a),n))/vlen(n);
}
double linetoline(point3 u1,point3 u2,point3 v1,point3 v2){point3 n=xmult(subt(u1,u2),subt(v1,v2));return fabs(dmult(subt(u1,v1),n))/vlen(n);
}//两直线夹角cos值
double angle_cos(line3 u,line3 v){return dmult(subt(u.a,u.b),subt(v.a,v.b))/vlen(subt(u.a,u.b))/vlen(subt(v.a,v.b));
}
double angle_cos(point3 u1,point3 u2,point3 v1,point3 v2){return dmult(subt(u1,u2),subt(v1,v2))/vlen(subt(u1,u2))/vlen(subt(v1,v2));
}//两平面夹角cos值
double angle_cos(plane3 u,plane3 v){return dmult(pvec(u),pvec(v))/vlen(pvec(u))/vlen(pvec(v));
}
double angle_cos(point3 u1,point3 u2,point3 u3,point3 v1,point3 v2,point3 v3){return dmult(pvec(u1,u2,u3),pvec(v1,v2,v3))/vlen(pvec(u1,u2,u3))/vlen(pvec(v1,v2,v3));
}//直线平面夹角sin值
double angle_sin(line3 l,plane3 s){return dmult(subt(l.a,l.b),pvec(s))/vlen(subt(l.a,l.b))/vlen(pvec(s));
}
double angle_sin(point3 l1,point3 l2,point3 s1,point3 s2,point3 s3){return dmult(subt(l1,l2),pvec(s1,s2,s3))/vlen(subt(l1,l2))/vlen(pvec(s1,s2,s3));
}int main() {
//    freopen("data3.txt", "r", stdin);int T;cin >> T;while (T--) {point3 a[2], b[2];scanf("%lf%lf%lf", &a[0].x, &a[0].y, &a[0].z);scanf("%lf%lf%lf", &a[1].x, &a[1].y, &a[1].z);scanf("%lf%lf%lf", &b[0].x, &b[0].y, &b[0].z);scanf("%lf%lf%lf", &b[1].x, &b[1].y, &b[1].z);point3 dir = xmult(subt(a[0], a[1]), subt(b[0], b[1])); //叉积与两直线垂直的向量double dis = 0;if (vlen(dir) > eps) { //判断共面double d = 1 / vlen(dir);dir.x = dir.x * d;dir.y = dir.y * d;dir.z = dir.z * d;dis = dmult(subt(a[0], b[0]), dir); //与单位向量的点积为这个方向的长度dir.x = dir.x * dis;dir.y = dir.y * dis;dir.z = dir.z * dis;}printf("%.6f\n", fabs(dis));point3 pa = intersection(a[0], a[1], plusv(b[0], dir),plusv(b[1], dir));point3 pb = subt(pa, dir);printf("%lf %lf %lf %lf %lf %lf\n", pa.x, pa.y, pa.z, pb.x, pb.y, pb.z);}return 0;
}

转载于:https://www.cnblogs.com/updateofsimon/p/3423730.html

和项目组研究计算几何相关推荐

  1. [转]清华梦的粉碎—写给清华大学的退学申请

        原文地址:http://74.207.235.141/discussion/4652/      今天无意中看到的转来大家看看~~王垠 四川大学97级本科毕业,保送到清华大学计算机系直博.期间 ...

  2. 国际物流杰信项目总结与面试

    国际物流杰信项目面试总结 1.面试时如何讲解项目? 讲出三个层次, 1)讲项目的背景,讲特色的业务 2)讲业务的复杂度 3)从业务角度牵扯出技术亮点 每一层都要挖陷阱,让面试官问问题. 从面试题中找出 ...

  3. 清华梦的粉碎——转自王垠

    小时候,妈妈给我一个梦.她指着一个大哥哥的照片对我说,这是爸爸的学生,他考上了清华大学,他是我们中学的骄傲.长大后,你也要进入清华大学读书,为我们家争光.我不知道清华是什么样子,但是我知道爱迪生和牛顿 ...

  4. 【转】清华梦的粉碎 - 写给清华大学的退学申请

    清华梦的粉碎-写给清华大学的退学申请(转自王垠Blog) 清华梦的诞生 小时候,妈妈给我一个梦.她指着一个大哥哥的照片对我说,这是爸爸的学生,他考上了清华大学,他是我们中学的骄傲.长大后,你也要进入清 ...

  5. GIS概念介绍和对webgis的理解

    1.什么是GIS 地理信息系统(Geographical Information System,GIS),它是一种计算机系统,具有对空间数据与属性数据进行输入.管理.查询和分析及输出等功能: 地理信息 ...

  6. 适用于实验室的新型能量回收污水处理铜板蚀刻机设计

    适用于实验室的新型能量回收污水处理铜板蚀刻机设计 设计者:范裕莹 董慧鑫 李纳瑞 陆婧 李亮 彭岩 罗豪 指导教师:阳晓宇 链接地址:https://poetryscience.ink/2018/09 ...

  7. 清华梦的粉碎—写给清华大学的退学申请 2005.9.22

    侵删 清华梦的诞生 小时候,妈妈给我一个梦.她指着一个大哥哥的照片对我说,这是爸爸的学生,他考上了清华 大学,他是我们中学的骄傲.长大后,你也要进入清华大学读书,为我们家争光.我不知道清华是什么样子, ...

  8. TDCS刺激强度对健康受试者工作记忆的影响

    经颅直流电刺激(tDCS)可以改善健康受试者工作记忆(WM)的表现.然而,不同文献得到的结论并不一致,也有认为没有效果,且对这些结果的解释是混乱的,包括tDCS强度(电流强度)和伪刺激的差异. 目的: ...

  9. [转]清华梦的粉碎——写给清华大学的退学申请

    [转]清华梦的粉碎--写给清华大学的退学申请 读了全文,感同身受,全文转载. By 王垠(2005.09.22) 作者王垠,非常有思想的一个人,川大计算机系97级本科,2001年毕业后直博保送清华大学 ...

最新文章

  1. 【学习笔记】block、inline(替换元素、不可替换元素)、inline-block的理解
  2. 2018-12-04-Python全栈开发-day92-自动登录
  3. bzoj2006 NOI2010 数据结构+堆维护区间和最大
  4. 细说angular Form addControl方法
  5. memcached的最佳实践方案(转)
  6. C语言实现录入学生信息并按分数排序输出
  7. javaioIOException - Cannot run program javac error 2 No such file or direct
  8. 地图 c-suite_C-Suite的模型
  9. 【51单片机快速入门指南】4.3: I2C读取MPU6050陀螺仪的原始数据
  10. websocket + node.js聊天系统
  11. WebFlux响应式编程基础之 4 reactive stream 响应式流
  12. UNIX环境高级编程——线程
  13. unity byte[]的压缩和解压
  14. Spring Cloud Eureka详解
  15. Linear Regression 和 Logistic Regression的不同(对比)
  16. 简单开发一个java 插件式demo
  17. 安装多个 PHP 版本(PHP7, PHP5)
  18. python做鼠标自动移动_Python实现鼠标自动在屏幕上随机移动功能
  19. 怎样实习才能成长最快
  20. 正常查看网页中压缩的js代码

热门文章

  1. python服务器环境搭建(2)——安装相关软件
  2. Eclipse的详细安装步骤
  3. 学生宿舍管理系统--需求说明、概要设计、详细设计
  4. 进程间通信(IPC)之内存映射mmap和共享内存shm
  5. matlab将x排序 y随之变化,在MATLAB中:XData和YData如何用更改的行数更新?
  6. kafka拉取mysql数据库_kafka里信息用flink获取后放入mysql
  7. Java基础与数据库对应数据--Java基础2阶段
  8. display:none的表单也会被提交
  9. qs.js 更好的处理url参数
  10. php preg_match_all匹配正则,字符串过长时出错