一種顧及障礙物的歐氏距離變換方法*

張青年

(中山大學(xué)地理科學(xué)與規(guī)劃學(xué)院，廣東廣州 510275）

距離變換是一種二值圖像處理技術(shù)，被廣泛地應(yīng)用于圖像處理和模式識(shí)別等領(lǐng)域。它將二值圖像中的像素區(qū)分為特征 (源)和背景兩種類型，計(jì)算背景像素的距離值而產(chǎn)生一幅距離圖像［1］。在距離變換相關(guān)研究中提出了大量的距離變換算法，這些算法主要集中在提高距離變換的精確性、簡(jiǎn)單性和時(shí)間效率等方面［2－13］。這些算法通常都沒有考慮障礙物的影響，其距離變換是在一個(gè)沒有障礙物的空間中進(jìn)行的。但在實(shí)際的地理空間中，通常有河流、山丘等地物存在，對(duì)距離的傳遞起到阻隔作用，即障礙物兩側(cè)的兩個(gè)點(diǎn)之間的通行距離并不是直線距離。因此，考慮障礙物影響的距離變換算法才可能得到實(shí)際通行距離，該類算法的應(yīng)用領(lǐng)域包括可視域分析、目標(biāo)分割、距離/厚度測(cè)量、路徑規(guī)劃、影響區(qū)域范圍等［13］。目前只有少數(shù)學(xué)者研究了顧及障礙物的距離變換算法［13－17］，而且他們的算法也較復(fù)雜，都采用了由源向外逐個(gè)圈層傳遞距離的方式和較復(fù)雜的像元可見性判斷方法。文獻(xiàn) ［14－15］最早在距離變換中考慮障礙物，采取逐圈層傳遞距離方式。Coeurjolly等［16］提出了一種基于可見性分析的距離變換方法，基于角度排序判斷像元可視性。Cárdenes等［13］提出了一種基于阻隔點(diǎn) (Occlusion points)探測(cè)的距離變換算法，依據(jù)阻隔點(diǎn)相對(duì)于障礙物的最大最小遮蔽角判斷后續(xù)距離傳遞過程中像元的可見性。這些算法復(fù)雜性高，軟件實(shí)現(xiàn)困難。除此之外，ArcGIS軟件也實(shí)現(xiàn)了一種距離變換算法，通過累計(jì)距離傳遞路徑上各個(gè)像元中心點(diǎn)之間距離的和得到距離值。這種算法原理簡(jiǎn)單，但若源像元與目標(biāo)像元不在同一行、同一列或同一條對(duì)角線上，其累計(jì)距離值與二者間的歐氏距離有一定偏差，并且誤差隨傳遞距離的增加而增大。本文從提高算法簡(jiǎn)單性角度出發(fā)，提出了一種顧及障礙物的距離變換算法，基于柵格掃描方式傳遞距離，并采用改良的直線柵格化方式判斷像元可見性，算法簡(jiǎn)單，具有線性時(shí)間復(fù)雜度，所得到的距離準(zhǔn)確度較高。

1 基于柵格掃描的歐氏距離變換原理

以8鄰居距離變換為例，基于柵格掃描的歐氏距離變換基本過程是［2］:首先自上而下掃描柵格平面，分別依據(jù)左、左上、上和右上側(cè)鄰居像元的距離值計(jì)算目標(biāo)像元的距離值，用4個(gè)新距離值和目標(biāo)像元的當(dāng)前距離值中的最小值更新目標(biāo)像元的距離值;然后自下而上掃描柵格平面，依據(jù)右、右下、下和左下側(cè)鄰居像元的距離值計(jì)算目標(biāo)像元的距離值，用4個(gè)新距離值和目標(biāo)像元的當(dāng)前距離值中的最小值更新目標(biāo)像元的距離值。

在距離變換過程中，可以采取累計(jì)距離的方式計(jì)算目標(biāo)像元到源的距離 (圖1)，可表示為:

圖1 柵格掃描與距離增量Fig.1 Distance offset in raster scan algorithm

其中，dij為目標(biāo)像元Pij的距離，duv為鄰居像元Puv的距離，Δd為Pij與Puv之間的距離。目標(biāo)像元與其上、下、左、右4個(gè)鄰居像元之間的距離為1，與左上、右上、左下、右下4個(gè)鄰居像元之間的距離為。

但這種累計(jì)距離計(jì)算方式只在源的上、下、左、右及左上、右上、左下、右下8個(gè)方向上的距離傳遞路徑是直線，得到準(zhǔn)確的直線距離;在其他方向上的距離傳遞路徑都是一條從源到目標(biāo)像元的折線，該折線長(zhǎng)度比二者之間的直線距離大。設(shè)某位置到源的行距和列距的最大值為k，其最大誤差方向的誤差為 0.089k［17］。

在計(jì)算空間中存在障礙物的情況下，只需限定障礙物像元不參與距離比較，該算法仍可有效地進(jìn)行距離變換。但隨著離源的距離增加，其累積誤差不斷增大。這種算法被ArcGIS用于計(jì)算成本距離，適用于小幅圖像的距離計(jì)算。

更準(zhǔn)確距離計(jì)算方法是依據(jù)偏移量計(jì)算像元到源的直線距離。在距離變換過程中，記錄各像元到源的行偏移量和列偏移量，可以視為一個(gè)向量，依據(jù)距離向量計(jì)算直線距離。目標(biāo)像元到源的偏移量依據(jù)其鄰居像元的偏移量來(lái)計(jì)算，即

其中，aij為目標(biāo)像元Pij的距離向量，auv為鄰居像元 Puv的距離向量，為從 Pij到 Puv的距離向量。如圖2所示，目標(biāo)像元到左、左上、上和右上側(cè)鄰居像元的距離向量分別為 (1，0)、 (1，1)、(0，1)和 (－1，1)，目標(biāo)像元到右、右下、下和左下側(cè)鄰居像元的距離向量分別為 (－1，0)、(－1， －1)、(0， －1)和 (1， －1)。采用距離向量方式能夠得到準(zhǔn)確的歐氏距離。

圖2 柵格掃描與距離矢量Fig.2 Distance vector in raster scan algorithm

采用偏移量計(jì)算距離的方式不適合于存在障礙物的情形。在目標(biāo)像元到源的連線上有障礙物時(shí)，源到目標(biāo)像元的距離傳遞路徑應(yīng)為一條繞過障礙物的折線，此時(shí)直線距離所代表的直線傳遞路徑顯然不符合實(shí)際情況。

2 顧及障礙物的歐氏距離變換方法

2.1 算法框架

如前所述，在目標(biāo)像元到特征像元 (源)的連線上存在障礙物時(shí)，目標(biāo)像元被障礙物阻隔，源到目標(biāo)像元的距離傳遞路徑不再是一條直線。本文在基于距離矢量的柵格掃描算法中引入“虛擬源”，允許距離傳遞路徑發(fā)生轉(zhuǎn)折，在被障礙物阻隔的區(qū)域中繼續(xù)進(jìn)行距離傳遞。虛擬源是被障礙物阻隔的一個(gè)非障礙物像元。為區(qū)別起見，將原有特征像元稱為初始源。

算法的主要過程如下:

1)初始化。假定柵格平面大小為M×N，設(shè)置4個(gè)二維數(shù)組 e［M，N］，a［M，N］，b［M，N］和 s［M，N］。ei，j記錄像元 Pi，j是否以自身為源，aij和 bij分別記錄像元Pi，j的源 (初始源或虛擬源)的列坐標(biāo)和行坐標(biāo)，si，j記錄像元 Pi，j到初始源的距離，1 ≤i≤M，1≤j≤N。如果像元Pi，j是初始源，令aij=i，bij=j，ei，j=1，si，j=0;否則令 ei，j=0，si，j=2×max(M，N)。

