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

 

地理信息系统软件开发中经常需要求取点、线、面之间的交点、交线、封闭区域面积和闭合集等结果,采用以矢量点乘和叉乘为基础的求取算法符合实际工作中已给出点位置和法向量等条件的情况,效率较高。

首先给出基本公式的推导。

矢量的结合率和交换率:

U+(V+W) = (U+V) + W;

U+V=V+U;

而设线段的起点和终点为P,Q;则中间任意一点R可以表示为:

R = P + r(Q - P); (0<r<1)

由P0,P1和P2三点构成的共面参数方程可表示为如下一种形式:

P(s,t) = (1-s-t)P0+sP1+tP2,其中s,t为参数

矢量的长度定义为:

设 V = (a,b); 则 |V| = sqrt(a*a + b*b); (即各轴数据平方和的开方)

而向量的规范化是 U = V/|V|,使得向量的长度为单位1。

矢量的点乘定义:

V–W = v1w1+v2w2,V= (v1,v2) , W = (w1,w2);

此定义可从二维空间类推到三维和其他高次空间。

我们知道,二维平面上面一点P(x,y)在坐标轴转动某角度α后,新坐标为P(xnew,ynew)

设P点距离原点为r,此点与原X轴夹角为δ,则

x = rcosδ, y = rsinδ;

可求得新坐标为:xnew = rcos(δ-α)=rcosδcosα + rsinδsinα = xcosα + ysinα;

ynew =rsin(δ-α)=rsinδcosα –rcosδsinα = ycosα –xsinα;

下面推导广泛运用的一个公式:

V–W = |V||W|cosϕ,ϕ为两向量的交角;

若V = (v1,0),既此点落在X轴上,则V–W = v1w1 + 0w2 = v1w1;

因为V点落在轴上,则 |V| = v1,|W| cosϕ = w1;所以 V–W = |V||W|cosϕ;

若为一般情况,即V= (v1,v2),则我们可旋转原始坐标轴角度δ后使得V点落在新X轴上,

由于: (v1w1 + v2w2)/|V| = w1cosδ+ w2sinδ = W点在新坐标系下的X轴数值 = |W| cosϕ

因此 v1w1 + v2w2 = |V||W|cosϕ;

在平时的运用中,计算cosϕ = V–W / |V||W| ,可对两向量的交角大小进行判断,若V–W为0,则说明两向量互相垂直;若V–W大于0,则|ϕ|<90度,说明交角为锐角。

设向量V(a,b),对其旋转90度后,得到U= (acos90 + bsin90, ycos90-xsin90) = (b,-a);可由此方法得到向量的垂直向量。

由V–W = |V||W|cosϕ公式的推导中可以类推|V||W|sinϕ等于什么呢?推导如下:

|W|sinϕ = W点在新坐标系下的Y轴数值 = w2cosδ-w1sinδ= w2v1/|V|-w1v2/|V|=

(v1w2-v2w1)/|V|,即

v1w2-v2w1 =  |V||W|sinϕ,ϕ为两向量的交角;

       在平时的运用中,可采用此公式计算三角形和四边形的面积。也能判断两向量的位置关系,如果计算v1w2-v2w1为0,则说明两向量交角为0度或180度,即共线。如果计算结果大于0,则说明两向量交角0度到180度之间。

求取向量之间位置关系和三角形面积的参考代码:

//    Input: three points P0, P1, and P2

//    Return: >0 for P2 left of the line through P0 and P1

//            =0 for P2 on the line

//            <0 for P2 right of the line

double CDEMAlgorithm::isLeft( Point P0, Point P1, Point P2 )

{

    return ( (P1.x - P0.x) * (P2.y - P0.y)

        - (P2.x - P0.x) * (P1.y - P0.y) );       //叉乘

}

double CDEMAlgorithm::orientation2D_Triangle( Point V0, Point V1, Point V2 )

{

    return isLeft(V0, V1, V2);                        //判断位置关系

}

double CDEMAlgorithm::area2D_Triangle( Point V0, Point V1, Point V2 )

{

    return isLeft(V0, V1, V2) / 2.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);        // always >= 0
 double    b = Dot(u,v);
 double    c = Dot(v,v);        // always >= 0
 double    d = Dot(u,w);
 double    e = Dot(v,w);
 double    D = a*c - b*b;       // always >= 0
 double    sc, tc;

 // compute the line parameters of the two closest points
 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);  // = L1(sc) - L2(tc)

 return sqrt(Dot(dP,dP));   // return the closest distance

}
附图:

继续上一节的内容,本节主要讲解三维空间中射线、线段与平面及三维物体的交点及距离的计算,它们在碰撞检测和可见性剔除等应用中是必不可少的。首先给出3D空间下点乘和叉乘的定义与定理的推导,再谈如何应用到程序编码的工作中。

