计算点、线、面等元素之间的交点、交线、封闭区域面积和闭合集(续2)

 

在前面两节探讨的基础上,利用已知空间顶点坐标情况下，以矢量点乘和叉乘为基础，继续讲解三维空间中两条直线之间的距离,二维平面上任意凸多边形的面积, 三维空间中凸多边形的面积等方面的应用和程序编写工作。

三维空间中两条直线之间的距离

       
       两条直线分别设置为P₀P₁，P₂P₃；

       找到直线上的两个点P_t，P_s；P_t= P₀ + t(P₁ – P₀) ；

                                           P_s = P₂ + s(P₃ – P₂) ；

两条直线之间的距离定义为它们之间的最短距离，可得知P_tP_s垂直于P₀P₁和P₂P₃。由上两节讲解得到的结论：两矢量互相垂直，则点乘为0，可得：

                                                   (P_s –P_t) •(P₁ – P₀) = 0；

                                                   (P_s –P_t) •(P₃ – P₂) = 0；

令U = P₁ – P₀；V = P₃ – P₂；则    (P₂ + sU– P₀ – tV) •U = 0；

                                                   (P₂ + sU– P₀ – tV) •V = 0；

再另W = P₂ – P₀；               则 W•U + sU•U – tV•U = 0；

W•V + sU•V – tV•V = 0；

      

                                    考虑到：U•U = |U|²

                                                           V•V = | V |²

                                                            V•U = U•V

解方程组，得系数：

                                            s = ((W•V)( U•V) – (W•U) | V |²) / ( |U|²| V |² – (U•V)²)

                                            t = ((W•V) |U|² – (W•U)( U•V) )/ (|U|²| V |² – (U•V)²)

分母|U|²| V |² – (U•V)² 若等于0，则说明两直线的夹角a的余弦为1，即两直线的夹角为0度，为互相平行。此时，不采用此公式计算，只需要在一条直线上任意取得一点，计算此点到另一直线的距离，即为两直线之间的距离。

参考代码：

double CDEMAlgorithm::dist3D_Line_to_Line( Line L1, Line L2)

{

    Vector   u = L1.P1 - L1.P0;

    Vector   v = L2.P1 - L2.P0;

    Vector   w = L1.P0 - L2.P0;

    double    a = Dot(u,u);        // 点乘的结果

    double    b = Dot(u,v);

    double    c = Dot(v,v);       

    double    d = Dot(u,w);

    double    e = Dot(v,w);

    double    D = a*c - b*b;      

    double    sc, tc;

    // 在我们已经推导后，计算参数

    if (D < EPS) {         // 两条线平行

        sc = 0.0;

        tc = (b>c ? d/b : e/c); 

    }

    else {

        sc = (b*e - c*d) / D;

        tc = (a*e - b*d) / D;

    }

    //sc,tc分别指示了在两条直线上的比例参数

    Vector   dP = w + (sc * u) - (tc * v); //得到两顶点构成的矢量

    return sqrt(Dot(dP,dP));   // 返回两顶点直接的距离

}

二维平面上任意凸多边形的面积

    由上两节讲解得到的结论：P₀x P₁y – P₀y P₁x = | P₀O | | P₁O |sina，将此数值除以2，即得到三角形P₀O P₁的面积，在顶点依照顺或逆时针排好序后，依次计算每两个顶点与原点构成的三角形面积，由于面积带正负号（由方向决定的），求总和（将抵消一部分），最后得到多边形的面积。

参考代码：

//    输入:  int n = 多边形的顶点个数

//            Point* V = 记录 n+2 个顶点的数组

//            with V[n]=V[0] and V[n+1]=V[1] //请注意在计算面积时，设置的顶点最后两个与最前两个重复

//    返回: 多边形的面积数值

double CDEMAlgorithm::area2D_Polygon( int n, Point* V )

{

    double area = 0;

    int   i, j, k;     // 索引

                     //(x₀y₁ - x₁y₀) + (x₁y₂ - x₂y₁) + (x₂y₃ - x₃y₂) ....

    for (i=1, j=2, k=0; i<=n; i++, j++, k++) {    //通过合并中间的项，可推出下面使用的公式计算面积（带符号）

        area += V[i].x * (V[j].y - V[k].y);

    }

    return area / 2.0;

}

三维空间中凸多边形的面积

    
     将多边形投影到某坐标面上，转换成2D平面的多边形面积问题后，将求得的投影面积再除以cosa，a为原始多边形与投影多边形的夹角。

   参考代码：

// 参数：顶点N指示平面的法向量

double CDEMAlgorithm::area3D_Polygon( int n, Point* V, Point N )

{

    double area = 0;

    double an, ax, ay, az; 

    int   coord;        

    int   i, j, k;      

    // 选择法向量中的最大数值

    ax = (N.x>0 ? N.x : -N.x);     //绝对值

    ay = (N.y>0 ? N.y : -N.y);     //绝对值

    az = (N.z>0 ? N.z : -N.z);     //绝对值

    coord = 3;                  

    if (ax > ay) {

        if (ax > az) coord = 1;    

    }

    else if (ay > az) coord = 2;   //平面的法向量的哪个轴数值最大，则准备在计算面积时候，投影到另外两根轴构成的平面上。

    for (i=1, j=2, k=0; i<=n; i++, j++, k++)

        switch (coord) {             //首先计算二维投影平面上的面积

        case 1:

            area += (V[i].y * (V[j].z - V[k].z));

            continue;

        case 2:

            area += (V[i].x * (V[j].z - V[k].z));

            continue;

        case 3:

            area += (V[i].x * (V[j].y - V[k].y));

            continue;

    }

    // 还原到原始面积

    an = sqrt( ax*ax + ay*ay + az*az);

    switch (coord) {

    case 1:

        area *= (an / ax);

        break;

    case 2:

        area *= (an / ay);

        break;

    case 3:

        area *= (an / az);  

    }

    return area/2;

}