【转帖】算法导论系列--计算几何(computational geometry)

huangyhg

Computational Geometry是个十分有趣的课题。很遗憾笔者不是数学/计算机专业，对这个领域认识很少。第一次接触这个词是在topcoder的算法教程里面，被其算法的巧妙所吸引，可惜未有机会进一步学习。认识所限，本文介绍仅能最简单的2D几何算法跟大家分享。对于从事CAD软件开发，图形开发，统计，游戏开发的人，计算几何都是必修课。这个领域有一本很著名的入门书：Mark de berg(荷)写<Computational Geometry: Algorithms and application 2nd edition>(国内有出版，计算几何—算法与应用第二版，邓俊辉译) 。
计算几何中的对象可以是点集，线段集，面集等。如果一条线段的端点是有次序之分的，我们把这种线段成为有向线段 (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 ) 。一般在不加说明的情况下，下述算法中所有的点都看作矢量，两点的加减法就是矢量相加减，而点的乘法则看作矢量叉积。

　　二维叉积的一个非常重要性质是可以通过它的符号判断两矢量相互之间的顺逆时针关系：若P×Q > 0 , 则P在Q的顺时针方向。若小于0， 为逆时针。若等于0， 可能同向/反向。
（这段话还真是误导人，三维情况下，叉积的结果是个矢量。二维叉积的结果是个标量是个特例）

判断点在线段上

设点为 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
因此两条线段可确定两条直线方程
A1x + B1y = C1
A2x + B2y = C2
则解直线方程即可求出直线交点
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 AC*sin∠CAB = 1/2 |AC×AB| . 因此距离CH = |AC×AB| / |AB| 。
[url="http://www.dimcax.com/resource/math/image/o_linedist.gif"][url="http://www.dimcax.com/resource/math/image/o_linedist.gif"]http://www.dimcax.com/resource/math/image/o_linedist.gif[/url]







点在线段的垂足

接上一个话题，寻找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共线。

多边形面积

给出一个点集序列，按顺序把各点连起来组合一个多边形, 求多边形面积。一个多边形可以分割成若干个三角形面积的和。上面求点到直线距离时提到三角形面积可以用向量叉乘计算。由于线段拐向改变使叉乘的正负号改变，而闭合图形的拐向改变次数必为偶数，因此该方法对于凹多边形同样适用(正负抵消，读者可自行证明)。
假设多边形由点集Ai顺序连成，则多边形面积 S = | ∑(Ai+1-Ai)×（Ai+2-Ai+1） | / 2 . 即选取第一个点开始，按顺序把相邻向量叉乘(注意向量方向，后一个点减前一个点)再求和。
for ( i from 0~N )
{
vecP <-- A[i+1] – A
vecQ <-- A[i+2] – A[i+1]
cross <-- vecP × vecQ
area <-- area + cross
}
http://www.dimcax.com/resource/math/image/o_polygon.gif








[/url]



[url="http://p.blog.csdn.net/images/p_blog_csdn_net/daly888/272736/o_polygon.gif"][/url]
判断点在多边形内

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

Q <-- random()*1000+1000
for( i from 0~N)
{
vec = A[i+1] – A
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 的凸包顶点按逆时针排列的点序列。需要注意的是，我们对点按极角逆时针排序时，并不需要真正求出极角，只需要求出任意两点的次序就可以了。而这个步骤可以用前述的矢量叉积性质实现。

http://www.dimcax.com/resource/math/image/o_convex2.gif









另外一些论题

1。正交区域查找(数据库查询)。数据库貌似与几何风马牛不相及。其实我们可以把数据库记录字段对应一个维度，每一条记录看作是K维空间中的一个点。假设数据表有两个字段，则构成2D空间，正交区间查询就是给定一个条件区间，查找符合条件的记录。其实质是按条件生成的正方形，求正方形内的点(如果有3D空间，则是一个立方体). 其方法是构造一个查找树，具体见计算几何书籍。
2。三角剖分。把区域剖分成三角形集，有着广泛应用。例如用于无线传感网络(sensor network)中的路由算法。
http://blog.csdn.net/daly888/archive/2007/01/17/1486085.aspx