三维几何之求空间俩线段的公垂线,以及分数类

来源:互联网 发布:手机淘宝阿里旺旺在哪 编辑:程序博客网 时间:2024/05/19 19:31

分析:

首先判断线段俩直线是否平行(或重合),如果是的话直接求。考虑4个端点到另外一条线段的距离,取最小值即可。

如果不平行或重合,说明俩条直线是异面直线,这时最短距离既可能是某端点到另外一条线段的距离,也可能是异面直线的最短距离。

如何求异面直线的最短距离?假设俩条直线分别为l1=(p1,v1)和l2=(p2,v2),那么最短距离会在某个q1=p1+sv1和q2=p2+tv2上取到,其中q1和q2分别在l1和l2上,且q1q2是这俩条异面直线的公垂线。

向量q1q1=q2-q1=p2-p1+tv2-sv1垂直于V1,因此Dot(p2-p1+tv2-sv1,v1)=0,根据分配率,有Dot(p2-p1,v1)+t*Dot(v2,v1)-s*Dot(v1,v1)=0.注意这里的三个点积都是可以直接算出来的,因此实际上得到的是一个关于t和s的一次方程。根据q1q2垂直于v2还可以得到一个一次方程,联立求解即可。下面的代码直接使用了经过复杂化简以后的结果。

注意本题比较特殊,要求以分数方式输出。所以可以考虑定义一个Rat类(代表Rational)来保存和计算有理数,并且重载加法、减法、和乘法,然后把本题用到的俩个关键函数:点到线段距离Distace2Tosegment和异面直线的最小距离LineDistace3D改写成返回有理数的版本。

代码:

#include<iostream>#include<cstdio>#include<cstring>#include<cstdlib>#include<cctype>#include<cmath>#include<string>#include<map>#include<set>#include<vector>#include<queue>#include<stack>#define LL long longusing namespace std;const double eps=1e-6;int dcmp(double x){    return fabs(x)<eps?0:(x>0?1:-1);}struct Point3{    LL x,y,z;    Point3(LL x,LL y,LL z):x(x),y(y),z(z){}    Point3() {}    void in()    {        cin>>x>>y>>z;    }};typedef Point3 Vector3;Vector3 operator +(Vector3 A,Vector3 B){    return Vector3(A.x+B.x,A.y+B.y,A.z+B.z);}Vector3 operator -(Vector3 A,Vector3 B){    return Vector3(A.x-B.x,A.y-B.y,A.z-B.z);}Vector3 operator * (Vector3 A,int p){    return Vector3(A.x*p,A.y*p,A.z*p);}Vector3 operator / (Vector3 A,int p){    return Vector3(A.x/p,A.y/p,A.z/p);}bool operator == (Vector3 A,Vector3 B){    return A.x==B.x&&A.y==B.y&&A.z==B.z;}int Dot(Vector3 A,Vector3 B){    return A.x*B.x+A.y*B.y+A.z*B.z;}double Length(Vector3 A){    return sqrt(Dot(A,A));}int Length2(Vector3 A){    return Dot(A,A);}Vector3 Cross(Vector3 A,Vector3 B){    return Vector3(A.y*B.z-A.z*B.y,A.z*B.x-A.x*B.z,A.x*B.y-A.y*B.x);}LL gcd(LL a,LL b){return b?gcd(b,a%b):a;}LL lcm(LL a,LL b){return a/gcd(a,b)*b;}struct Rat{    LL a,b;    Rat(LL a=0):a(a),b(1) {}    Rat(LL x,LL y):a(x),b(y){        if(b<0) a=-a,b=-b;        LL d=gcd(a,b);        if(d<0) d=-d;        a/=d;        b/=d;    }};Rat operator + (const Rat &A,const Rat &B){    LL x=lcm(A.b,B.b);    return Rat(A.a*(x/A.b)+B.a*(x/B.b),x);}Rat operator - (const Rat &A,const Rat &B){    return  A+Rat(-B.a,B.b);}Rat  operator *(const Rat &A,const Rat &B){    return Rat(A.a*B.a,A.b*B.b);}void updatemin(Rat &A,const Rat &B){    if(A.a*B.b>B.a*A.b) A.a=B.a,A.b=B.b;}Rat Rat_Distace2ToSegment(const Point3 &P,const Point3 &A,const Point3 &B){    if(A==B) return Length2(P-A);    Vector3 v1=B-A,v2=P-A,v3=P-B;    if(Dot(v1,v2)<0) return Length2(v2);    else if(Dot(v1,v3)>0) return Length2(v3);    else return Rat(Length2(Cross(v1,v2)),Length2(v1));}bool Rat_LineDistace3D(const Point3 &p1,const Point3 &u,const Point3 &p2,const Vector3 &v, Rat &s){    LL b =(LL)Dot(u,u)*Dot(v,v)-(LL)Dot(u,v)*Dot(u,v);    if(b==0) return false;    LL a=(LL)Dot(u,v)*Dot(v,p1-p2)-(LL)Dot(v,v)*Dot(u,p1-p2);    s=Rat(a,b);    return true;}void Rat_GetPointOnLine(const Point3 &A,const Point3 &B,const Rat &t,Rat &x,Rat &y,Rat &z){    x=Rat(A.x)+Rat(B.x-A.x)*t;    y=Rat(A.y)+Rat(B.y-A.y)*t;    z=Rat(A.z)+Rat(B.z-A.z)*t;}Rat Rat_Distance2(const Rat &x1,const Rat &y1,const Rat &z1,const Rat &x2,const Rat &y2,const Rat &z2){    return (x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)+(z1-z2)*(z1-z2);}int main(){    int T;    cin>>T;    while(T--)    {        Point3 A,B,C,D;        A.in();        B.in();        C.in();        D.in();        Rat s,t;        bool ok=false;        Rat ans=Rat(1000000000);        if(Rat_LineDistace3D(A,B-A,C,D-C,s))            if(s.a>0&&s.a<s.b&&Rat_LineDistace3D(C,D-C,A,B-A,t))                if(t.a>0&&t.a<t.b)        {            ok=true;            Rat x1,y1,z1,x2,y2,z2;            Rat_GetPointOnLine(A,B,s,x1,y1,z1);            Rat_GetPointOnLine(C,D,t,x2,y2,z2);            ans=Rat_Distance2(x1,y1,z1,x2,y2,z2);        }        if(!ok)        {            updatemin(ans,Rat_Distace2ToSegment(A,C,D));            updatemin(ans,Rat_Distace2ToSegment(B,C,D));            updatemin(ans,Rat_Distace2ToSegment(C,A,B));            updatemin(ans,Rat_Distace2ToSegment(D,A,B));        }        printf("%lld %lld\n",ans.a,ans.b);    }    return 0;}

0 0
原创粉丝点击