hdu 4741 Save Labman No.004 计算几何

直接贴模版！！！

求异面直线的距离及坐标！！

代码如下：

#include<iostream>
#include<stdio.h>
#include<algorithm>
#include<iomanip>
#include<cmath>
#include<cstring>
#include<vector>
#define I(p) scanf("%lf%lf%lf",&p.x,&p.y,&p.z)
#define eps 1e-20
#define zero(x) (((x)>0?(x):-(x))<eps)
using namespace std;
struct point
{
      double x,y,z;
      point operator+(point a)
      {
            point c;
            c.x=x+a.x;
            c.y=y+a.y;
            c.z=z+a.z;
            return c;
      }
};
struct line{point a,b;}l1,l2;
struct plane
{
      point a,b,c;
      plane(point _a,point _b,point _c)
      {
            a=_a;
            b=_b;
            c=_c;
      }
};
point xmult(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 dmult(point u,point v)
{
      return u.x*v.x+u.y*v.y+u.z*v.z;
}
point subt(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;
}
//取平面法向量

point  pvec(plane s)
{

    return  xmult(subt(s.a,s.b),subt(s.b,s.c));

}

point  pvec(point  s1,point  s2,point  s3)
{
    return  xmult(subt(s1,s2),subt(s2,s3));

}

double vlen(point p)
{
      return sqrt(p.x*p.x+p.y*p.y+p.z*p.z);
}
double distances(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));
}

//判两直线平行

int parallel(line u,line v)
{

    return vlen(xmult(subt(u.a,u.b),subt(v.a,v.b)))<eps;
}

int parallel(point  u1,point u2,point  v1,point  v2)
{

    return vlen(xmult(subt(u1,u2),subt(v1,v2)))<eps;

}

//判两平面平行

int parallel(plane u,plane v)
{

    return vlen(xmult(pvec(u),pvec(v)))<eps;

}

int parallel(point  u1,point u2,point  u3,point  v1,point  v2,point  v3)
{
    return vlen(xmult(pvec(u1,u2,u3),pvec(v1,v2,v3)))<eps;

}

//判直线与平面平行

int parallel(line l,plane s)
{
    return zero(dmult(subt(l.a,l.b),pvec(s)));

}

int parallel(point  l1,point  l2,point  s1,point  s2,point  s3)
{

    return zero(dmult(subt(l1,l2),pvec(s1,s2,s3)));

}

//判两直线垂直

int perpendicular(line u,line v)
{

    return zero(dmult(subt(u.a,u.b),subt(v.a,v.b)));

}
int perpendicular(point  u1,point u2,point  v1,point  v2)
{

    return zero(dmult(subt(u1,u2),subt(v1,v2)));

}

//判两平面垂直

int perpendicular(plane u,plane v)
{

    return zero(dmult(pvec(u),pvec(v)));

}

int perpendicular(point  u1,point u2,point  u3,point  v1,point  v2,point  v3)
{
    return zero(dmult(pvec(u1,u2,u3),pvec(v1,v2,v3)));

}

//判直线与平面平行

int perpendicular(line l,plane s)
{
    return vlen(xmult(subt(l.a,l.b),pvec(s)))<eps;

}

int perpendicular(point  l1,point  l2,point  s1,point  s2,point  s3)
{

    return vlen(xmult(subt(l1,l2),pvec(s1,s2,s3)))<eps;

}

//计算两直线交点,注意事先判断直线是否共面和平行!
//线段交点请另外判线段相交(同时还是要判断是否平行!)
point  intersection(point  u1,point u2,point  v1,point  v2)
{

    point  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;

}
point  intersection(point  l1,point  l2,point  s1,point  s2,point  s3)
{

    point  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;
}
//计算两平面交线,注意事先判断是否平行,并保证三点不共线!
line intersection(plane u,plane v)
{

    line 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;
}
line intersection(point  u1,point  u2,point u3,point  v1,point  v2,point  v3)
{

    line 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;

}
//计算两直线交点,注意事先判断直线是否共面和平行!
//线段交点请另外判线段相交(同时还是要判断是否平行!)
point  intersection(line u,line v)
{

    point  ret=u.a;
    double x, y;
    x= ((u.a.x-v.a.x)*(v.a.y-v.b.y)-(u.a.y-v.a.y)*(v.a.x-v.b.x));
    y= ((u.a.x-u.b.x)*(v.a.y-v.b.y)-(u.a.y-u.b.y)*(v.a.x-v.b.x));
    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));
    if(fabs(x)>eps){
    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;

}
int main()
{
      int t;
      double dis,len1,len2;
      scanf("%d",&t);
      while(t--){
            I(l1.a);
            I(l1.b);
            I(l2.a);
            I(l2.b);
            point n=xmult(subt(l1.a,l1.b),subt(l2.a,l2.b)),a,b;
            a=l1.a+n;
            b=l2.a+n;
            plane p1(l1.a,l1.b,a),p2(l2.a,l2.b,b);
            line m=intersection(p1,p2);
            point ans1=intersection(m,l1);
            point ans2=intersection(m,l2);
            printf("%.6lf
",distances(ans1,ans2));
            printf("%.6lf %.6lf %.6lf %.6lf %.6lf %.6lf
",ans1.x,ans1.y,ans1.z,ans2.x,ans2.y,ans2.z);
      }
      return 0;
}