2)以柵格平面的左上角為起點(diǎn)，行號(hào)遞增，列號(hào)遞增，自上而下順序訪問各像元 Pi，j。若 Pi，j不是障礙物像元，則依據(jù)其左、左上、上和右上側(cè)鄰居像元依次更新其源坐標(biāo)和距離 aij、bij、ei，j、si，j。不失一般性，記鄰居像元為U。

首先判斷U是否為障礙物像元，若是則直接跳過，不參與距離比較和更新。

圖3 依據(jù)鄰居像元U更新目標(biāo)像元P的距離Fig.3 Update distance of P according to neighbor U

然后判斷U的源src(U)的上級(jí)源src(src(U))是否可直達(dá)目標(biāo)像元P(圖3b)。如果P到src(src(U))的連線上沒有任何障礙物像元，則src(src(U))可直達(dá)P，依據(jù)它們的行列坐標(biāo)計(jì)算二者之間的直線距離snew。若snew小于P的距離值 sij，則以 src(src(U))為源，更新 eij，aij，bij，sij。令 eij=0，aij=srcx(src(U))，bij=srcy(src(U))，sij=snew。其中srcx(src(U))和srcy(src(U))分別表示src(src(U))的列坐標(biāo)和行坐標(biāo)。

若上級(jí)源src(src(U))不能直達(dá)P，則判斷src(U)是否可直達(dá)P(圖3a)。若src(U)可直達(dá)P，依據(jù)它們的行列坐標(biāo)計(jì)算二者之間的直線距離snew。若snew小于 P的距離值 sij，則以 src(U)為源，更新 eij，aij，bij，sij。令 eij=0，aij=srcx(U)，bij=srcy(U)，sij=snew。

