异面直线求最短距离及坐标点

        先将直线方程化为对称式，分别得到两直线方向向量AB=（x1,y1,z1), CD=(x2,y2,z2)，
再将两向量AB, CD叉乘得到其公垂向量N=（x,y,z）,在AB, CD两直线上分别选取点E,F(任意)，得到向量M，求向量M在向量N方向的投影即为两异面直线间的距离了（就是最短距离啦）。

     最短距离的求法：d=|向量N向量M|/|向量N|（上面是两向量的数量积，下面是取模）。*
设两直线与距离的交点分别为S,T，可带入公垂线N的对称式中得到第一个方程，又因为S,T两点分别满足直线AB和CD的方程，所以得到关于S（或T）的第二个方程，联立两个方程分别解出来即可！

 public class Vector
    {
        public double x { get; set; }
        public double y { get; set; }
        public double z { get; set; }
        public Vector(double xx, double yy, double zz)
        {
            x = xx;
            y = yy;
            z = zz;
        }

        private static double Dot(Vector a, Vector b)
        {
            return a.x * b.x + a.y * b.y + a.z * b.z;
        }

        /// <summary>
        /// 模长
        /// </summary>
        /// <param name="a"></param>
        /// <returns></returns>
        private static double GetLength(Vector a)
        {
            return Math.Sqrt(Dot(a, a));
        }

        private static Vector Cross(Vector a, Vector b)
        {
            return new Vector(a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x);
        }


        public static LineCrashAttribute CalculateLineDistance(Vector a1, Vector b1, Vector a2, Vector b2)
        {
            if (a1 != null && b1 != null && a2 != null && b2 != null)
            {

                Vector v1 = (a1 - b1);
                Vector v2 = (a2 - b2);

                //向量外积 = 0，则平行
                double value = GetLength(Cross(v1, v2));
                if (value == 0)
                {
                    Vector Q = (a1 - a2);
                    Vector CrossValue = Cross(v1, Q);
                    double Nor = GetLength(CrossValue);
                    double NorOne = GetLength(v1);
                    double dis = Nor / NorOne;

                    return new LineCrashAttribute(null, null, dis, LineState.InterfaceParallel);
                }
                else
                {
                    Vector N = Cross(v1, v2);
                    Vector ab = (a1 - a2);
                    double ans = Dot(N, ab) / GetLength(N);
                    Vector p1 = a1;
                    Vector p2 = a2;
                    Vector d1 = b1 - a1, d2 = b2 - a2;
                    Vector ans1, ans2;
                    double t1, t2;
                    t1 = Dot((Cross(p2 - p1, d2)), Cross(d1, d2));
                    t2 = Dot((Cross(p2 - p1, d1)), Cross(d1, d2));
                    double dd = GetLength((Cross(d1, d2)));
                    t1 /= dd * dd;
                    t2 /= dd * dd;
                    ans1 = (a1 + (b1 - a1) * t1);
                    ans2 = (a2 + (b2 - a2) * t2);
                    ans = Math.Abs(ans);

                    if (ans != 0)
                    {
                        return new LineCrashAttribute(new Point3d(ans1.x, ans1.y, ans1.z), new Point3d(ans2.x, ans2.y, ans2.z), ans, LineState.DifferentSurface);
                    }
                    else
                    {
                        double L1 = Math.Sqrt(Math.Pow(a1.x - a2.x, 2) + Math.Pow(a1.y - a2.y, 2) + Math.Pow(a1.z - a2.z, 2));
                        double L2 = Math.Sqrt(Math.Pow(a1.x - b2.x, 2) + Math.Pow(a1.y - b2.y, 2) + Math.Pow(a1.z - b2.z, 2));
                        double L3 = Math.Sqrt(Math.Pow(b1.x - b2.x, 2) + Math.Pow(b1.y - b2.y, 2) + Math.Pow(b1.z - b2.z, 2));
                        double L4 = Math.Sqrt(Math.Pow(b1.x - a2.x, 2) + Math.Pow(b1.y - a2.y, 2) + Math.Pow(b1.z - a2.z, 2));

                        if (L1 < L2 && L1 < L3 && L1 < L4)
                        {
                            return new LineCrashAttribute(new Point3d(a1.x, a1.y, a1.z), new Point3d(a2.x, a2.y, a2.z), L1, LineState.InterfaceNoParallel);
                        }

                        if (L2 < L1 && L2 < L3 && L2 < L4)
                        {
                            return new LineCrashAttribute(new Point3d(a1.x, a1.y, a1.z), new Point3d(b2.x, b2.y, b2.z), L2, LineState.InterfaceNoParallel);
                        }

                        if (L3 < L1 && L3 < L2 && L3 < L4)
                        {
                            return new LineCrashAttribute(new Point3d(b1.x, b1.y, b1.z), new Point3d(b2.x, b2.y, b2.z), L3, LineState.InterfaceNoParallel);
                        }

                        if (L4 < L1 && L4 < L3 && L4 < L2)
                        {
                            return new LineCrashAttribute(new Point3d(b1.x, b1.y, b1.z), new Point3d(a2.x, a2.y, a2.z), L4, LineState.InterfaceNoParallel);
                        }

                        return null;
                    }
                }
            }
            else
            {
                return null;
            }


        }

        public static Vector operator +(Vector a, Vector b)
        {
            return new Vector(a.x + b.x, a.y + b.y, a.z + b.z);
        }


        public static Vector operator -(Vector a, Vector b)
        {
            return new Vector(a.x - b.x, a.y - b.y, a.z - b.z);
        }
        public static Vector operator *(Vector a, double p)
        {
            return new Vector(a.x * p, a.y * p, a.z * p);
        }

        public static Vector operator /(Vector a, double p)
        {
            return new Vector(a.x / p, a.y / p, a.z / p);
        }

    }

    public class Point3d
    {
        public double x { get; set; }
        public double y { get; set; }
        public double z { get; set; }

        public Point3d(double xx, double yy, double zz)
        {
            x = xx;
            y = yy;
            z = zz;
        }

    }

    public class LineCrashAttribute
    {
        public LineState LineState { get; set; }
        public Point3d ShortestP1 { get; set; }
        public Point3d ShortestP2 { get; set; }
        public double Length { get; set; }

        public LineCrashAttribute(Point3d p1, Point3d p2, double L, LineState S)
        {
            ShortestP1 = p1;
            ShortestP2 = p2;
            Length = L;
            LineState = S;
        }
    }

    public enum LineState
    {
        /// <summary>
        /// 共面平行
        /// </summary>
        InterfaceParallel = 1,
        /// <summary>
        /// 共面不平行
        /// </summary>
        InterfaceNoParallel = 2,
        /// <summary>
        /// 异面直线
        /// </summary>
        DifferentSurface = 3,

    }