以障礙物凸角點為中繼源的歐氏距離變換算法*

張青年

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

距離變換是一種計算背景像元的距離的過程，其計算結(jié)果是一幅距離圖像，標識各個背景像元到最近的特征像元(源)的距離[1]。常用的距離有棋盤距離、街區(qū)距離、倒角距離和歐氏距離。其中，棋盤距離、街區(qū)距離、倒角距離等非歐氏距離是對歐氏距離的一種近似。已有的歐氏距離變換算法在提高距離變換的精確性、簡單性和時間效率等方面進行了大量研究[2-12]。但這些距離變換算法通常都沒有考慮障礙物的阻隔作用，所得距離為兩點間的直線距離。在存在河流、山丘等障礙物的地理空間中，兩點間的實際可通行路徑需要繞過障礙物，其路徑長度大于直線距離。因此，普通距離變換得到的距離值與實際可通行距離存在一定差距。

如果一個地理空間中存在若干個障礙物，該地理空間是一個非凸域。非凸域中兩點之間的距離被稱為“測地距離(geodesic distance)”，也稱為“大地距離”[13]。如圖1所示，面目標A中有一個洞。對于面目標A而言，A中任意兩個相異點P1和P2之間的測地距離為包含于A中的該兩點間最短路徑的長度。相應(yīng)地，非凸域中的距離變換被稱為“測地距離變換(geodesic distance transform)”[14]。

目前已有少數(shù)學者研究了測地距離變換算法[14-18]。其中，文獻[14-17]采用了以源為中心的圈層傳播算法，文獻[18]采用了光柵掃描算法。文獻[17-18]都引入了中繼源來計算測地距離。但由于中繼源并不在最短距離傳播路徑上，因而存在較大的計算誤差。

圖1 測地距離Fig.1 Geodesic distance

本文將障礙物凸角點作為中繼源，采用蠻力算法進行測地距離變換。由于障礙物凸角點位于計算測地距離的最短路徑上，因此距離偏移量為0。本文首先討論測地距離變換中的中繼源選擇問題，然后利用障礙物凸角點設(shè)計了新的測地距離變換算法，再后分析了新算法的實驗結(jié)果。

1 距離傳播路徑上的中繼源

距離變換可以看作是一個將光線從源逐步傳遞到周圍的背景像元并計算源到背景像元的距離的過程。為簡便起見，本文將這一過程稱為距離傳播。在障礙空間中，從源出發(fā)的光線必須繞過障礙物才能到達位于障礙物背后的像元。在前人的研究中，都是將障礙物附近的某個背景像元作為中繼源，從源出發(fā)沿直線路徑到達中繼源，再經(jīng)其他中繼源沿直線路徑到達目標像元。例如，文獻[17]將障礙物背后的某個背景像元設(shè)置為中繼源，稱為局部最近的不可見像元(Locally Nearest Hidden Pixel, LNHP)。LNHP需同時滿足5個條件：屬于背景像元、有一個8-鄰居為障礙物像元、被障礙物遮擋、有一個8-鄰居為背景像元，到源的距離小于8個鄰居像元的距離。

類似地，以障礙物前方或后方的其他背景像元為中繼源也會存在計算誤差[18]。產(chǎn)生距離誤差的原因在于從源到障礙物后方的背景像元的最短距離傳播路徑不經(jīng)過作為中繼源的背景像元的中心點。

實際上，最短距離傳播路徑經(jīng)過障礙物的凸角點，在凸角點處發(fā)生轉(zhuǎn)折。例如，在圖2中，從源O到背景像元Q的最短可通行路徑在障礙物像元B的角點0和角點3處發(fā)生轉(zhuǎn)折。障礙物像元B位于一個條帶狀障礙物的端點處，而其角點0和3是像元B上的兩個凸角。因此，本文以障礙物的凸角點作為中繼源，使距離傳播路徑緊貼障礙物繞轉(zhuǎn)到障礙物背后的不可見像元。

圖2 LNHP和距離誤差Fig.2 LNHP and distance error

2 測地距離變換新算法

新算法將障礙物凸角點作為中繼源，在可見性判斷基礎(chǔ)上采用蠻力算法計算障礙物凸角點和各個像元到源的測地距離。

2.1 障礙物凸角點的識別

一種方法是將障礙物像元轉(zhuǎn)換為矢量圖形，再利用多邊形頂點凸凹性識別算法提取障礙物凸角點?？蓪⒄系K物像元表示成矢量單元格，再將相鄰的矢量單元格合并成一個障礙物多邊形，即可分析多邊形邊界上各個頂點的凸凹性，從而提取障礙物凸角點。

