張青年

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2.1 算法框架

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

算法的主要過程如下:

1)初始化。假定柵格平面大小為M×N，設置4個二維數(shù)組 e［M，N］，a［M，N］，b［M，N］和 s［M，N］。ei，j記錄像元 Pi，j是否以自身為源，aij和 bij分別記錄像元Pi，j的源 (初始源或虛擬源)的列坐標和行坐標，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)以柵格平面的左上角為起點，行號遞增，列號遞增，自上而下順序訪問各像元 Pi，j。若 Pi，j不是障礙物像元，則依據(jù)其左、左上、上和右上側(cè)鄰居像元依次更新其源坐標和距離 aij、bij、ei，j、si，j。不失一般性，記鄰居像元為U。

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

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

然后判斷U的源src(U)的上級源src(src(U))是否可直達目標像元P(圖3b)。如果P到src(src(U))的連線上沒有任何障礙物像元，則src(src(U))可直達P，依據(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))的列坐標和行坐標。

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

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

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

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

2.2 可見性判斷

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

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

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

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

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

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

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

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

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

4)其他情況下繼續(xù)向前進行檢測。

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

2.3 算法效率分析

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

3 實驗結(jié)果及分析

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

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

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

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

圖7 誤差分析Fig.7 Error analysis

4 結(jié)論

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

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

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

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

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

5)本文算法的時間復雜度為O(n)，計算效率高。

［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］.計算機學報，1995，18(8):611－616.

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

［8］任勇勇，潘泉，張紹武，等.基于圍線分層掃描的完全歐氏距離變換算法［J］.中國圖象圖形學報，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］徐達麗，任洪娥，徐海濤，等.基于鏈碼技術的距離變換改進算法［J］.計算機工程與應用，2009，45(25):176－178.

［11］陸宗騏，朱煜.用帶形狀校正的腐蝕膨脹實現(xiàn)Euclidean距離變換［J］.中國圖象圖形學報，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版.武漢:武漢大學出版社，2006.