《python算法教程》Day11 - 分治法求解平面凸包问题

这是《python算法教程》的第11篇读书笔记,笔记主要内容是使用分治法求解凸包。

平面凸包问题简介

在一个平面点集中,寻找点集最外层的点,由这些点所构成的凸多边形能将点集中的所有点包围起来。
如下图所示,红色的点能将点集中所有的点包围起来。


convexHull.png

分治法求解思路

按照暴力法的思路(求出所有由点集任意两点的直线,再获取使得点集剩余的点在该直线的一侧的直线)去求解凸包问题,显然算法复杂度达到了n^3,这并不是在时间复杂度上可以接受的算法。
因此,可考虑使用分治法去求解凸包。大体思路如下:
1.找出由横坐标最大、最小的两个点p1p2所组成的直线。用该直线将点集分成上下两set1,set2部分。
2.分别从set1、set2找出与线段p1p2构成的面积最大的三角形的点p3,p4。
3.从set1找出在直线p1p3左侧的点集leftset1、在直线p3p2右侧的点集[图片上传中...(行列式.JPG-bc60bb-1525191974104-0)]
rightset1。
4将leftset1,leftset2重复2、3步骤,直至找不到在直线更外侧的点。
5.从set2找出在直线p1p4左侧的点集leftset2、在直线p3p4右侧的点集rightset2。
6.将leftset1,leftset2重复2、3步骤,直至找不到在直线更外侧的点。

点与直线的位置判断

可通过以下行列式的正负值判断直线与点之间的位置关系,同时数值为点与线段所围成的三角形的面积:


image.png

下图表明了若点在直线外围(图中用线段表示直线),上述行列式的值的正负性。
有一点需要注意,下图成立的前提条件是组成直线的两个点(x1,y1)和(x2,y2)必须满足x1<x2,点(x3,y3)必须是判断与直线关系的点。


position.jpg

代码示例

下面的代码示例中加入了绘制散点图的代码,便于观察每一步的情况以及查看最终结果。

#递归法求解凸包
import random
import matplotlib.pyplot as plt


#通过计算三角形p1p2p3的面积(点在直线左边结果为正,直线右边结果为负)来判断 p3相对于直线p1p2的位置
def calTri(p1,p2,p3):
    size=p1[0]*p2[1]+p2[0]*p3[1]+p3[0]*p1[1]-p3[0]*p2[1]-p2[0]*p1[1]-p1[0]*p3[1]
    return size


#找出据直线最远的点(该点与直线围成的三角形的面积为正且最大)
def maxSize(seq,dot1,dot2,dotSet):
    maxSize=float('-inf')
    maxDot=()
    online=[]
    maxSet=[]
    for u in seq:
        size=calTri(dot1,dot2,u)
        #判断点u是否能是三角形u dot1 dot2 的面积为正
        if size<0:
            continue
        elif size==0:
            online.append(u)
        #若面积为正,则判断是否是距离直线最远的点
        if size>maxSize:
            if len(maxDot)>0:
                maxSet.append(maxDot)
            maxSize=size
            maxDot=u
        else:
            maxSet.append(u)
    #结果判断
    #maxSet为空
    if not maxSet:
        #没找到分割点,同时可能有点落在直线dot1 dot2上
        if not maxDot:
            dotSet.extend(online)
            return [],()
        #有分割点
        else:
            dotSet.append(maxDot)
            return [],maxDot
    #maxSet不为空
    else:
        dotSet.append(maxDot)
        return maxSet,maxDot


#找出据直线最远的点(该点与直线围成的三角形的面积为负数且最大)
def minSize(seq,dot1,dot2,dotSet):
    minSize=float('inf')
    minDot=()
    online=[]
    minSet=[]
    for u in seq:
        size=calTri(dot1,dot2,u)
        #判断点u是否能是三角形u dot1 dot2 的面积为负
        if size>0:
            continue
        elif size==0:
            online.append(u)
         #若面积为负,则判断是否是距离直线最远的点
        if size<minSize:
            if len(minDot)>0:
                minSet.append(minDot)
            minDot=u
            minSize=size
        else:
            minSet.append(u)
    #结果判断
    #maxSet为空
    if not minSet:
        #没找到分割点,同时可能有点落在直线dot1 dot2上
        if not minDot:
            dotSet.extend(online)
            return [],()
        #有分割点
        else:
            dotSet.append(minDot)
            return [],minSet
    #maxSet不为空
    else:
        dotSet.append(minDot)
        return minSet,minDot

    
    