设三维空间中任意两矢量V(v1,v2,v3)和W(w1,w2,w3)

                                                                                                                                    
               由三角形的余弦公式:

                                      cosα = (a2 + b2 – c2)/2ab

               则两矢量夹角α的余弦 : cosα = {(v1*v1 + v2*v2 + v3*v3) + (w1*w1+w2*w2+w3*w3)- [ (w1 - v1)2+( w2 - v2)2+( w3 - v3)2] } / 2 |V| |W| = (v1w1+ v2w2 + v3w3)/ |V| |W|

               由于将矢量的点乘定义为 : V·W = v1w1+ v2w2 + v3w3 ;

               所以 |V| |W| cosα = V·W;

               由上面的结论引申,可知道如果两矢量的点乘为0,则说明夹角α的余弦为0,即矢量的夹角是90度,互相垂直。

               再来看两矢量的叉乘,与点乘的计算结果定义为一数值不同,叉乘的结果定义为另外一个矢量 U (u1,u2,u3) ,其中:

                                             u1 = v2w3 - v3w2;

                                             u2 = v3w1 - v1w3;

                                             u3 = v1w2 - v2w1;

               由这个定义我们来推导:

               V·U = v1(v2w3 - v3w2) + v2(v3w1 - v1w3) + v3(v1w2 - v2w1) = v1v2w3 - v1 v3w2 + v2v3w1-v2 v1w3 + v3v1w2-v3 v2w1 = 0;两矢量的点乘为0,说明它们互相垂直。

               W·U = w1(v2w3 - v3w2) + w2(v3w1 - v1w3) + w3(v1w2 - v2w1) = 0,说明它们也是互相垂直的。

               在三维空间中,U与V、W都垂直,说明U是垂直与V与W构成的平面,也就是这个平面的法向量。

               U虽然是矢量,但其模|U|依旧是一个数值,表明其长度。

               |U|2 = u1*u1 + u2*u2 + u3*u3 = (v2w3 - v3w2)2+( v3w1 - v1w3)2+( v1w2 - v2w1)2;

               同时考虑,在α是V、W两矢量的夹角条件下,(|V| |W|sinα)2= |V|2|W|2(1-cos2α)= |V|2|W|2 - |V|2|W|2cos2α = |V|2|W|2 - (V·W)2= (v1*v1 + v2*v2 + v3*v3)( w1*w1+w2*w2+w3*w3) - (v1w1+ v2w2 + v3w3)2

               通过化解,可以得到一个数学等式,即:

               (bz-cy)2+(cx-az)2+(ay-bx)2=(a2+b2+c2)(x2+y2+z2)-(ax+by+cz)2

               所以 (v2w3 - v3w2)2+( v3w1 - v1w3)2+( v1w2 - v2w1)2  = (v1*v1 + v2*v2 + v3*v3)( w1*w1+w2*w2+w3*w3) - (v1w1+ v2w2 + v3w3)2;

               即                          |U|2 = (|V| |W|sinα)2

                                    |U| = |V| |W|sinα

            在完成三维空间的点乘和叉乘定义与公式推导后,可以直接运用结论到程序的编写中,这些结论包括利用点乘来求得夹角的余弦,利用点乘为0来求得垂直于某向量的另一向量,利用两向量的叉乘来求得平面的法向量,利用叉乘的模来求得三角形的面积等。

            三维空间中点到直线、射线和线段的距离。即从点P作垂直于线段(P0,P1)的垂线,求得垂心T后,计算两点之间的距离。

                                                   

          设置T点为参数表达方式: T = P0 + t(P1 –P0);将T点看作是从P0开始往P1移动的动态点,则矢量PT为PP0与P0T两个矢量之和,当PT垂直于P0P1时,利用两个矢量的点乘为0列等式,求得参数t的数值。

                       (PP0 + P0T) ·P1P0 = 0;

                              即 (P-P0)·(P1 –P0) + t(P1 –P0) ·(P1 –P0) = 0;

                       t = - (P-P0)·(P1 –P0) / (P1 –P0) ·(P1 –P0)

        可以看出分子是需要计算两个矢量的点乘,而分母是矢量P0P1点乘自身,即为模的平方。

double CDEMAlgorithm::dist_Point_to_Line( Point P, Line L)

{

    Vector v = L.P1 - L.P0;                                     //顶点相减得到矢量

    Vector w = P - L.P0;                                    //顶点相减得到矢量

    double c1 = Dot(w,v);                                   //两矢量的点乘

    double c2 = Dot(v,v);                                   //两矢量的点乘

    double b = c1 / c2;                                     //求得参数

    Point Pb = L.P0 + b * v;                                //解得参数后,得到垂心位置

    return d(P, Pb);                                        //返回两点之间的距离

}

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

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

       
       两条直线分别设置为P0P1,P2P3;

       找到直线上的两个点Pt,Ps;Pt= P0 + t(P1 – P0) ;

                                           Ps = P2 + s(P3 – P2) ;

两条直线之间的距离定义为它们之间的最短距离,可得知PtPs垂直于P0P1和P2P3。由上两节讲解得到的结论:两矢量互相垂直,则点乘为0,可得:

                                                   (Ps –Pt) •(P1 – P0) = 0;

                                                   (Ps –Pt) •(P3 – P2) = 0;

令U = P1 – P0;V = P3 – P2;则    (P2 + sU– P0 – tV) •U = 0;

                                                   (P2 + sU– P0 – tV) •V = 0;

再另W = P2 – P0;               则 W•U + sU•U – tV•U = 0;

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

      

                                    考虑到:U•U = |U|2

                                                           V•V = | V |2

                                                            V•U = U•V

解方程组,得系数:

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

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

分母|U|2| V |2 – (U•V)2 若等于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));   // 返回两顶点直接的距离

}

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


    由上两节讲解得到的结论:P0x P1y – P0y P1x = | P0O | | P1O |sina,将此数值除以2,即得到三角形P0O P1的面积,在顶点依照顺或逆时针排好序后,依次计算每两个顶点与原点构成的三角形面积,由于面积带正负号(由方向决定的),求总和(将抵消一部分),最后得到多边形的面积。

参考代码:

//    输入:  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;     // 索引

         

  • 0
    点赞
  • 16
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值