若src(U)不能直達(dá)P(圖3a)，則計(jì)算經(jīng)過U到達(dá)P的距離snew。令snew=s(U)+Δd。其中，Δd為P與U之間的直線距離。若snew小于P的距離值sij，則以src(U)為上級(jí)源，更新 eij，aij，bij，sij。令 eij=1，aij=i，bij=j，sij=snew。即引入 P 為虛擬源。

3)以柵格平面的右下角為起點(diǎn)，行號(hào)遞減，列號(hào)遞減，自下而上順序訪問各像元 Pi，j。若 Pi，j不是障礙物像元，則依據(jù)其右、右下、下和左下側(cè)鄰居像元依次更新ei，j、它的源的坐標(biāo)aij和bij及距離 si，j。

4)重復(fù)執(zhí)行步驟2)和3)若干次，直到距離沒有變化為止。

2.2 可見性判斷

上述算法的一個(gè)核心過程就是可見性判斷，即源到障礙物之間是否有障礙物。Coeurjolly et al.［16］基于障礙物的角度排序進(jìn)行可見性測(cè)試;Cárdenes et al.［13］則記錄每個(gè)LNHP對(duì)應(yīng)的障礙物的遮蔽角范圍并將落在遮蔽角范圍內(nèi)的像元判定為不可見。本文基于直線柵格化算法直接判斷目標(biāo)的可見性。

直線柵格化有八方向柵格化、恒密度柵格化、全路徑柵格化等多種算法。其中八方向柵格化得到的直線最細(xì)，每行或每列只有1個(gè)像元被涂黑;恒密度柵格化得到的直線粗細(xì)均勻，相鄰行或列之間有程度相似的重復(fù);全路徑柵格化則將直線經(jīng)過的每個(gè)像元都涂黑。

圖4 全路徑柵格化Fig.4 Full path rasterization

顯然，可以依據(jù)全路徑柵格化算法判定目標(biāo)像元的可見性。即，在全路徑柵格化過程中，從源開始檢測(cè)位于路徑上的每一個(gè)像元，只要有一個(gè)像元是障礙物像元，則目標(biāo)像元不可見;否則繼續(xù)檢測(cè)下一個(gè)像元，如果路徑上的所有像元都不是障礙物像元，則目標(biāo)像元位于源的可見區(qū)域內(nèi)。

圖5 確定像元P的可見性Fig.5 Check the visibility of pixel P

本文基于全路徑柵格化原理提出一種高效的可見性判斷算法，不需要檢測(cè)直線經(jīng)過的每一個(gè)像元是否為障礙物。其算法原理是從目標(biāo)像元P開始向源src反向柵格化，判斷最初兩行 (或列)像元的可見性即可。設(shè)該直線的跨度為m行n列，如果m＜n，逐列檢測(cè)像元是否需要涂黑，否則逐行檢測(cè)像元是否需要涂黑。如果需要涂黑的多個(gè)像元的源就是src，則可直接判定目標(biāo)像元P在src的可見域內(nèi)，而不必遍歷直線上的所有像元并檢測(cè)其是否為障礙物。以圖5為例，從P到src的連線跨6行13列，從P開始對(duì)連線上的像元進(jìn)行檢測(cè)，每個(gè)列最多檢測(cè)2個(gè)像元。具體算法過程如下:

