唐磊的个人博客

点集的凸包问题

凸的意思就是,点集当中任意两个点的连线仍然在这个集合当中。凸包就是说用一个足够大的橡皮筋把这些点包围住,这个最小的橡皮筋就构成了这些点集的凸包。

求出凸包是计算几何当中一个经典的问题,也是贯穿这整个计算几何的学科。如下图所示,最外层的那个圈就是这些点集的凸包。凸包上的那些点称为极点,如果找出了所有的极点,那么,凸包也就求出来了。极点的可以用这种方式来进行判断,点集P中任意一个点s,如果存在另外的三个不同的点p,q,r构成的三角形将s包围住,那么这个点s就不是极点,反之就是极点。于是也就有了第一个算法:

凸包

Mark All points in Set P as extreme;
For each triangle (p,q,r):
  For each k in Set P and k ! = (p,q,r):
       If k in triangle(p,q,r)
Mark k as NON-extreme.

判断一个点k是否在点p,q,r构成的三角形中可以用toleft测试。如上图中的例子,如果k在有向线段pq的左边 and k 在 qr的左边and k在rp的左边,那么k就在pqr构成的的三角形里面了。核心代码如下:

//O(n^4)
vector<CP_Vector2D> cal_extreme_points()
{
vector<CP_Vector2D> extreme_points;
vector<MyPoint> my_points;
unsigned int size = points.size();
//inital all points to extreme
for (unsigned int i = 0; i < size; i++)
my_points.push_back(MyPoint(points[i], true));
for (unsigned int p = 0; p < size; p++)
{
for (unsigned int q = p+1 ; q < size; q++)
{
for (unsigned int r = q+1; r < size; r++)
{
for (unsigned int k = 0; k < size; k++)
{
if (k == p || k == q || k == r || !my_points[k].isExtreme)
continue;
if (in_triangle(my_points[k].point, my_points[p].point, my_points[q].point, my_points[r].point))
my_points[k].isExtreme = false;
}
}
}
}
for (unsigned int i = 0; i < size; i++)
{
if(my_points[i].isExtreme)
extreme_points.push_back(my_points[i].point);
}
//sort the extrme_points
sort(extreme_points.begin(), extreme_points.end(), compare);
leftest_and_lowest = extreme_points[0];
//sort the extrme_points by left_and_lowest polar
sort(extreme_points.begin()+1, extreme_points.end(), compareByAngle);
return extreme_points;
}

这是一个O(n4)的算法,找出所有极点之后,还得用O(nlgn)排一下序将所有的极点连接起来(可以按照极点与最左下角的点和沿水平X轴的夹角排序)。

稍微有些改进的方法是找出所有的极边(构成凸包的多边形的每一条边都是极边),可以这么来计算极边,除了极边两个点外,其他的点都在这个边的同一侧。如下算法:

let EE = null
for each directed segment pq
if points in P\{p,q} lie to the SAME side of pq then let EE = EE+{pq}

找出所有的极边需要O(n3),然后再连接起来,整个算法复杂度为O(n3).

//O(n^3)
vector<CP_Vector2D> cal_extreme_edges()
{
vector<ExtremeEdge> edgs;
unsigned int size = points.size();
for (unsigned int p = 0; p < size; p++)
{
for (unsigned int q = p+1 ; q < size; q++)
{
bool all_left = true, all_right = true;
for (unsigned int k = 0; k < size; k++)
{
if (k == p || k == q || p == q)
continue;
if(to_left(points[k], points[p], points[q]))
all_left = false;
else
all_right = false;
}
//let the direction be all the same
if (all_left)
edgs.push_back(ExtremeEdge(points[p], points[q]));
if (all_right)
edgs.push_back(ExtremeEdge(points[q], points[p]));
}
}
vector<CP_Vector2D> convexPoints;
////O(n^2), should better sort the edges by angle.
//convexPoints.push_back(edgs[0].end);
//while(convexPoints.size() != edgs.size())
//{
// for (unsigned int i = 0; i < edgs.size(); i++)
// {
// if (convexPoints[convexPoints.size() - 1] == edgs[i].start)
// convexPoints.push_back(edgs[i].end);
// }
//}
//O(nlogn)
sort(edgs.begin(), edgs.end(), compareEdgeByAngle);
for (unsigned int i=0; i < edgs.size(); i++)
{
convexPoints.push_back(edgs[i].start);
}
return convexPoints;
}

再说一下Gift Wrapping算法,算法思路如下:因为所有的极边都通过端点连接而成围成一个圈,一旦找到第一条极边,可以通过利用O(n)搜索剩下的点构成接下来的极边,圈围成后,凸包即计算出来。整个算法复杂度O(n2)。其中第一个点可以通过O(n)的时间找到最左下角的点。就像其名字一样,整个算法就像在包装gift。

find an arbitrary EP, denoted as p0
let i = 0
repeat
find the EE pipi+1
i++
until pi = p0