#上包的递归划分
def divideUp(seq,dot1,dot2,dot3,dot4,dotSet=None):
    print(dot1,dot2,dot3,dot4)
    #初始化第一次运行时的参数
    if len(seq)==0:
        return dotSet
    if dotSet is None:
        dotSet=[]
    if len(seq)==1:
        dotSet.append(seq[0])
        return dotSet
    leftSet,rightSet=[],[]
    #划分上包左边的点集
    leftSet,maxDot=maxSize(seq,dot1,dot2,dotSet)
    
    #绘图检测---------------------------------------------------------------
    plt.title('up_left')
    #plt.axis([-20,20,-20,20])
    plt.axis([-1100,1100,-1100,1100])
    #plt.scatter([d[0] for d in seq0],[d[1] for d in seq0],color='black')
    plt.scatter([d[0] for d in seq],[d[1] for d in seq],color='blue')
    plt.scatter([dot1[0],dot2[0]],[dot1[1],dot2[1]],color='orange')
    if maxDot:
        plt.scatter(maxDot[0],maxDot[1],color='red')
    plt.show()
    #----------------------------------------------------------------------
    
    #对上包左包的点集进一步划分
    if leftSet:
        divideUp(leftSet,dot1,maxDot,maxDot,dot2,dotSet)
    
    #划分上包右边的点集
    rightSet,maxDot=maxSize(seq,dot3,dot4,dotSet)
        
    #绘图检测---------------------------------------------------------------
    plt.title('up_right')
    #plt.axis([-20,20,-20,20])
    plt.axis([-1100,1100,-1100,1100])
    #plt.scatter([d[0] for d in seq0],[d[1] for d in seq0],color='black')
    plt.scatter([d[0] for d in seq],[d[1] for d in seq],color='blue')
    plt.scatter([dot3[0],dot4[0]],[dot3[1],dot4[1]],color='orange')
    if maxDot:
        plt.scatter(maxDot[0],maxDot[1],color='red')
    plt.show()
    #----------------------------------------------------------------------
    
    #对上包右包的点集进一步划分
    if rightSet:
        divideUp(rightSet,dot3,maxDot,maxDot,dot4,dotSet)
    
    return dotSet


#下包的递归划分
def divideDown(seq,dot1,dot2,dot3,dot4,dotSet=None):
    #初始化第一次运行时的参数
    if len(seq)==0:
        return dotSet
    if dotSet is None:
        dotSet=[]
    if len(seq)==1:
        dotSet.append(seq[0])
        return dotSet
    leftSet,rightSet=[],[]
    #划分下包左边的点集
    leftSet,minDot=minSize(seq,dot1,dot2,dotSet)

    
    #绘图检测---------------------------------------------------------------
    plt.title('down_left')
    #plt.axis([-20,20,-20,20])
    plt.axis([-1100,1100,-1100,1100])
    #plt.scatter([d[0] for d in seq0],[d[1] for d in seq0],color='black')
    plt.scatter([d[0] for d in seq],[d[1] for d in seq],color='blue')
    plt.scatter([dot1[0],dot2[0]],[dot1[1],dot2[1]],color='orange')
    if minDot:
        plt.scatter(minDot[0],minDot[1],color='red')
    plt.show()
    #----------------------------------------------------------------------
    
    #对下包的左包进行进一步划分
    if leftSet:
        divideDown(leftSet,dot1,minDot,minDot,dot2,dotSet)
    
    #划分下包右包的点集
    rightSet,minDot=minSize(seq,dot3,dot4,dotSet)
        
    #绘图检测---------------------------------------------------------------
    plt.title('down_right')
    #plt.axis([-20,20,-20,20])
    plt.axis([-1100,1100,-1100,1100])
    #plt.scatter([d[0] for d in seq0],[d[1] for d in seq0],color='black')
    plt.scatter([d[0] for d in seq],[d[1] for d in seq],color='blue')
    plt.scatter([dot3[0],dot4[0]],[dot3[1],dot4[1]],color='orange')
    if minDot:
        plt.scatter(minDot[0],minDot[1],color='red')
    plt.show()
    #----------------------------------------------------------------------
    
    #对下包的右包进一步划分
    if rightSet:
        divideDown(rightSet,dot3,minDot,minDot,dot4,dotSet)
    
    return dotSet
    
    