1)如果連線上某個(gè)列有2個(gè)像元并且它們的源都是src，則 P在src的可見域內(nèi) (如圖5a所示)，退出檢測(cè)過程;

2)如果連線上某個(gè)列只有1個(gè)像元并且它的源是src，繼續(xù)檢測(cè)連線上的下一列的所有像元的源是否為src，如果下一列全部像元的源都是src，則P在src的可見域內(nèi) (如圖5b所示)，退出檢測(cè)過程;

3)如果連線上遇到一個(gè)障礙物像元，則P不在src的可見域內(nèi)，退出檢測(cè)過程;

4)其他情況下繼續(xù)向前進(jìn)行檢測(cè)。

可以預(yù)見，如果P在src的可見域內(nèi)，其連線上靠近P處一定有源為src的像元，因此檢測(cè)過程很快就會(huì)結(jié)束。即使連線長(zhǎng)、經(jīng)過的像元多，也只需檢測(cè)靠近P的數(shù)個(gè)像元。

2.3 算法效率分析

分析上述算法框架可以發(fā)現(xiàn)，其基本過程是柵格掃描，每個(gè)像元被掃描2次，每次與4個(gè)鄰居像元對(duì)應(yīng)的距離值進(jìn)行比較。設(shè)像元個(gè)數(shù)為n，則距離比較和賦值的次數(shù)為8n。參與比較的距離值是依據(jù)鄰居像元計(jì)算出來(lái)的，在距離計(jì)算過程中要執(zhí)行可見性判斷。本文采取直線柵格化算法判定可見性，它是一個(gè)雙重循環(huán)算法，但實(shí)際上并不需要檢測(cè)連線上的所有像元是否為障礙物，而只需檢測(cè)靠近目標(biāo)像元一端的有限個(gè)像元，因此該可見性判斷算法是常量階的。另外，在目標(biāo)像元被障礙物遮擋層次較深的情況下，掃描過程需要重復(fù)多次才能完成距離傳播，但重復(fù)的次數(shù)可以視為一個(gè)常數(shù)。因此，該距離變換算法總體上的時(shí)間復(fù)雜度為O(n)。

3 實(shí)驗(yàn)結(jié)果及分析

為了驗(yàn)證上述方法的可行性，本文構(gòu)造了一幅255行318列的圖像進(jìn)行了距離變換實(shí)驗(yàn)，計(jì)算結(jié)果見圖6a。圖6a中有3個(gè)源，其中O1和O2為點(diǎn)狀源，O3為面狀源;有3個(gè)面狀障礙物B1、B2和B3。從圖6a可以看出，從源向外，距離逐漸增大，距離等值線之間的間距相等 (5個(gè)柵格距離)。點(diǎn)源周圍的距離等值線渾圓，呈現(xiàn)為同心圓結(jié)構(gòu);面源周圍的等值線也是一種圈層結(jié)構(gòu)，其延伸趨勢(shì)與面源輪廓線一致。距離等值線圈層結(jié)構(gòu)在遇到障礙物B1和B2之后遭到破壞，在障礙物背后繼續(xù)進(jìn)行距離傳播，并在障礙物端點(diǎn)處形成新的圈層結(jié)構(gòu)(例如在B1背后靠近兩端處)。在遇到第二重障礙物B3之后，圈層結(jié)構(gòu)再次被破壞，在B3背后靠近兩端處形成新的圈層結(jié)構(gòu)。