分析下GiftWrapping算法,每条极边的确定需要O(N)的时间,算法所耗费的时间依赖与极边的数量,因此GiftWrapping算法在最坏的情况下需要O(n2),最好的情况下仅仅需要O(N),更确切地讲,需要O(n*h),h是凸包包含的点数量,这个算法是一个output sensitive的算法。

//O(n^2),in fact O(n*h),h = convex hull.size, output sensitive
vector<CP_Vector2D> giftwraping()
{
int ltl = 0;//find the lowest-and-leftmost point
for (unsigned int i = 1; i < points.size(); i++)
{
if(points[i].m_y < points[ltl].m_y ||
(points[i].m_y == points[ltl].m_y && points[i].m_x < points[ltl].m_x))
ltl = i;
}
vector<CP_Vector2D> convexPoints;
convexPoints.push_back(points[ltl]);
for(int p = ltl; ;)
{
int q = -1;
for(unsigned int k = 0; k < points.size(); k++)
{
if( (k != p) && (q < 0 || !to_left(points[k], points[p], points[q])))
q = k;//update q if k lies in right of pq
}
if(ltl == q)
break;// has find out the convex hull circle
//find q
convexPoints.push_back(points[q]);
p = q;
}
return convexPoints;
}

是否能够继续优化下去,将算法复杂度再次提高?或者Convex Hull的解法是否有一个下界?下面的算法就是一个总体复杂度O(nlgn)的算法,已经有人证明,这个是凸包问题的下界。不能再继续优化了。下面就看看著名的O(nlgn)的Graham Scan算法。


算法思路如下:首先仍然是找到最左下角的点,然后按照极角排序,如上图所示,编上号1——n,然后从第三个点开始,依次看该点是否在前两个点构成的直线的左边(3是否在有向线段12的左边,4是否在有向线段23的左边,依次类推),上例中,到第5个点时,发现5不在34的左边,则回退一个,连接35,到第6个点时,发现6也不在前两个点35(此时4已经被排除)构成的有向线段的左边,再回退,看6仍然不在23的左边,继续回退,直到发现6在12的左边,此时连接26……就按照这样走下去,完成了整个凸包的构成。从具体实现上可以用两个栈来实现,算法如下所示:

find the LTL point p1
sort all points w.r.t. p1 by polar angle: {p2, p3, ..., pn}
create stack S; S.push(p1) and S.push(p2)
create stack T; for i = n down to 3 do T.push(pi)
while (!T.empty()) do
while (!toLeft(S[1], S[0], T[0])) do
S.pop()
S.push(T.pop())

最后S中的点就是凸包上的极点,并且是有序的。注意上述算法中,每一个从S中pop出来的点都不是极点,因为会被其他三角形包围住。而所有的点都会在p1和p2的左边。这个代码用stl的stack有些别扭啊。因为stack不能直接查看栈顶的第二个元素….不要紧,意思意思即可。:)

//O(nlgn)
vector<CP_Vector2D> graham_scan()
{
int ltl = 0;//find the lowest-and-leftmost point
for (unsigned int i = 1; i < points.size(); i++)
{
if(points[i].m_y < points[ltl].m_y ||
(points[i].m_y == points[ltl].m_y && points[i].m_x < points[ltl].m_x))
ltl = i;
}
leftest_and_lowest = points[ltl];
//compareByAngle should use the ltl point
sort(points.begin(), points.end(), compareByAngle);
vector<CP_Vector2D> convexPoints;
stack<CP_Vector2D> S, T;
S.push(points[0]);S.push(points[1]);
for (int i = points.size()-1; i >= 2;i--)
{
T.push(points[i]);
}
CP_Vector2D s0, s1;
while (!T.empty())
{
s0 = S.top();
S.pop();//s0 is gone
if(S.empty())
{
S.push(T.top());
break;
}
s1 = S.top();//get s1
S.push(s0);
while (!to_left(T.top(), s1, s0))
{
S.pop();
s0 = S.top();
S.pop();//s0 is gone
s1 = S.top();//get s1
S.push(s0);
}
S.push(T.top());
T.pop();
}
int size = S.size();
for(int i = 0; i < size; i++)
{
convexPoints.push_back(S.top());
S.pop();//size is change, cannot use i < S.size()
}
return convexPoints;
}

来分析下这个算法的复杂度。刚开始用O(N)的时间找p1,然后排序O(nlgn),最复杂的就是这个Scan,表面上看有两层while循环,仔细分析发现:toleft测试的数量=T.pop+s.push =O(N),而每一次toleft测试,总有一个点会从S或者T中pop出来,每个点从S中push进去一次且最多被pop出来一次。因此,整个while循环那里实际只耗费了O(N)。所以整个算法复杂度就是O(nlgn)。

最后自己实现了下以上几个算法,用MFC+OPENGL如图。注意实现时细节的处理,比如按照极角排序等,若考虑算法鲁棒性,得处理多点共线、重点等情况(本文未考虑)。

参考资料:邓俊辉 《计算几何》

tanglei wechat
欢迎扫描二维码关注我的微信公众号