這是《python算法教程》的第11篇讀書(shū)筆記醉蚁，筆記主要內(nèi)容是使用分治法求解凸包。

平面凸包問(wèn)題簡(jiǎn)介

在一個(gè)平面點(diǎn)集中，尋找點(diǎn)集最外層的點(diǎn)煌珊，由這些點(diǎn)所構(gòu)成的凸多邊形能將點(diǎn)集中的所有點(diǎn)包圍起來(lái)私杜。
如下圖所示蚕键，紅色的點(diǎn)能將點(diǎn)集中所有的點(diǎn)包圍起來(lái)。

convexHull.png

分治法求解思路

按照暴力法的思路（求出所有由點(diǎn)集任意兩點(diǎn)的直線(xiàn)衰粹，再獲取使得點(diǎn)集剩余的點(diǎn)在該直線(xiàn)的一側(cè)的直線(xiàn)）去求解凸包問(wèn)題锣光，顯然算法復(fù)雜度達(dá)到了n^3，這并不是在時(shí)間復(fù)雜度上可以接受的算法铝耻。
因此誊爹，可考慮使用分治法去求解凸包。大體思路如下：
1.找出由橫坐標(biāo)最大瓢捉、最小的兩個(gè)點(diǎn)p1p2所組成的直線(xiàn)频丘。用該直線(xiàn)將點(diǎn)集分成上下兩set1，set2部分泡态。
2.分別從set1搂漠、set2找出與線(xiàn)段p1p2構(gòu)成的面積最大的三角形的點(diǎn)p3,p4。
3.從set1找出在直線(xiàn)p1p3左側(cè)的點(diǎn)集leftset1某弦、在直線(xiàn)p3p2右側(cè)的點(diǎn)集[圖片上傳中...(行列式.JPG-bc60bb-1525191974104-0)]
rightset1状答。
4將leftset1,leftset2重復(fù)2、3步驟刀崖，直至找不到在直線(xiàn)更外側(cè)的點(diǎn)惊科。
5.從set2找出在直線(xiàn)p1p4左側(cè)的點(diǎn)集leftset2、在直線(xiàn)p3p4右側(cè)的點(diǎn)集rightset2亮钦。
6.將leftset1,leftset2重復(fù)2馆截、3步驟，直至找不到在直線(xiàn)更外側(cè)的點(diǎn)蜂莉。

點(diǎn)與直線(xiàn)的位置判斷

可通過(guò)以下行列式的正負(fù)值判斷直線(xiàn)與點(diǎn)之間的位置關(guān)系蜡娶，同時(shí)數(shù)值為點(diǎn)與線(xiàn)段所圍成的三角形的面積：

image.png

下圖表明了若點(diǎn)在直線(xiàn)外圍（圖中用線(xiàn)段表示直線(xiàn)），上述行列式的值的正負(fù)性映穗。
有一點(diǎn)需要注意窖张，下圖成立的前提條件是組成直線(xiàn)的兩個(gè)點(diǎn)（x1，y1）和（x2蚁滋，y2）必須滿(mǎn)足x1<x2宿接，點(diǎn)（x3赘淮，y3）必須是判斷與直線(xiàn)關(guān)系的點(diǎn)。

position.jpg

代碼示例

下面的代碼示例中加入了繪制散點(diǎn)圖的代碼睦霎，便于觀(guān)察每一步的情況以及查看最終結(jié)果梢卸。

#遞歸法求解凸包
import random
import matplotlib.pyplot as plt


#通過(guò)計(jì)算三角形p1p2p3的面積（點(diǎn)在直線(xiàn)左邊結(jié)果為正，直線(xiàn)右邊結(jié)果為負(fù)）來(lái)判斷 p3相對(duì)于直線(xiàn)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


#找出據(jù)直線(xiàn)最遠(yuǎn)的點(diǎn)（該點(diǎn)與直線(xiàn)圍成的三角形的面積為正且最大）
def maxSize(seq,dot1,dot2,dotSet):
    maxSize=float('-inf')
    maxDot=()
    online=[]
    maxSet=[]
    for u in seq:
        size=calTri(dot1,dot2,u)
        #判斷點(diǎn)u是否能是三角形u dot1 dot2 的面積為正
        if size<0:
            continue
        elif size==0:
            online.append(u)
        #若面積為正副女，則判斷是否是距離直線(xiàn)最遠(yuǎn)的點(diǎn)
        if size>maxSize:
            if len(maxDot)>0:
                maxSet.append(maxDot)
            maxSize=size
            maxDot=u
        else:
            maxSet.append(u)
    #結(jié)果判斷
    #maxSet為空
    if not maxSet:
        #沒(méi)找到分割點(diǎn),同時(shí)可能有點(diǎn)落在直線(xiàn)dot1 dot2上
        if not maxDot:
            dotSet.extend(online)
            return [],()
        #有分割點(diǎn)
        else:
            dotSet.append(maxDot)
            return [],maxDot
    #maxSet不為空
    else:
        dotSet.append(maxDot)
        return maxSet,maxDot