由于在障礙空間中逐像元計(jì)算距離真值難度很大，無(wú)法將本文算法計(jì)算結(jié)果與距離真值進(jìn)行對(duì)比。這里改用與ArcGIS軟件計(jì)算結(jié)果比較的方式進(jìn)行驗(yàn)證。目前顧及障礙物的距離變換算法很少且算法復(fù)雜，在實(shí)際工作中通常利用ArcGIS軟件進(jìn)行距離變換，因此ArcGIS軟件計(jì)算結(jié)果具有一定代表性。從前述分析可知，ArcGIS累計(jì)距離算法中的距離傳遞路徑和本文算法中的距離傳遞路徑都繞過了障礙物，計(jì)算出的距離值都大于實(shí)際距離值。因此，計(jì)算出的距離較小的算法更有效 (其結(jié)果更接近距離真值)。將ArcGIS計(jì)算的距離與本文算法計(jì)算的距離相減，得到距離差值圖，見圖6b。可以發(fā)現(xiàn)，本文算法計(jì)算的距離值更小。據(jù)第1節(jié)歐氏距離變換原理分析可知，累計(jì)距離在源周圍的上、下、左、右、左上、左下、右上和右下共8個(gè)方向上是沒有誤差的，而在其他方向上都有距離誤差，而且在這8個(gè)方向中的相鄰兩個(gè)方向的角平分線上誤差最大。在圖6b中，兩種計(jì)算結(jié)果在點(diǎn)狀源O1和O2周圍的“米”字形方向上距離差為0，在面狀源O3周圍的疊加“米”字形方向上距離差也為0，表明本文算法在源周圍的8個(gè)方向上沒有誤差。在源周圍的“米”字形圖案的相鄰兩條分岔的角平分線方向上，ArcGIS計(jì)算的距離值相對(duì)于本文算法的距離值增加最多，這與ArcGIS距離累計(jì)算法的誤差在各方向上的分布一致，表明本文算法在8方向之外其他方向上的距離值也接近于真值。在障礙物背后，距離值最小的像元作為新的“源”繼續(xù)傳遞距離，兩種算法的距離差在障礙物背后形成新的不完整的“米”字形圖案。由于兩種算法在新的“源”位置上存在距離差，“源”周圍的8個(gè)方向上也存在與“源”位置上相等的距離差，其他方向上的距離差更大。此外，兩種算法計(jì)算結(jié)果的距離差隨著與源的距離的增加而變大，兩種算法的最大距離差為15.20。其原因是，ArcGIS距離累計(jì)算法的誤差與到源的距離正相關(guān)，而本文算法的誤差與空間距離無(wú)關(guān)。

圖6 基于柵格掃描的歐氏距離變換實(shí)驗(yàn)Fig.6 Case study of raster scan based euclidean distance transformation

雖然ArcGIS計(jì)算結(jié)果不是距離真值圖，無(wú)法分析本文算法計(jì)算結(jié)果的最大誤差，但我們可以依據(jù)算法原理對(duì)最大誤差作一個(gè)初步分析，見圖7。依據(jù)本文算法，源s在傳遞距離到p時(shí)遇到障礙物，經(jīng)過虛擬源s'之后繼續(xù)傳遞距離到q。這條路徑上，障礙物前方的sp段和后方的s'q段為直線距離，都沒有誤差。但繞過障礙物的ps'q段沒有緊貼障礙物的邊沿uv(uv為障礙物底邊的兩個(gè)角)，有一定的誤差。這里分spu和us'q兩段來(lái)分析。當(dāng)s離p很遠(yuǎn)時(shí)，折線spu與直線su幾乎重合，其長(zhǎng)度差可忽略不計(jì);反之，則有一定差距，當(dāng)s和p為對(duì)角線相鄰時(shí)兩者的長(zhǎng)度差最大，為1+－=0.126。類似地，如果s'q位于uv的同一側(cè) (例如目標(biāo)像元為 q')，則 us'q'的長(zhǎng)度小于uvq'，應(yīng)取us'q'長(zhǎng)度作為真距離，此時(shí)us'q'沒有誤差;反之，當(dāng)s'和q位于uv的兩側(cè)時(shí)，us'q的長(zhǎng)度大于uvq，存在一定誤差，并且q與s'對(duì)角線相鄰時(shí)，誤差達(dá)到最大，為 2－1=0.414。因此，障礙物兩側(cè)的最大誤差和為0.540。將障礙物遮擋的區(qū)域稱為陰影區(qū)，陰影區(qū)內(nèi)最大誤差為0.540，非陰影區(qū)內(nèi)無(wú)誤差。若距離傳遞路徑經(jīng)過n個(gè)障礙物，第n個(gè)障礙物背后的陰影區(qū)內(nèi)的最大誤差不超過0.540 n。

圖7 誤差分析Fig.7 Error analysis

4 結(jié)論

通過以上分析和實(shí)驗(yàn)，可以得出以下結(jié)論:

