转载:算法导论系列--计算几何(computational geometry)

转自:http://blog.csdn.net/daly888/article/details/1486085  ,略作补充

计算几何中的对象可以是点集，线段集，面集等。如果一条线段的端点是有次序之分的，我们把这种线段成为有向线段 (directed segment) 。如果有向线段 p1p2 的起点 p1 在坐标原点，我们可以把它称为向量/矢量 (vector)p2 。

矢量运算：

   在计算几何问题中，经常要用到矢量的运算。设二维矢量 P = ( x1, y1 ) ， Q = ( x2 , y2 ) ，则矢量加法定义为： P + Q = ( x1 + x2 , y1 + y2 ) ，同样的，矢量减法定义为： P - Q = ( x1 - x2 , y1 - y2 ) 。显然有性质 P + Q = Q + P ， P - Q = - ( Q - P ) 。

矢量点乘(dot product) 也叫做向量内积(inner product). 设三维空间中两个向量，P(a1,a2,a3) , Q(b1,b2,b3)， P·Q ＝ a1b1 ＋ a2b2 ＋ a3b3. 几何意义是|P||Q| cos(PQ夹角)。
     两向量点乘为0时，即这两向量互相垂直

矢量叉积(cross product)  设矢量 P = ( x1, y1 ) ， Q = ( x2, y2 ) ，则矢量叉积定义为由 (0,0) 、 p1 、 p2 和 p1+p2 所组成的平行四边形的带符号的面积，即： P × Q = x1*y2 - x2*y1 ，其结果是一个标量(更高维数的叉积参阅数学书)，几何意义为|P||Q| sin(PQ夹角)。显然有性质 P × Q = - ( Q × P ) 和 P × ( - Q ) = - ( P × Q ) 。一般在不加说明的情况下，下述算法中所有的点都看作矢量，两点的加减法就是矢量相加减，而点的乘法则看作矢量叉积。

    设点为 Q ，线段为 AB ，判断点 Q 在直线上的依据是∠QAB = 0或180： ( Q - A ) × ( B - A ) = 0。再判断点是否在线段内，min(xA,xB) <= xQ <= max(xA,xB)且min(yA,yB) <= yQ <= max(yA,yB)。(考虑垂直和水平线段，所以要两个条件同时成立).

线段/直线交点

   设线段两点坐标P1(x1,y1), P2(x2,y2). 此线段构成的直线方程可表示为:

    Ax + By = C

    其中

A = y2-y1
B = x1-x2
C = A*x1+B*y1

因此两条线段可确定两条直线方程

A₁x + B₁y = C₁

A₂x + B₂y = C₂

则解直线方程即可求出直线交点

det = A1*B2 – A2*B1

if (det == 0)

{ 线段平行 }

    else   // 求出交点坐标

      { x = B2*C1 – B1*C2 / det

         y = A1*C2 – A2*C1 / det

      }

    对于线段，只需判断交点坐标是否都在两线段中即可。

点到线段距离

    设点为C，线段AB，C到AB距离为三角形CAB的高，假设C在AB或延长线上的投影为H。我们知道三角形面积有两种计算方法： (1) 1/2 （AB*CH）, 还可以1/2 AB*AC*sin∠CAB = 1/2 |AC×AB| . 因此距离CH = |AC×AB| / |AB| 。

点在线段的垂足

   接上一个话题，寻找H的坐标，线段AB确定直线Px＋Qy = R, 则其垂线方程为-Qx+Py = T. 把C点坐标带入算出T的值得出方程，再求出直线交点即可。若想求C关于AB的轴对称点S，则 S = H – (C – H)

点到线段最近点

    点C到线段AB的最近点，首先应判断C是否在线段AB上, 若是则最近点为C。再判断C在AB的哪一侧。可用点乘判断夹角的方法。

If (AB·BC >= 0)   // ∠ABC 钝角

 { 在B侧，最近点为B }

if (BA·AC >= 0 ) // ∠BAC钝角

 { 在A侧，最近点为A }

else   //H在AB内

 { 最近点为垂足H，按上面的方法求H }

 

线段拐向问题

   设折线段ABC，求拐向。用叉乘可以判断方向，若AC×AB 大于0，则线段右拐，小于0左拐，等于0共线。