#递归主函数
def mainDivide(seq):
    #将序列中的点按横坐标升序排序
    seq.sort()
    res=[]
    #获取横坐标做大、最小的点及横坐标中位数
    dot1=seq[0]
    dot2=seq[-1]
    seq1=[]
    maxSize=float('-inf')
    maxDot=()
    seq2=[]
    minSize=float('inf')
    minDot=()
    #med_x=(seq[len(seq)//2][0]+seq[-len(seq)//2-1][0])/2
    #对序列划分为直线dot1 dot2左右两侧的点集并找出两个点集的距直线最远点
    for u in seq[1:-1]:
        size=calTri(dot1,dot2,u)
        if size>0:
            if size>maxSize:
                if len(maxDot)>0:
                    seq1.append(maxDot)
                maxSize=size
                maxDot=u
                continue
            else:
                seq1.append(u)
        elif size<0:
            if size<minSize:
                if len(minDot)>0:
                    seq2.append(minDot)
                minSize=size
                minDot=u
                continue
            else:
                seq2.append(u)
    print('seq1',seq1,maxDot)
    print('seq2',seq2,minDot)
    #调用内建递归函数
    res1=divideUp(seq1,dot1,maxDot,maxDot,dot2)
    res2=divideDown(seq2,dot1,minDot,minDot,dot2)
    if res1 is not None:
        res.extend(res1)
    if res2 is not None:
        res.extend(res2)
    for u in [dot for dot in [dot1,dot2,maxDot,minDot] if len(dot)>0]:
        res.append(u)
    return res
    

seq0=[(random.randint(-1000,1000),random.randint(-1000,1000)) for x in range(20)]
seq0=list(set(seq0))
res=mainDivide(seq0)
print('res',sorted(res))
plt.axis([-1100,1100,-1100,1100])
plt.title("overview")
plt.scatter([dot[0] for dot in seq0],[dot[1] for dot in seq0],color='black')
plt.scatter([dot[0] for dot in res],[dot[1] for dot in res],color='red')
plt.show()

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 159,015评论 4 362
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 67,262评论 1 292
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 108,727评论 0 243
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 43,986评论 0 205
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 52,363评论 3 287
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 40,610评论 1 219
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 31,871评论 2 312
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 30,582评论 0 198
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 34,297评论 1 242
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 30,551评论 2 246
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 32,053评论 1 260
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 28,385评论 2 253
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 33,035评论 3 236
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 26,079评论 0 8
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 26,841评论 0 195
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 35,648评论 2 274
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 35,550评论 2 270

推荐阅读更多精彩内容

  • 1.概念 凸包(Convex Hull)是一个计算几何(图形学)中的概念。用不严谨的话来讲,给定二维平面上的点集,...
    三三de酒阅读 3,665评论 0 1
  • 归去来兮。 1.1 说明 本篇为《挑战程序设计竞赛(第2版)》[http://www.ituring.com.cn...
    尤汐Yogy阅读 14,114评论 0 160
  • Divide-and-Conquer算法的设计 设计过程分为三个阶段: Divide:整个问题划分为多个子问题 C...
    三三de酒阅读 3,131评论 0 4
  • 判断点是否在线段上、判断两条线段是否相交 这里采用向量的解法。有2个概念:向量的内积和外积。内积又称为点积dot ...
    fruits_阅读 830评论 0 2
  • 非常感谢蒙蒙提出的这个接龙的建议,演讲群也用到了这种,群的参与度一下提升了很多。每个人都留下痕迹,知道自己的位置。...
    24K超超老师阅读 522评论 0 0