两异面直线公垂线垂足位置的计算方法

 本文参考下面这篇文章的计算方法，使用C++实现了两异面直线公垂线垂足位置的计算

http://wenku.baidu.com/view/44c4aa07bd64783e09122b91.html?re=view

这篇文章是1990年刊登的，奇怪的是我在搜索过程中发现的这篇文章的第三作者在2002年随便改了一下标题，内容几乎一字不差的发表在另一篇期刊上

http://www.doc88.com/p-3425370084181.html，但是我比较了这两篇文章，发现后者在依葫芦画瓢的过程中抄错了很多字母和符号，所以我们以第一篇文章为准。

在算法实现过程中需要用到大量行列式的计算，最高阶不超过三阶，可以自己编写求解行列式值的程序，在这里我使用的是Eigen开源矩阵库，这个库不需要安装，只要下载下来将路径添加到附加包含目录即可。

#include<iostream>
#include <Eigen/Dense>


void main()
{

    //参考文献《异面直线公垂线的垂足位置》

    float x1, y1, z1, m1, n1, p1;//直线1上一点及其方向向量
    float x2, y2, z2, m2, n2, p2;//直线2上一点及其方向向量

    x1 = 0; y1 = 1, z1 = 0; m1 = 0; n1 = 1; p1 =0;

    x2 = 1; y2 = 1; z2 = 1; m2 = 1; n2 =1; p2 = 0;

    float G[3], H[3];//存储两个垂足的坐标
    Eigen::MatrixXf A1(2, 2), B1(2, 2), C1(2, 2), A2(2, 2), B2(2, 2), C2(2, 2);
    float a1, b1, c1, a2, b2, c2;//存储上述矩阵的行列式的值

    A1 << n1, p1, n2, p2;
    B1 << p1, m1, p2, m2;
    C1 << m1, n1, m2, n2;

    a1 = A1.determinant();
    b1 = B1.determinant();
    c1 = C1.determinant();

    A2 << n2, b1, p2, c1;
    B2 << p2, c1, m2, a1;
    C2 << m2, a1, n2, b1;

    a2 = A2.determinant();
    b2 = B2.determinant();
    c2 = C2.determinant();

    Eigen::MatrixXf W1(3, 3);
    W1 << a1, b1, c1, a2, b2, c2, n1, -m1, 0;

    float w1 = W1.determinant();
    float d1 = a2*(x2 - x1) + b2*(y2 - y1) + c2*(z2 - z1);

    G[0] = x1 - d1*m1 * c1 / w1;
    G[1] = y1 - d1*n1 * c1 / w1;
    G[2] = z1 + d1*(a1*m1 + b1*n1) / w1;

    float a3, b3, c3;
    a3 = n1*(m1*n2 - m2*n1) + p1*(m1*p2 - p1*m2);
    b3 = p1*(n1*p2 - n2*p1) + m1*(n1*m2 - n2*m1);
    c3 = m1*(p1*m2 - p2*m1) + n1*(p1*n2 - p2*n1);

    float w2 = n2*(b1*c3 - b3*c1) + m2*(a1*c3 - a3*c1);
    float d2 = a3*(x1 - x2) + b3*(y1 - y2) + c3*(z1 - z2);

    H[0] = x2 - d2*m2*c1 / w2;
    H[1] = y2 - d2*n2*c1 / w2;
    H[2] = z2 + d2*(a1*m2 + b1*n2) / w2;

    std::cout << "垂足坐标G:" << G[0] << " " << G[1] << " " << G[2] << std::endl;
    std::cout << "垂足坐标H:" << H[0] << " " << H[1] << " " << H[2] << std::endl;

}