三维几何之求空间俩线段的公垂线，以及分数类-白红宇的个人博客

发布日期：2021-10-08 15:48:57 浏览次数：25 分类：技术文章

本文共 3206 字，大约阅读时间需要 10 分钟。

分析：

首先判断线段俩直线是否平行（或重合），如果是的话直接求。考虑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
   
    #include
    
     #include
     
      #include
      
       #include
       
        #include
        
         #include
         
          #include
          #include
           
            #include
            
             #include
             
              #include
              
               #define LL long longusing namespace std;const double eps=1e-6;int dcmp(double x){ return fabs(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
                
                 0&&t.a

转载地址：https://blog.csdn.net/ONE_PIECE_HMH/article/details/45795527 如侵犯您的版权，请留言回复原文章的地址，我们会给您删除此文章，给您带来不便请您谅解！

上一篇：三维几何之判断俩个一个六面题的重心稳定性

下一篇：简单三维几何，判断俩个三角形是否相交

发表评论

关于作者

喝酒易醉，品茶养心，人生如梦，品茶悟道，何以解忧？唯有杜康！

-- 愿君每日到此一游！

发表评论

最新留言

关于作者

推荐文章