求直线上最近的二点(在公垂线上)

直线是用一个点和一个方向向量(已单位化)构成.现在要求二条直线最近的点(即公垂线与两直线的交点)

class Cray//射线类

{

protected:

    CVector_3D m_ptStart;//原点

    CVector_3D m_ptDir;//方向

……

}

如果二个直线平行,则两条直线r1,r2的方向向量的内积必为1或-1,此时可计算S2在r1上的投影(点乘),在r1上计算投影点为Sproject,则要求的两点分别为Sproject和S2.如下图所示

如果二直线不平行,可用数学公式求解:

 

设所求公垂线上的二点分别为M(m0, m1, m2), N(n0, n1, n2).二个起点分别为S1(S1x, S1y, S1z), S2(S2x, S2y, S2z),二个方向向量分别为D1(D1x, D1y,D1z)则根据直线的参数方程法:

(M0- S1x)/ D1x = (M1- S1y)/ D1y = (M2- S1z)/ D1z = K1;

(N0- S2x)/ D2x = (N1- S2y)/ D2y = (N2- S2z)/ D2z = K2;

M(m0, m1, m2)可表示为(D1x* K1+ S1x, D1y* K1+ S1y, D1z* K1+ S1z)

N(n0, n1, n2) 可表示为(D2x* K2+ S2x, D2y* K2+ S2y, D2z* K2+ S2z)

而MN点乘D1结果为0，MN点乘D1结果为0。所以，可得：

-a+ K1-t.K2=0（D1.D1=1，D2.D2=1，可设D1.D2=t, D1 .(S1-S2)=-a）

       b+t.K1-K2=0

解之，可得：

K1=(a+ b* dot)/ (1-t*t)

K2=(b+ a* dot)/ (1-t*t)

即M为(D1x* K1+ S1x, D1y* K1+ S1y, D1z* K1+ S1z),N为(D2x* K2+ S2x, D2y* K2+ S2y, D2z* K2+ S2z)求解完毕。

 

代码如下:

double dot= CVector_3D::Dot( m_ptDir, r.Get_Dir() );//求二个方向向量的内积

    double t2= 1- dot* dot;

    CVector_3D v;

 

    if ( abs( t2 ) < INFINITY_LITTLE)//如果二射线平行

    {

       double dProject= CVector_3D::Dot( m_ptDir, m_ptStart- r.Get_Start()) / dot;//获取投影

 

       pt1= m_ptStart;

       CVector_3D::Add( r.Get_Start(), CVector_3D::Multiply( r.Get_Dir(), dProject, v ), pt2);

 

 

    }

    else

    {

       double a= CVector_3D::Dot( m_ptDir, CVector_3D::Subtract( r.Get_Start(), m_ptStart, v));

       double b= CVector_3D::Dot( r.Get_Dir(), CVector_3D::Subtract( m_ptStart, r.Get_Start(), v));

 

       double k1= ( a+ b* dot)/ t2;//二个直线方向参数

       double k2= ( b+ a* dot)/ t2;

 

       CVector_3D::Add( m_ptStart , CVector_3D::Multiply( m_ptDir, k1, v ), pt1);

       CVector_3D::Add( r.Get_Start() , CVector_3D::Multiply( r.Get_Dir(), k2, v ), pt2);

    }