和项目组研究计算几何
一，凸包
二维凸包：可用极角排序后每次倍增一个点看是否在栈顶两个点内侧，在内侧就加入栈，否则弹出一个。优化是在省去了计算极角的计算量，采用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)); } //三角形面积*2 double area(Point a, Point b, Point c) { return vlen((b - a) * (c - a)); } //四面体a到其他三点有向距离有向体积*6 double 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; }
