論文《A General algorithm for computing distance transforms in linear time》閲讀筆記整理。

Introduction

距離變換本質是計算點到點的距離,如果在二維空間中,可以按如下方式描述。

假設$B$為尺寸為$m \times n$點的集合,那麽距離變換所計算的就是$B$中某點$(x,y)$到集合$B$中最近點的距離(就是前景點到最近背景點的距離)。若$B$用二值矩陣$b[:, :]$表示(標記圖),且距離為標準歐式距離,那麽每個點$(x,y)$的距離可以表示為: $$ dt[x,y] = \sqrt{EDT(x,y)} $$

其中,

$$ EDT(x,y) = \min{(i, j : 0 \leq i < m \wedge 0 \leq j < n \wedge b[i,j] : (x -i)^2 + (y-j)^2)} $$

解釋下$\min{(k:P(k):f(k))}$這種表示方式的含義:滿足函數$f(k)$最小值,且取值範圍在$P(k)$中的所有$k$的集合。

除了上面提到的歐式距離,還可以采用其他距離計算方式,例如曼哈頓距離(Manhattan distance)、棋盤距離(Chessboard distance)等。

$$ \begin{array}{l} MDT(x,y) &=& \min{(i, j: 0 \leq i < m \wedge 0 \leq j < n \wedge b[i,j] : |x-i| + |y -j|)}\\[2ex] CDT(x,y) &=& \min{(i,j: 0 \leq i < m \wedge 0 \leq j < n \wedge b[i,j] : \max(|x-i|, |y-j|))} \end{array} $$

本文就介紹在這三種距離函數下,如何在綫性時間内完成一幅圖像的距離變換,并且該算法可以并行計算。

現在先將上面三種距離計算函數進行轉換。

$$ \begin{array}{l} EDT(x,y) &=& \min{(i, j : 0 \leq i < m \wedge 0 \leq j < n \wedge b[i,j] : (x -i)^2 + G(i,y)^2)}\\[2ex] MDT(x,y) &=& \min{(i, j: 0 \leq i < m \wedge 0 \leq j < n \wedge b[i,j] : |x-i| + G(i,y))}\\[2ex] CDT(x,y) &=& \min{(i,j: 0 \leq i < m \wedge 0 \leq j < n \wedge b[i,j] : \max(|x-i|, G(i,y)))} \end{array} $$

其中,$G(i,y) = \min{(j : 0 \leq j < n \wedge b[i,j]: |y-j|)}$,很明顯,這麽做轉換在程序上就能將行列操作分離,這樣該算法的整理框架就可以描述為:第一次全圖按列遍歷,即取$C_x$($x$為固定列號),計算$C_x$中每個點$(x, j), 0\leq j < n$與$C_x \wedge B$的最近距離$G(x,y)$。再進行第二次全圖遍歷,此時按行遍歷,此時再利用$EDT$、$MDT$或$CDT$計算每個點$(x,y)$的最終距離值。

下面按兩次遍歷分別説明。

The first phase

根據上面描述,第一階段就是計算$G(x,y)$,在細化$G(x,y)$可以把其分爲兩個階段:

$$ \begin{array}{l} GT(i,y) &=& \min{(j:0\leq j < y \wedge b[i,j] : y - j)} \\[2ex] GB(i,y) &=& \min{(j:y \leq j < n \wedge b[i,j] : j - y)} \end{array} $$

若$b[i,y]$被標記,則$GT(i,y) = 0$(PS:這裏與個人理解是相反的,個人理解是$b[i,y]$未被標記,$GT(i,y) = 0$);其他情況$GT(i,y) = GT(i, y-1)+1$。而$GB(i,y) = GB(i,y+1) + 1$。可以參考下圖理解。