本文不進行柵-矢轉(zhuǎn)換，而是直接依據(jù)障礙物像元的鄰接關(guān)系提取障礙物的凸角點。顯然，障礙物的凸角點位于障礙物邊緣，其所在的障礙物像元必須與背景像元(或源)相鄰。進一步分析，凸角點位于障礙物邊緣的凸起處，也就是該角點位于1個障礙物像元和3個背景像元的交點處。設(shè)存儲圖像的數(shù)組為P[M, N]，像元的左上角、右上角、右下角和左下角分別編號為0、1、2和3，對于任意像元P[i,j]， 采用以下準則進行判斷：

if (P[i,j]是障礙物)

{

if(P[i,j-1]、P[i-1,j-1]和P[i-1,j]都不是障礙物)

{P[i,j]的角點0為凸角點; }

if(P[i-1,j]、P[i-1,j+1] 和P[i,j+1]都不是障礙物)

{P[i,j]的角點1為凸角點; }

if(P[i,j+1]、P[i+1,j+1]和P[i+1,j]都不是障礙物)

{P[i,j]的角點2為凸角點; }

if(P[i,j-1]、P[i+1,j-1]和P[i+1,j]都不是障礙物)

{P[i,j]的角點3為凸角點; }

}

如圖3所示，像元A的角點0和3、像元D的角點3、像元E的角點2、像元I的角點1和2將被識別為障礙物凸角點。

圖3 障礙物的凸角點Fig.3 Convex corners on obstacles

2.2 可見性判斷方法

在距離傳播計算過程中，必須檢測源(或中繼源)相對于目標像元是否可見。若源(或中繼源)可見，則沿直線路徑傳播距離到目標像元；否則，需引入新的中繼源，距離傳播路徑在障礙物附近發(fā)生轉(zhuǎn)折，并經(jīng)過中繼源繼續(xù)傳播距離到目標像元。

可見性判斷算法有文獻[16]基于障礙物的角度排序的測試算法，文獻[17]基于LNHP的障礙物遮蔽角的測試算法，以及文獻[18]基于直線全路徑柵格化的公共源檢測方法。由于文獻[18]方法的檢測對象可以是像元中心點、角點和其他任意點的點對，本文采用該方法判斷可見性。其算法原理是從目標像元P開始向源src反向柵格化，檢測P和src的連線是否經(jīng)過障礙物像元。如果沿直線檢測過程中遇到障礙物像元，則src對于P而言不可見；否則，沒有遇到障礙物，則判定像元src相對于P可見。該算法并不需要檢測P和src連線上的所有像元，其時間復雜度為常數(shù)階[18]。如圖4(a)所示，從背景像元P開始朝著源src方向檢測，在第一個轉(zhuǎn)折處的兩個像元的源都是src，則可判定src對于P可見，退出檢測過程。如圖4(b)所示，如果第一個轉(zhuǎn)折處的兩個像元不是障礙物但其源不都是src，則繼續(xù)向前檢測，沒有遇到障礙物，并且在第二個轉(zhuǎn)折處的兩個像元的源都是src，可判定src對于P可見，退出檢測過程。

圖4 判斷目標像元的可見性Fig.4 Check visibility of target pixel

2.3 距離計算步驟

首先計算各個中繼源的測地距離，再逐個背景像元計算其距離。

1)遍歷中繼源序列，計算其測地距離。

構(gòu)造一個鏈表srcList，其元素為結(jié)構(gòu)類型sturct{ IPoint pt; double dist; }。插入所有源和中繼源，以下統(tǒng)稱為“源”。其中，源的距離dist為0，中繼源的距離dist初值為極大值D。對srcList進行多輪遍歷，更新源的距離值，直到距離沒有變化為止。算法的基本流程如下：

do{

置flag為false;

依次取srcList的元素srcList[i];

{

依次取srcList的元素srcList[j];

{

計算newDist = srcList[j].dist + srcList[j].pt與srcList[i].pt之間的直線距離;

if (newDist>=srcList[i].dist) continue;

判斷srcList[j]對于srcList[i]是否可見;

如果可見

{

srcList[i].dist = newDist;

置flag為true;

}

}

}

} while (flag)

2)采取光柵掃描方式逐個像元計算其測地距離。算法的基本流程是：

逐行逐列遍歷像元陣列元素P[i, j];

{

依次取srcList的元素srcList[k];

{

if (P[i, j]是障礙物像元) continue;

計算newDist = srcList[k].dist + srcList[k].pt與P[i, j]之間的直線距離;

if (newDist>= P[i, j].dist) continue;

判斷srcList[k]對于P[i, j]是否可見;

如果可見, 置P[i, j].dist = newDist;

}

}

2.4 時間復雜度分析

新算法包括三個過程：識別中繼源、計算中繼源的測地距離、計算背景像元的測地距離。設(shè)像元總個數(shù)為n，障礙物像元個數(shù)為m，中繼源個數(shù)為k。顯然，m<

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