多边形面积

    给出一个点集序列，按顺序把各点连起来组合一个多边形, 求多边形面积。一个多边形可以分割成若干个三角形面积的和。上面求点到直线距离时提到三角形面积可以用向量叉乘计算。由于线段拐向改变使叉乘的正负号改变，而闭合图形的拐向改变次数必为偶数，因此该方法对于凹多边形同样适用(正负抵消，读者可自行证明)。

    假设多边形由点集A_i顺序连成，则多边形面积 S = | ∑(A_i+1-Ai)×（A_i+2-A_i+1） | / 2 . 即选取第一个点开始，按顺序把相邻向量叉乘(注意向量方向，后一个点减前一个点)再求和。

   for ( i from 0~N )

   { 

 vecP <-- A[i+1] – A[i]

         vecQ <-- A[i+2] – A[i+1]

       cross <-- vecP × vecQ

       area <-- area + cross

}

判断点在多边形内

    判断点P是否在顺序点集Ai确定的多边形内。过点P作一条足够长(远长于任一条边)的线段模拟过P点的射线，可随机产生一个足够大的数Q。 判断 PQ与多边形边相交的次数，若偶数则P在多边形外，奇数在多边形内(读者自己证明)。还要判断点P是否在边上。

    Q <-- random()*1000+1000

    for( i from 0~N)

    {

       vec = A[i+1] – A[i]

       if ( IsOnSegment(P, vec) == TRUE) //在线段上

         { return 边上 }

       if ( Intersect ( PQ, vec) == TRUE) //相交

         { cnt <-- cnt + 1 }

     } //结束循环

      if ( cnt 是偶数) { return 外面 }

      if (cnt 是奇数) { return 里面 }

凸包(convex hull)

  一个有序点集Ai, 求出最小的凸多边形，把所有点都包含在内，这个多边形叫做凸包(convex hull)。求凸包的方法很多，这里仅介绍一种。

    对于一个有三个或以上点的点集 Q ， Graham 扫描法的过程如下：

　　令 p0 为 Q 中 Y-X 坐标排序下最小的点   
　　设 <p1,p2,...pm> 为对其余点按以 p0 为中心的极角逆时针排序所得的点集（如果有多个点有相同的极角，除了距 p0 最远的点外全部移除 
　　压 p0 进栈 S 
　　压 p1 进栈 S 
　　压 p2 进栈 S 
    for i ← 3 to m 
      do while 由 S 的栈顶元素的下一个元素、 S 的栈顶元素以及 pi 构成的折线段不拐向左侧 
        对 S 弹栈 
       压 pi 进栈 S 
    return S;

　　此过程执行后，栈 S 由底至顶的元素就是 Q 的凸包顶点按逆时针排列的点序列。需要注意的是，我们对点按极角逆时针排序时，并不需要真正求出极角，只需要求出任意两点的次序就可以了。而这个步骤可以用前述的矢量叉积性质实现。

 

 

另外一些论题

   1。正交区域查找(数据库查询)。数据库貌似与几何风马牛不相及。其实我们可以把数据库记录字段对应一个维度，每一条记录看作是K维空间中的一个点。假设数据表有两个字段，则构成2D空间，正交区间查询就是给定一个条件区间，查找符合条件的记录。其实质是按条件生成的正方形，求正方形内的点(如果有3D空间，则是一个立方体). 其方法是构造一个查找树，具体见计算几何书籍。

2。三角剖分。把区域剖分成三角形集，有着广泛应用。例如用于无线传感网络(sensor network)中的路由算法。

补充:

任意多边形面积=每条边向量的叉乘之和的绝对值的1/2:

选取第一个点开始，按顺序把相邻向量叉乘(注意向量方向，后一个点减前一个点)再求和。一般题目给出的是顺时针或逆时针的点坐标,单纯地给出点坐标并不能唯一地确定多边形.

struct Point 
{ 
 double x, y;
}p[10000];
double polygon_area(int n)
{ 
 double area=0; 
 for(int i=0;i<n-1;++i)
 {
 area+=p[i].x*p[i+1].y-p[i].y*p[i+1].x;
 }
 area+=p[n-1].x*p[0].y-p[n-1].y*p[0].x; //首尾相连
 return fabs(area/2.0); //返回绝对值
}