1. 程式人生 > 實用技巧 >經典分治問題,平面N個點求最近點對

經典分治問題,平面N個點求最近點對

大家好,我們今天來看一道非常非常經典的演算法題——最近點對問題

這個問題經常在各種面試當中出現,難度不低,很少有人能答上來。說實話,我也被問過,因為毫無準備,所以也沒有答上來。是的,這道題有點神奇,沒有準備的人往往答不上來。

題意

我們先來看下題意吧,題意很簡單,在一個平面當中分佈著n個點。現在我們知道這n個點的座標,要求找出這n個點當中距離最近的兩個點的間距。

我不確定這個問題是否出自於天文學,但是把它放到天文的背景當中非常合適。想象一下在浩瀚的宇宙當中,存在著無數的星辰,我們想要找到其中距離最近的兩顆天體。它們有可能是雙子星,也有可能是伴星系……這麼想想,有沒有覺得很浪漫?

我們來分析一下問題,會發現一個矛盾之處。矛盾的地方在於如果我們要求出每兩個點之間的距離,那麼複雜度一定是

,因為n個點取兩個點一個有種可能。如果存在更快的演算法,那麼勢必我們不能求出所有點對之間的距離,但如果我們連所有的距離都沒有列舉過,如何可以判斷我們找到的一定是對的呢?

我當時在面試的時候就是這麼回答的,雖然我們現在知道這個說法是錯的,但是如果沒有這一層資訊,你還能判斷嗎?

分治法

如果我們仔細思考一下,會發現這個問題和排序其實非常類似。因為我們在排序的時候,表面上來看每兩個點之間都存在大小關係,我們要排序似乎也要獲得這些關係。但實際上,我們都知道,無論是快速排序還是歸併排序都可以做到的時間內完成排序。

無論是快速排序還是歸併排序,本質上都是利用的分治法。那麼這道題是否也可以使用分治法求解呢?

答案當然是可以的,既然是使用分治法,那麼我們首先要做的就是拆分,將整個的資料拆成兩個部分,使用遞迴分別完成兩個部分,然後再合併得到完整的結果。在這個問題當中,我們要拆分資料非常簡單,只需要按照x軸座標對所有點進行排序,然後選擇中點進行分割即可,分割之後我們得到的結果如下:

拆分結束之後,我們只需要分別統計左邊部分的最近點對、右邊部分的最近點對,以及一個點在左邊一個點在右邊的最近點對即可。對於前面兩種情況都很好解決,我們只需要遞迴就可以搞定了,但對於第三種情況應該怎麼辦?這也是本題的難點所在。

要分析清楚這個問題不是非常容易,需要深入的思考,首先我們通過遞迴呼叫可以獲得左邊部分SL的最短距離D1以及右邊部分SR的最短距離D2。我們取

,也就是左右兩邊最短距離的最小值,這個應該很好理解。

求出了D之後,我們就可以用它來限定一個點在SL一個點在SR這種情況的點對的範圍了,不然的話我們要比較兩邊各有n/2個點的情況,依然計算複雜度很大。

我們來分析一下問題,我們在左側隨便選擇一個點p,我們來想一個問題,對於點p而言,SR一側所有的點都有可能與它構成最近點對嗎?當然不是,有一些離得遠的是明顯不可能的,對於這些點我們沒有必要一一遍歷,直接都可以批量忽略。要想和p點構成最近點對,必須在下圖這個虛線框起來的範圍內

這個虛線構成的框是一個長方形,它的寬是D,長是2D。這是怎麼來的呢?其實很簡單,對於p點來說,要想和他構成全域性的最近點對,那麼距離它的距離一定要小於目前的最優解D。既然距離要小於D,那麼意味著它們的橫縱座標之差的絕對值必須也要小於D。

當然這個框只是我們直觀看到的,在實際演算法執行的時候是沒有這個框的,我們只能根據座標軸自己去進行判斷某個點在不在框裡。

有了這個框之後,我們產生了另外一個問題,那就是這個框到底可以起到多大的作用呢?有了這個框就可以降低演算法複雜度了嗎?會不會出現右側所有點都在框裡的極端情況呢?

其實我們只需要簡單分析一下就會發現這是不可能的,不僅可以判斷出這是不可能的,而且我們還可以得出另外一個非常非常驚人的結論。

首先,我們來論證一下為什麼右側所有的點都落在這個虛線框裡是不可能的。我們先來看最極端的情況,最極端的情況就是我們選中的p點就在分割線上。那麼以它畫出來的框應該全部都落在SR區域,畫成圖大概是這樣的:

但是我們簡單想一下會發現一個問題,就是這個虛線框裡的點的數量不可能是無限的。因為對於框裡的點我們有一個基本的要求,就是在這個框裡並且在SR區域內的點,兩兩之間的距離不得小於D。如果小於D了就和我們剛才得到D是左右兩側最小距離的結論矛盾了。那麼上面圖中的情況其實是不可能的,因為這麼多點聚集在一起明視訊記憶體在兩個點的距離小於D了,這就矛盾了。