采用一幅256*256大小的圖像進行了實驗。該圖像有1個源，3個寬度為3個像元的障礙物。其中，源的坐標為(35, 175)，位于圖像右上角。位于圖像上半部的障礙物兩端均為鋸齒狀，中部有U形轉(zhuǎn)折。位于圖像左下部的障礙物呈直角形狀，右下部的障礙物是兩個直角的組合物，它們的端點是平直的。

新算法識別了障礙物的20個凸角。其中，上半部、左下部和右下部的障礙物上各有12、3和5個凸角點。圖像上的全部障礙物凸角點都被正確識別，沒有漏識別和錯識別的凸角點。經(jīng)過解析方法計算比較，各個障礙物凸角點的距離計算誤差為0。具體情況見表1所示。

表1 新算法凸角計算準確性

新算法計算的距離圖像如圖5所示。在圖像右上部，等距離線在源附近呈圓環(huán)狀分布。在圖像上部障礙物的兩端，分別形成一個環(huán)狀結(jié)構(gòu)，并在圖像左下部匯合。圖像左下角和右下角距離值最大。整幅圖像的最大距離值為570.131，平均距離值為277.934。

圖5 新算法計算的距離圖像Fig.5 The distance image by the new algorithm

采用解析方法計算參考距離圖，并采用文獻[17]中的GDT-LNHP算法和文獻[18]中的光柵掃描法進行測地距離變換。其中，GDT-LNHP算法是現(xiàn)有精度最高的測地距離變換算法[17]。將新算法與文獻[17-18]算法分別與參考距離圖進行比較，統(tǒng)計距離誤差，見表2所示。已有算法的平均誤差小于0.5個像元，但最大誤差都接近或超過1個像元。新算法的最大誤差和平均誤差均為0?？梢姡滤惴ǖ木嚯x準確性優(yōu)于已有算法。

表2 GDT-LNHP和光柵掃描法的距離計算誤差

文獻[17-18]算法的距離誤差的空間分布見圖6。新算法不存在距離誤差。從圖6可以看出，文獻[17-18]的算法的最大距離誤差方向略有差異。其中，文獻[18]算法在緊靠障礙物背后位置的距離誤差較大。

圖6 GDT-LNHP和光柵掃描算法的誤差分布Fig.6 Absolute differences for GDT-LNHP and raster scan algorithm

圖7 新算法在非凸域內(nèi)的距離變換Fig.7 Distance transform in a non-convex domain by the new algorithm

此外，還采用新算法在一些復雜非凸域中進行了距離變換實驗，部分結(jié)果見圖7。在圖7中，有一個帶分岔的開敞空間，其外圍是封閉區(qū)域。開敞空間為背景區(qū)域，封閉區(qū)域為障礙物。在開敞空間中的兩個分岔附近設(shè)一個源，計算開敞空間中各處到源的測地距離。采用新算法計算的最大距離為497.284，平均距離為230.446，如圖7(a)所示；將距離值按間距10分段，各距離段的可視化結(jié)果如圖7(b)所示。從圖7(a)看，開敞空間中離源越遠的位置的距離越大。從圖7(b)看，距離等值線在直線段開敞空間呈較規(guī)則的環(huán)狀分布，在彎曲段和分岔口則呈現(xiàn)為扭曲的環(huán)狀。該圖的距離值與實際相符。

4 結(jié) 論

根據(jù)上述距離變換研究和實驗，可以得到以下結(jié)論：① 新算法以障礙物凸角點為中繼源，中繼源位于最短距離傳播路徑上，不存在距離偏差，從而避免了中繼源設(shè)置不當引入的附帶誤差。② 新算法計算的距離值準確，其距離誤差是0個像元。距離值的精度僅與像元的大小有關(guān)。③ 新算法的時間復雜度接近于O(n)。其中，線性階n有一個與障礙物凸角點個數(shù)有關(guān)的系數(shù)。當障礙物凸角點的個數(shù)較多時，該系數(shù)較大。④ 新算法可應(yīng)用于二維障礙空間中的距離變換。其距離準確性不受障礙物數(shù)量、位置和形狀的影響。

參考文獻：

[1] ROSENFELD A, PFALTZ J. Sequential operations in digital pictures 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(3): 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(11): 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] 徐達麗,任洪娥,徐海濤,等. 基于鏈碼技術(shù)的距離變換改進算法[J]. 計算機工程與應(yīng)用, 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(2): 252-257.

[13] 郭仁忠. 空間分析[M]. 武漢測繪科技大學出版社， 2000.

[14] LANTUEJOUL C, MAISONNEUVE F. Geodesic methods in quantitative image analysis[J], Pattern Recognition, 1984, 17(2): 177-187.

[15] PIPER J, GRANUM E. Computing distance transformations in convex and non-convex 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(5): 561-570.

[18] 張青年. 一種顧及障礙物的歐氏距離變換方法[J]. 中山大學學報：自然科學版, 2013, 52(1): 130-135.