给出两直线经过的两点,没有平行的情况,求两直线最短距离和最短距离的两点。
#include <iostream> #include <string.h> #include <stdio.h> #include <math.h> using namespace std; const double eps = 1e-8; //三维空间点 struct Point { double x, y, z; Point(double x=0,double y=0,double z=0): x(x),y(y),z(z){} Point(const Point& a) { x = a.x; y = a.y; z = a.z; return; } void Print() { printf("%lf %lf %lf\n", x, y, z); } Point operator + (Point &t) { return Point(x+t.x, y+t.y, z+t.z); } }; //空间直线 struct Line { Point a,b; }; //空间平面 struct Plane { Point a,b,c; Plane(){} Plane(Point a, Point b, Point c):a(a),b(b),c(c){} void showPlane() { a.Print(); b.Print(); c.Print(); return; } }; double dcmp(double t) { if(fabs(t) < eps) return 0; return t < 0 ? -1 : 1; } //三维叉积 Point cross(Point u,Point v) { Point 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; } //三维点积 double multi(Point u,Point v) { return u.x * v.x + u.y * v.y + u.z * v.z; } //矢量差 Point sub(Point u,Point v) { Point ret; ret.x = u.x - v.x; ret.y = u.y - v.y; ret.z = u.z - v.z; return ret; } //两点距离 double dist(Point p1, Point 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 VectorLength(Point p) { return sqrt(p.x*p.x + p.y*p.y + p.z*p.z); } //空间直线距离 double LineToLine(Line u,Line v,Point &tmp ) { tmp = cross(sub(u.a,u.b),sub(v.a,v.b)); return fabs(multi(sub(u.a,v.a),tmp))/VectorLength(tmp); } //取平面法向量 Point normalVector(Plane s) { return cross(sub(s.a,s.b),sub(s.b,s.c)); } //空间平面与直线的交点 Point Intersection(Line l,Plane s) { Point ret = normalVector(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; } /************以上为模板*************/ void work(Line A, Line B) { Point normal; double d = LineToLine( A, B, normal ); printf("%.6lf\n",d); Plane alpha = Plane(A.a, A.b, A.a + normal); Plane beta = Plane(B.a, B.b, B.a + normal); Point u = Intersection(B,alpha); Point v = Intersection(A,beta); printf("%.6lf %.6lf %.6lf %.6lf %.6lf %.6lf\n", v.x, v.y, v.z, u.x, u.y, u.z ); } int main() { int T; scanf("%d",&T); while(T--) { Line A,B; scanf("%lf%lf%lf", &A.a.x, &A.a.y, &A.a.z); scanf("%lf%lf%lf", &A.b.x, &A.b.y, &A.b.z); scanf("%lf%lf%lf", &B.a.x, &B.a.y, &B.a.z); scanf("%lf%lf%lf", &B.b.x, &B.b.y, &B.b.z); work(A,B); } return 0; }