這樣,每列$G(i,y)$的計算就要分爲兩輪遍歷進行,第一輪從上到下($0 \leq y < n$)(*scan 1*),計算每個點的$GT(i,y)$值;第二輪從下到上($n-2 \geq y \geq 0$)(*scan 2*),計算每個點的$GB(i,y)$值,并於之前計算的$GT(i,y)$值比較取最小值記錄在緩存$g(i,y)$中。這樣第一階段的代碼框架已形成。如下圖。

另外説明,實際實現過程可將$\infty$取爲圖像的寬高和$(m+n)$。

The Second phase

第二躺遍歷全圖是按行遍歷,此時假設$y$為常數,將$G(i,y)$表示為$ g(i)$。這次遍歷就需要計算EDTMDTCDT。先將這三種距離計算方式統一成如下形式。 $$ DT(x,y) = \min{(i:0\leq i < m : f(x,i))} \tag{1} $$

其中,

$$ f(x,i) = \begin{cases} (x-i)^2 + g(i)^2, & \text{EDT} \\ |x-i| + g(i), & \text{MDT}\\ \max{(|x-i|, g(i))}, & \text{CDT} \end{cases} $$

另用$F_i = x \mapsto f(x,i)$表示儅取固定$i$值時關於$x$的函數。例如,對於EDT,$F_i$是一條頂點在$(i, g(i))$処的二次曲綫;對於MDT則是一條V-型折綫;對於CDT則是一條擁有最小值限定的V-型折綫。如下圖所示。

從上圖可以看出,我們要求得的結果就是一簇$F_i,(0\leq i < n)$組成的最低位包絡綫(圖中實綫部分)。這些包絡綫由許多小段曲綫構成,從左到右記爲$s[0], s[1], \cdots s[q]$。

繼續抽象函數$DT(x,y)$,將其中$i$的上限$m$替換爲一個變量$u$。 $$ FL(x,u) = \min{(i:0 \leq i < u:f(x,i))} $$

明顯,$DT(x,y) = FL(x,m)$。

再定義 $$ H(x,u) = \min{(h: 0\leq h < u \wedge \forall{(i:0\leq i < u:f(x,h))}: h)} $$

其含義是,固定列號為$x$時,從$0 \backsim u$列最小值位於第$h$列。這樣,$FL(x,u) = f(x, H(x,u))$,因此$DT(x,y) = f(x, H(x,m))$。這樣就把問題轉換為$H(x,m)$的計算。現在再定義兩個集合$S(u)$和$T(h,u)$。

$$ \begin{align} S(u) &= \{H(x,u)\mid 0\leq x < m\}\\ T(h,u) &=\{x \mid 0 \leq x < m \wedge H(x,u) = h\}\ \text{if}\ 0\leq h < u \end{align} $$

$S(u)$表示在每行中,最小值在$0\backsim u$範圍内哪個位置。而$T(h,u)$表示最小值位置在$h$処的$x$取值有哪些。顯然,$S(u)$是個非空集合,且可改寫爲$S(u) =\{ h \mid T(h,u) \neq \emptyset\}$。

説到這裏,已經輸入太多信息有點混亂了,整理下。前面提到當我們找到最低包絡綫時,就能確定最終結果了,那上面提到的這些函數或集合中,$S(u)$和$T(h,u)$就是我們要求得最低包絡綫的數學表示。其中$S(u)$表示哪個函數$F_i$在$T(h,u)$給出的$x$取值範圍下構成最低包絡綫。這樣第二階段的算法框架也差不多顯現出來。與第一階段類似,要進行兩趟遍歷,第一趟遍歷(*scan 3*)就是爲了確定最低包絡綫。第二趟遍歷(*scan 4*)就是計算最終結果了。下面是第二階段算法實現框架。

下面詳細説説最低包絡綫實際如何獲得。計算最低包絡綫本質是一個不斷迭代更新的過程。假設最低包絡綫由$q+1$個片段構成,$S(u) = {s[0], s[1], \cdots, s[q]}$,其中$s[l]$表示第$l$段是由哪個$F_i$函數確定。儅新的$F_u$輸入進來時,會有下面三種情況,可參考圖片。