1)在距離傳播過程中進(jìn)行可見性檢測(cè)，并在目標(biāo)不可見時(shí)引入“虛擬源”，確保了在障礙物背后繼續(xù)進(jìn)行有效的距離傳播，所得距離是從初始源出發(fā)繞過障礙物的折線距離。

2)在可見性檢測(cè)確定上級(jí)源可見的情況下，直接以當(dāng)前的源的上級(jí)源為中繼向外進(jìn)行距離傳播，從而避免了傳播路徑經(jīng)過當(dāng)前的源而發(fā)生不必要的轉(zhuǎn)折，使所得距離值更接近于最短可通行距離。

3)本文算法的距離準(zhǔn)確性較高，優(yōu)于ArcGIS距離累計(jì)算法，在距離傳遞路徑經(jīng)過的第n個(gè)障礙物背后的陰影區(qū)內(nèi)最大誤差為0.540 n，非陰影區(qū)內(nèi)無(wú)誤差。

4)本文算法簡(jiǎn)單，不涉及桶排序等任何復(fù)雜的數(shù)據(jù)結(jié)構(gòu)，易于理解和實(shí)現(xiàn)，適合于點(diǎn)、線、面3種形態(tài)的源和障礙物的距離變換。

5)本文算法的時(shí)間復(fù)雜度為O(n)，計(jì)算效率高。

［1］ROSENFELD A，PFALTZ J.Sequential operations in digital picutures processing［J］.Journal of the ACM，1966，13(4):471－494.

［2］BORGEFORS G.Distance transformations in digital images［J］.Computer Vision，Graphics and Image Processing，1986，34:344 －371.

［3］PAGLIERONI D W.A unified distance transform algorithm and architecture［J］.Machine Vision and Applications，1992，5(1):47－55.

［4］SAITO T，TORIWAKI J.New algorithms for euclidean distance transformation of an n-dimensional digitized picture with applications［J］.Pattern Recognition，1994，27:1551－1565.

［5］FABBRI R，COSTA L D F，TORELLI J C，et al.2D Euclidean distance transform algorithms:A comparative survey［J］.ACM Computing Surveys，2008，40(1):1－44.

［6］陳崚.完全歐幾里德距離變換的最優(yōu)算法［J］.計(jì)算機(jī)學(xué)報(bào)，1995，18(8):611－616.

［7］王鉦旋，李文輝，龐云階.基于圍線追蹤的完全歐氏距離變換算法［J］.計(jì)算機(jī)學(xué)報(bào)，1998，21(3):217－222.

［8］任勇勇，潘泉，張紹武，等.基于圍線分層掃描的完全歐氏距離變換算法［J］.中國(guó)圖象圖形學(xué)報(bào)，2011，16(1):32－36.

［9］LUCET Y.New sequential exact Euclidean distance transform algorithms based on convex analysis［J］.Image and Vision Computing，2009，27(2):37－44.

［10］徐達(dá)麗，任洪娥，徐海濤，等.基于鏈碼技術(shù)的距離變換改進(jìn)算法［J］.計(jì)算機(jī)工程與應(yīng)用，2009，45(25):176－178.

［11］陸宗騏，朱煜.用帶形狀校正的腐蝕膨脹實(shí)現(xiàn)Euclidean距離變換［J］.中國(guó)圖象圖形學(xué)報(bào)，2010，15(2):294－300.

［12］GUSTAVSON S，STRAND R.Anti-aliased Euclidean distance transform［J］.Pattern Recognition Letters，2011，32:252－257.

［13］CáRDENES R，ALBEROLA-LóPEZ C，RUIZ-ALZOLA J.Fast and accurate geodesic distance transform by ordered propagation［J］.Image and Vision Computing，2010，28:307－316.

［14］LANTUEJOUL C，MAISONNEUVE F.Geodesic methods in quantitative image analysis［J］.Pattern Recognition，1984，17:177 －187.

［15］PIPER J，GRANUM E.Computing distance transformations in convex and nonconvex domains［J］.Pattern Recognition，1987，20(6):599 －615.

［16］COEURJOLLY D，MIGUET S，TOUGNE L.2D and 3D visibility in discrete geometry:an application to discrete geodesic paths［J］.Pattern Recognition Letters，2004，25:561－570.

［17］胡鵬，游鏈，楊傳勇，等.地圖代數(shù)［M］.2版.武漢:武漢大學(xué)出版社，2006.