也就是說由於存在這個距離的限制,能夠落在這個虛線框裡的點的數量是有限的,而且這個數量比大家想的也許要小得多,有多小呢?小到最多隻有6個,也就是下面這種情況:

在上圖當中,一共有6個點,這6個點兩兩之間的最短距離是D,這是最極端的情況。無論我們如何往其中加入點,都一定會產生兩個點之間的距離小於D。這是我們很直觀的感受,有沒有辦法證明呢?其實也是有的,我們可以把這個矩形進行六等分變成下圖這樣:

我們來分析一下,上圖的每一個小矩形的長是,寬是,它的對角線長度是。那麼根據鴿籠原理,如果我們放入超過6個點,必然會存在一個小矩形記憶體在兩個點。而小矩形內最大的距離小於D,也就是說這兩個點的距離必然也小於D,這就和我們之前的假設矛盾了,所以可以得出超過7個點的情況是不存在的。

也就是說對於SL側的點p,我們在SR側最多隻能找出6個點來可能構成最短點對,這樣我們需要篩查的點對數量就大大減小。並且對於SL側的點來說,並不是所有的點都需要考慮的,只有和中點O橫座標差值小於D的點才需要考慮

表面上看起來我們所有的分析都結束了,但實際上還有一個問題沒有解決。就是我們怎麼樣找到這6個點呢?顯然只根據橫座標是不行的,這個時候就需要考慮縱座標了。我們將點集分成左右兩個部分之後,對右側部分按照縱座標進行排序,對於左側的點(x, y)而言,我們只需要篩選出右側滿足縱座標範圍在(y - d, y + d)範圍內的點,這樣的點最多隻有6個。我們可以利用二分法找到縱座標大於 y - d的最小的點,然後依次列舉之後的6個點即可。

程式碼實現

在我們實現演算法之前,我們需要先生成測試資料,否則如何驗證我們的演算法是否有問題呢?而且這個演算法也是我從頭開發的,對於debug也有幫助。

在這道題當中,測試資料還是比較簡單的,只需要生成兩個隨機數作為座標即可。我們呼叫這個方法先生成200個點作為測試。

importrandom


defrandom_point():
x,y=random.uniform(0,1000),random.uniform(0,1000)
return(x,y)

points=[random_point()for_inrange(200)]

接著我們再實現暴力解法,用來檢測我們的演算法的正確性,這一段我想應該不用我多說,大家都能理解。

defdistance(x,y):
returnmath.sqrt((x[0]-y[0])**2+(x[1]-y[1])**2)


defbrute_force(points):
ret=sys.maxsize
a,b=None,None
n=len(points)
foriinrange(n):
forjinrange(i+1,n):
dis=distance(points[i],points[j])
ifdis<ret:
ret=dis
a,b=i,j
returnret,points[a],points[b]

最後是重頭戲了,其實演算法本身的程式碼量並不大,但是其中的細節不少,稍有不慎就可能翻車。

defdivide_algorithm(points):
n=len(points)
#特判只有一個點或者是兩個點的情況
ifn<2:
returnsys.maxsize,None,None
elifn==2:
returndistance(points[0],points[1]),points[0],points[1]

#對所有點按照橫座標進行排序
points=sorted(points)
half=(n-1)//2
#遞迴,這裡有一個問題,為什麼要先排序再遞迴?
d1,a1,b1=divide_algorithm(points[:half])
d2,a2,b2=divide_algorithm(points[half:])
d,a,b=(d1,a1,b1)ifd1<d2else(d2,a2,b2)

calibration=points[half][0]

#根據中間的位置將點集分成兩個部分
left,right=[],[]
foruinpoints:
ifcalibration-d<u[0]<calibration:
left.append(u)
elifcalibration<=u[0]<calibration+d:
right.append(u)

#右側點集按照縱座標排序
right=sorted(right,key=lambdax:(x[1],x[0]))

res=d

foruinleft:
#左開右閉的二分
l,r=-1,len(right)-1
whiler-l>1:
m=(l+r)>>1
ifright[m][1]<=u[1]-d:
l=m
else:
r=m

idx=r

#在範圍內最多隻有6個點能夠構成最近點對
forjinrange(7):
ifj+idx>=len(right):
break
ifdistance(u,right[idx+j])<res:
res=distance(u,right[idx+j])
a,b=u,right[idx+j]

returnres,a,b

演算法是實現完了,但是仍然有一些細節,比如說為什麼我們在分治的時候需要先排序再遞迴呢?直接分成兩個部分遞迴行不行?為什麼不行?比如我們二分的時候使用的是左閉右開的區間二分?

這兩個問題我先不給出答案,希望大家能夠自己嘗試著思考一下。如果有想法的話,歡迎在下方給我留言討論。

今天的文章就到這裡,衷心祝願大家每天都有所收穫。如果還喜歡今天的內容的話,請來一個三連支援吧~(點贊、關注、轉發

原文連結,求個關注