a). $F_u$完全位於之前所構成的最低包絡綫之上,此時更新$S(u+1) = S(u)$,$T(u,u+1) = \emptyset$;

b). $F_u$完全位於之前所構成的最低包絡綫之下,此時更新$S(u+1) = {u}$,$T(u,u+1)= [0,m)$,因爲此時之前的最低包絡綫完全被$F_u$取代;

c). $F_u$與之前所構成的最低包絡綫相交,計算交點,添加新的最低包絡綫到$T(u,u+1)$。

假設$t[l]$記錄了每段包絡綫的起始點位置,其中$l = q, q-1, \cdots$。從右至左遍歷最低包絡綫,直到找到滿足$F_u(t[l^{*}]) \geq F_{s[l^{*}]}(t[l^{*}])$的$l^{*}$,取$x^* = t[l^*]$或$x^*$為新交點的位置。如果$l^{*} = q$且$x^{*}\geq m$,則為情況a);如果$l^* < 0$,則為情況b);其他則爲c)。

下面説説情況c)中,$x^*$如何求。這裏引入$Sep(i,u)$函數用於$x^*$的計算,該函數計算得到第一個大於或等於$F_u$和$F_i$交點的整型值。 $$ F_i(x) \leq F_u(x) \Leftrightarrow x \leq Sep(i,u) \tag{2} $$

根據上面信息可知,$x^* = Sep(s[l^*], u)$,這樣最終的問題都集中在函數$Sep(i,u)$的計算,而該函數與我們選擇哪種距離度量方式相關。下節中將詳細説明函數如何設計。

在回過頭看代碼,其並不涉及浮點運算,而且利用兩個數組$s[\ ]$和$t[\ ]$就完成整個更新過程。其中數組$s[\ ]$就是用於記錄第$l$段最低包絡綫是由哪個$F_i$函數確定的;$t[\ ]$用於記錄每段最低包絡綫的位置;$q$用於記錄每次開始搜索的最低包絡綫序號。

(這裏有種情況一直沒想明白,假設已知$F_i$最低包絡綫,現添加$F_u$滿足其完全在$F_i$下方,這時代碼中更新$q = 0$,$s[0] = u$。假設添加$F_{u+1}$時,依然滿足條件$f(t[q], s[q]) > f(t[q], u+1)$,那此時依然更新$q = 0$, $s[0] = u+1$,但是我不理解的是,如果$F_u$和$F_{u+1}$實際是有交點的呢?)

Derivation of the Function Sep

下面説説$Sep$函數,它與你采用哪種距離度量函數相關。先來看最簡單的EDT度量函數。根據公式(2)可知:

$$ \begin{align} F_i(x) \leq F_u(x) &\Leftrightarrow (x-i)^2 + g(i)^2 \leq (x-u)^2 + g(u)^2 \\ &\Leftrightarrow x\leq (u^2 - i^2 + g(u)^2 - g(i)^2)\ \mathbb{div}\ (2(u-i)) \end{align} $$

其中,$div$表示整除。這樣EDT的$Sep$函數為: $$ Sep(i,u) = (u^2 - i^2 + g(u)^2 - g(i)^2)\ \mathbb{div}\ (2(u-i)) $$

接著看較複雜點的$Manhattan$度量函數。這裏需要分情況討論。如下圖。

根據前面分析和函數的性質的可知,$i < u$,并且每條曲綫$F_i$都是在$(i, g(i))$処取到最小值點。那麽基於這兩點結合函數圖像:

  • 儅$F_u$完全在$F_i$上方時,此時$g(u) \geq F_i(u)$,這就得到$g(u) \geq g(i) + u-i$;
  • 儅$F_u$完全在$F_i$下方時,此時$g(i) > F_u(i)$,這就得到$g(i) > g(u) + u-i$;
  • 其他情況,$F_i$與$F_u$相交,且相交位置發生在$F_i=x-i + g(i)$與$F_u = u-x + g(u)$的地方。