#找出據(jù)直線(xiàn)最遠(yuǎn)的點(diǎn)（該點(diǎn)與直線(xiàn)圍成的三角形的面積為負(fù)數(shù)且最大）
def minSize(seq,dot1,dot2,dotSet):
    minSize=float('inf')
    minDot=()
    online=[]
    minSet=[]
    for u in seq:
        size=calTri(dot1,dot2,u)
        #判斷點(diǎn)u是否能是三角形u dot1 dot2 的面積為負(fù)
        if size>0:
            continue
        elif size==0:
            online.append(u)
         #若面積為負(fù)蛤高，則判斷是否是距離直線(xiàn)最遠(yuǎn)的點(diǎn)
        if size<minSize:
            if len(minDot)>0:
                minSet.append(minDot)
            minDot=u
            minSize=size
        else:
            minSet.append(u)
    #結(jié)果判斷
    #maxSet為空
    if not minSet:
        #沒(méi)找到分割點(diǎn),同時(shí)可能有點(diǎn)落在直線(xiàn)dot1 dot2上
        if not minDot:
            dotSet.extend(online)
            return [],()
        #有分割點(diǎn)
        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)
    #初始化第一次運(yùn)行時(shí)的參數(shù)
    if len(seq)==0:
        return dotSet
    if dotSet is None:
        dotSet=[]
    if len(seq)==1:
        dotSet.append(seq[0])
        return dotSet
    leftSet,rightSet=[],[]
    #劃分上包左邊的點(diǎn)集
    leftSet,maxDot=maxSize(seq,dot1,dot2,dotSet)
    
    #繪圖檢測(cè)---------------------------------------------------------------
    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()
    #----------------------------------------------------------------------
    
    #對(duì)上包左包的點(diǎn)集進(jìn)一步劃分
    if leftSet:
        divideUp(leftSet,dot1,maxDot,maxDot,dot2,dotSet)
    
    #劃分上包右邊的點(diǎn)集
    rightSet,maxDot=maxSize(seq,dot3,dot4,dotSet)
        
    #繪圖檢測(cè)---------------------------------------------------------------
    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()
    #----------------------------------------------------------------------
    
    #對(duì)上包右包的點(diǎn)集進(jìn)一步劃分
    if rightSet:
        divideUp(rightSet,dot3,maxDot,maxDot,dot4,dotSet)
    
    return dotSet


#下包的遞歸劃分
def divideDown(seq,dot1,dot2,dot3,dot4,dotSet=None):
    #初始化第一次運(yùn)行時(shí)的參數(shù)
    if len(seq)==0:
        return dotSet
    if dotSet is None:
        dotSet=[]
    if len(seq)==1:
        dotSet.append(seq[0])
        return dotSet
    leftSet,rightSet=[],[]
    #劃分下包左邊的點(diǎn)集
    leftSet,minDot=minSize(seq,dot1,dot2,dotSet)

    
    #繪圖檢測(cè)---------------------------------------------------------------
    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()
    #----------------------------------------------------------------------
    
    #對(duì)下包的左包進(jìn)行進(jìn)一步劃分
    if leftSet:
        divideDown(leftSet,dot1,minDot,minDot,dot2,dotSet)
    
    #劃分下包右包的點(diǎn)集
    rightSet,minDot=minSize(seq,dot3,dot4,dotSet)
        
    #繪圖檢測(cè)---------------------------------------------------------------
    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()
    #----------------------------------------------------------------------
    
    #對(duì)下包的右包進(jìn)一步劃分
    if rightSet:
        divideDown(rightSet,dot3,minDot,minDot,dot4,dotSet)
    
    return dotSet
    
    
#遞歸主函數(shù)
def mainDivide(seq):
    #將序列中的點(diǎn)按橫坐標(biāo)升序排序
    seq.sort()
    res=[]
    #獲取橫坐標(biāo)做大、最小的點(diǎn)及橫坐標(biāo)中位數(shù)
    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
    #對(duì)序列劃分為直線(xiàn)dot1 dot2左右兩側(cè)的點(diǎn)集并找出兩個(gè)點(diǎn)集的距直線(xiàn)最遠(yuǎn)點(diǎn)
    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)
    #調(diào)用內(nèi)建遞歸函數(shù)
    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()