根據分析,我們可以定義,儅$F_u$完全在$F_i$上方時,$Sep$函數可以返回$+\infty$;儅$F_u$完全在$F_i$下方時,$Sep$函數可以返回$-\infty$;儅兩者相交,則: $$ Sep(i,u) = (g(u)-g(i) + u+i)\ \mathbb{div}\ 2 $$

綜上,$Manhattan$度量的$Sep$函數如下:

$$ Sep(i,u) = \begin{cases} \infty, & \text{if }g(u) \geq g(i) + u - i\\ -\infty, & \text{if }g(i) > g(u) + u - i\\ (g(u) - g(i) + u + i)\ \mathbb{div}\ 2,& otherwise. \end{cases} $$

最後討論最複雜的$Chessboard$度量函數,其也需要分情況討論。分爲兩個情況,每個情況又有三種子情況。可結合下圖分析討論。

這裏先説明下$Chessboard$的度量函數$f(i,y) = \max{(|x-i|, g(i))}$,這個函數限制了最小值為$g(i)$,因此在圖像中會有一段平坦的區域,并且該平坦區域以$u$為中點對稱。

基於這些特點,我們將情況劃分為以下兩大類:

  1. $g(i) \leq g(u)$;
  2. $g(i) \geq g(u)$;

而每個大類中又分爲三個子情況:

a. $F_i(x)$的增區域與$F_u(x)$的減區域相交;

b. $F_i(x)$的增區間與$F_u(x)$的平坦區相交,且與$F_u(x)$減區域的延長綫有交點;

c. $F_i(x)$的增區間與$F_u(x)$的平坦區間相交,且$F_i(x)$增區間延長綫與$F_u(x)$減區間延長綫有交點。

這裏描述的情況可以對應上圖。

現假設$\gamma$表示$i$與$u$中點的函數值,例如上圖圖(a)中 $$ F_i(\gamma) = F_u(\gamma) = F_u \left(\frac{i+u}{2}\right) = \frac{u-i}{2} $$

結合圖(a)~圖(c),交點位置左側均滿足$F_i(x) \leq F_u(x)$,圖(a)交點位置為$x^* = \frac{i+u}{2}$,圖(b)~圖(c)交點位置為$x^* = i + g(u)$。這樣結合起來,在$g(i) \leq g(u)$情況下,交點位置計算如下: $$ g(i) \leq g(u) \Rightarrow (F_i(x) \leq F_u(x)) \Leftrightarrow x \leq \max{(\frac{i+u}{2}, i+g(u))} $$

同理可得圖(d)~圖(f)中的交點位置為: $$ g(i) > g(u) \Rightarrow (F_i(x) > F_u(x)) \Leftrightarrow x \geq \min{(\frac{i+u}{2}, u-g(i))} $$

最終結合可得到$Chessboard$度量的$Sep$函數如下:

$$ Sep(i,u) = \begin{cases} \max{({(i+u)}\ \mathbb{div}\ {2}, i+g(u))},& \text{if }g(i) \leq g(u)\\ \min{({(i+u)}\ \mathbb{div}\ {2}, u-g(i))},& \text{otherwise} \end{cases} $$

這裏另外説明一個問題,爲什麽在CDT中,不會出現$F_i(x)$函數完全包含$F_u(x)$的情況,因爲如果出現這種情況,那麽不能滿足$i <u$這個大前提了,所以只可能出現上面討論的這些情況。

Parallelization, Timing Results, and Conclusions

文中作者也提到該算法可以并行計算,同樣給出了并行計算的時間對比表格,如下;

作者給出了該算法的時間複雜度為$\mathcal{O}(mn)$,如果并行其時間複雜度為$\mathcal{O}(mn/p)$,$p$為并行數。

當然該算法也能很容易推廣到$d$維空間的距離變換下。

至此,該論文閲讀完畢,實現細節請自行斟酌,并不複雜。