基于GVF力的拓?fù)浔３址指钅Ｐ图捌鋵ε妓惴?/h1>
2022-03-16 07:21:26沈夢潔潘振寬宋金濤魏偉波
沈夢潔 潘振寬 宋金濤 魏偉波
摘要：針對自排斥Snake模型對于狹窄圖像區(qū)域作用力不足，傳統(tǒng)加性算子分裂方法計(jì)算復(fù)雜，內(nèi)存用量也會(huì)隨著圖像大小的增加而迅速增長等問題，提出了在原模型的基礎(chǔ)上增加梯度矢量流有向力場，以加快輪廓線在圖像狹窄區(qū)域的演化速度，并為改進(jìn)的模型設(shè)計(jì)快速對偶算法以簡化算法設(shè)計(jì)，提高求解效率。數(shù)值實(shí)驗(yàn)表明，改進(jìn)模型及算法在計(jì)算效率方面較經(jīng)典模型及算法有較大提高。
關(guān)鍵詞：自排斥Snake模型;拓?fù)浔３址指?對偶算法;梯度矢量流;變分法
中圖分類號：TP391 文獻(xiàn)標(biāo)志碼：A
圖像分割是把圖像分成各具特性的區(qū)域并提取感興趣目標(biāo)的技術(shù)和過程，在醫(yī)學(xué)影像[1]、遙感影像[2]、深度學(xué)習(xí)[3]等多個(gè)領(lǐng)域得到了廣泛應(yīng)用。在某些特定情況下，需要在圖像分割的基礎(chǔ)上保持分割目標(biāo)的拓?fù)浣Y(jié)構(gòu)來實(shí)現(xiàn)更好的分割效果，比如在醫(yī)學(xué)圖像中，大腦皮層表面的分割結(jié)果必須與真實(shí)世界大腦皮層結(jié)構(gòu)一致[4]，再比如具有噪聲或遮擋的對象的分割，可以使用拓?fù)浼s束來防止過度分割或欠分割等。學(xué)者在這一領(lǐng)域提出了許多模型及方法，如基于拓?fù)湎闰?yàn)的最小切割/最大流算法[5];具有可變視場拓?fù)浼s束的自動(dòng)磁共振脊髓分割方法[6];能在多標(biāo)簽圖像分割中保持拓?fù)涞姆椒╗7];約束輪廓演化的拓?fù)浔３址椒╗8];非局部的拓?fù)浔３址指钅Ｐ秃突诜蔷植啃螤蠲枋龇姆指钅Ｐ蚚9-10];使用形狀的Beltrami表示進(jìn)行拓?fù)浔３謭D像分割的方法[11];基于超彈性正則化的拓?fù)浔３秩S圖像分割模型[12]等。為使測地線活動(dòng)輪廓模型[13]具有拓?fù)浣Y(jié)構(gòu)的適應(yīng)性，模型求解時(shí)往往采用變分水平集方法[14]。變分水平集方法由于具有自適應(yīng)復(fù)雜拓?fù)浣Y(jié)構(gòu)變化、二維和三維圖像分割模型表達(dá)一致、數(shù)值計(jì)算方法穩(wěn)定以及具有集成多模型信息能力等特點(diǎn)，在圖像分割領(lǐng)域被廣泛應(yīng)用[15]。變分水平集方法可以方便地追蹤物體的拓?fù)浣Y(jié)構(gòu)改變，但對于有粘連的物體分割效果不理想。為了在輪廓演化過程中保持拓?fù)浣Y(jié)構(gòu)，通常在變分公式中添加一個(gè)約束項(xiàng)，以防止輪廓自相交，即合并或分裂，如隱式水平集框架中的簡單點(diǎn)檢測方案，該方案通過引入算法，在每次迭代時(shí)都監(jiān)視水平集函數(shù)的符號變化并防止水平集函數(shù)改變符號不簡單的網(wǎng)格點(diǎn)[16];或?qū)⑼負(fù)浔３謫栴}重鑄為水平集的形狀優(yōu)化問題，使用了符號距離函數(shù)，避免了演化輪廓的窄帶重疊[17];利用基于結(jié)能量最小化的方法來實(shí)現(xiàn)同樣的效果[18];或在引入非局部正則化項(xiàng)時(shí)使用了類似的思想，并應(yīng)用于遙感圖像中細(xì)長物體的跟蹤[19];而后出現(xiàn)了自排斥Snake模型（Self-repelling Snake，SRS）[20]，該模型使用隱式水平集方法表示曲線，并在經(jīng)典測地活動(dòng)輪廓線模型中添加了非局部排斥項(xiàng)[13]。SRS模型對于拓?fù)浔３址指钚Ч^好，但也存在一些問題，如對于狹窄圖像區(qū)域作用力不足，梯度下降方程的推導(dǎo)復(fù)雜，需要使用迎風(fēng)差分離散格式，內(nèi)存用量也會(huì)隨著圖像大小的增加而迅速增長。文獻(xiàn)[21]使用了Split Bregman算法求解原SRS模型，相對于SRS模型，減少了內(nèi)存用量，縮短了迭代時(shí)間。本文針對狹窄圖像區(qū)域作用力不足，設(shè)計(jì)一種加入梯度矢量流有向力場[22]的模型，通過引入對偶變量，利用凸優(yōu)化模型的快速對偶算法[23]，使用投影方法進(jìn)行約束，并對其迭代求解，與原模型及方法在多幅圖像上開展對比實(shí)驗(yàn)，通過對比分割效果、迭代次數(shù)以及運(yùn)行時(shí)間，驗(yàn)證了所提模型及方法的有效性和高效性。
1 相關(guān)研究工作基礎(chǔ)
1.1 SRS模型
SRS模型是在測地活動(dòng)輪廓線模型[10]基礎(chǔ)上加入氣球力項(xiàng)和排斥力項(xiàng)，使用隱式水平集方法來表示曲面，其中分割輪廓隱式表示為符號距離函數(shù)[24]的零水平線，模型的定義為：令fx：Ω→R為標(biāo)量值圖像，x∈Ω，Ω為圖像域，Ω為其邊界，標(biāo)準(zhǔn)邊緣檢測函數(shù)gx
gx=11+ρGσfs（1）
其中，s=1或2，ρ是縮放參數(shù)，Gσ表示圖像的高斯卷積，標(biāo)準(zhǔn)差為σ，對象邊界C用水平集函數(shù)的零水平集表示
C=x∈Ωx=0（2）
水平集函數(shù)被定義符號距離函數(shù)
x=-dx，C? x inside C0????? x∈Cd（x，C）?? x outside C （3）
其中，dx，C是點(diǎn)x與輪廓C之間的歐氏距離。作為一個(gè)符號距離函數(shù)，滿足以下約束條件
=1（4）
為了表示圖像區(qū)域和輪廓，引入了Heaviside函數(shù)H和Dirac函數(shù)δ，由于原始 Heaviside 函數(shù)是不連續(xù)的，因此不可微，所以采用以下正則化方案
Hε=121+ε+1πsinπε ≤ε1??????????? >ε0??????????? <-ε （5）
δε=12ε1+cosπε??? ≤ε0??????????? >ε （6）
SRS模型的能量泛函為
E（）=γ∫ΩgxHεxdx+α∫Ωgx1-Hεxdx-
β∫Ω∫Ωe-x-y2d2x·yhεxhεydxdy（7）
其中，等式右邊第一項(xiàng)為等高線的測地線長度[13]，改寫為γ∫Ωgxxδεxdx;第二項(xiàng)為氣球力項(xiàng)，能夠推動(dòng)分割輪廓越過弱邊緣;第三項(xiàng)為排斥力項(xiàng)，在分割輪廓線較近時(shí)能夠限制拓?fù)浣Y(jié)構(gòu)的改變。γ，α，β為平衡三個(gè)條件的懲罰參數(shù)，e-x-y2d2是用來測量兩個(gè)點(diǎn)x和y的接近度，距離越遠(yuǎn)的點(diǎn)排斥力越小，hεx和hεy表示x點(diǎn)和y點(diǎn)周圍的窄帶（l為實(shí)數(shù)）
hεx=Hεx+l1-Hεx-l（8）
hεy=Hεy+l1-Hεy-l（9）
SRS模型的能量極值問題
minE=γ∫Ωgxxδεxdx+α∫Ωgx1-Hεxdx-????? β∫Ω∫Ωe-x-y2d2x·yhεxhεydxdy s.t.=1 （10）
應(yīng)用變分方法解以上三個(gè)能量項(xiàng)，演化方程為
t=δεxγ·gx+αgx+4βd2hεx∫Ωe-x-y2d2x-y·yhεydy x∈Ωx，0=0x???? t=0=0???????? x∈Ωs.t.=1 （11）
關(guān)于=1的約束，采用以下動(dòng)態(tài)初始化方案
t+signx-1=0x，0=x （12）
式（12）是一個(gè)典型的Hamilton-Jacobi方程，可以通過迎風(fēng)差分格式[24]進(jìn)行離散求解。
1.2 SRS模型的AOS（Additive Operator Splitting）方法
文獻(xiàn)[20]中應(yīng)用的是傳統(tǒng)AOS方法。用半點(diǎn)差分格式和諧波平均近似對式（11）中的第一個(gè)式子等式右邊的第一項(xiàng)進(jìn)行離散，后兩項(xiàng)采用迎風(fēng)方案。通過將圖像的行和列分別連接起來，構(gòu)造了兩個(gè)半隱式格式
1-2τAx1kvk+1=k+τT2k+T3k（13）
1-2τAx2kwk+1=k+τT2k+T3k（14）
其中，Ax1和Ax2為兩個(gè)串聯(lián)矩陣，v和w為中間變量，T2和T3是演化方程中的第一個(gè)式子等式右邊的第二項(xiàng)和第三項(xiàng)的迎風(fēng)離散化，對于每個(gè)All∈x1，x2
Alijk=2γοkiοkigi+οkjgj???? j∈Nli-Σm∈Nl（i）2γοkiοkigi+οkmgm j=i0????????????? else （15）
其中，i，j為圖像中的兩點(diǎn)，Nli是矩陣Al中i的鄰域
οkj=i+1，j-i-1，j22+i，j+1-i，j-122（16）
最后，k+1計(jì)算為
k+1=12vk+1+wk+1（17）
總體上說，AOS方法解決了SRS模型的計(jì)算問題，但是引入了新的參數(shù)，且算法計(jì)算復(fù)雜，內(nèi)存用量也會(huì)隨著圖像大小的增加而迅速增長。
2 ISRS（Improved Self-repelling Snake）模型及對偶算法
2.1 對偶算法
通過使用對偶算法[23]，可以簡化算法設(shè)計(jì)，提高求解效率。對偶算法的主要解決思路是引入對偶變量
γ∫ΩgxHεxdx=max：≤γg（x）∫ΩHεx·dx（18）
對原變量和對偶變量進(jìn)行快速迭代計(jì)算，得到SRS模型的交替優(yōu)化形式為（k為迭代次數(shù)）
k+1，k+1=argmin：=1max：≤γgx∫ΩHεx·dx+α∫Ωgx1-Hεxdx-
β∫Ω∫Ωe-x-y2d2x·yhεxhεydxdy? （19）
交替優(yōu)化方法通過固定其他變量，從而求解一個(gè)變量的極值問題，式（19）可轉(zhuǎn)化為以下子優(yōu)化問題
k+1=argmax E1=Ek，s.t.≤γgx （20）
k+1=argmin E2=E，k+1（21）
首先，固定變量k，由式（19）得到關(guān)于k+1的能量泛函極值問題
k+1=max≤γgxEk，=max≤γgx∫ΩHεkx·xdx （22）
求解式（22），得到的梯度下降方程
t=-Hεkxs.t.≤γgx （23）
得
x=kx+τkpHεkx（24）
求對偶變量的解，得
k+1x=kx+τkpHεkxs.t. k+1≤γgx （25）
由式（25）得的投影公式
k+1xγgxk+1xmaxk+1（x），γgx（26）
固定k+1，求解k+1，由式（19）得到關(guān)于的能量泛函極值問題
k+1=min：=1E，k+1=min：=1∫ΩHεx·k+1dx+α∫Ωgx1-Hεxdx-
β∫Ω∫Ωe-x-y2d2x·yhεxhεydxdy （27）
對式（27）求偏導(dǎo)，得
E+εηε=∫Ω·k+1δεxηxdx-α∫Ωgxδεxηxdx-
4βd2∫Ωhεx∫Ωe-x-y2d2x-y·yhεydyηxdx（28）
其中，ε為一個(gè)極小增量，η為一個(gè)與性質(zhì)相同的任意函數(shù)。
通過梯度降方法來求解，得
t=αgx-·k+1δεkx+?? 4βd2hεkx∫Ωe-x-y2d2x-y·kyhεkydy=1 （29）
整理，得
t=αgx-·k+1δεkx+?? 4βd2hεkx∫Ωe-x-y2d2x-y·kyhεkydyt+signx-1=0 （30）
2.2 GVF模型
本文提出的方法添加了GVF（Gradient Vector Flow）模型[22]，該模型可以有效地克服模型不能收斂到輪廓凹陷的問題，且該模型是基于外力矢量的能量函數(shù)構(gòu)造的，利用該能量函數(shù)可以產(chǎn)生約束輪廓向邊界演化的力。其中GVF力的力場為：x，y=ux，y，vx，y，能量泛函定義為
E=α∫ 2dx+∫f 2-f2dx（31）
其中，f為與圖像數(shù)據(jù)有關(guān)的外部能量項(xiàng)，f為外部作用力。α為懲罰參數(shù)用來平衡能量項(xiàng)中的第一項(xiàng)與第二項(xiàng)，α越大，該力場就越平滑。在圖像邊緣附近f2較大，外部能量項(xiàng)（第二項(xiàng)）起主要作用，此時(shí)趨近于f。而在灰度比較均勻的區(qū)域，式中第二項(xiàng)為零，只有第一項(xiàng)起作用，此時(shí)具有平滑的作用。
令Qk=αgx-·k+1，用代替δεx，并將GVF力·加入到式（30）當(dāng)中，得
t=Qk+·+4βd2hεkx∫Ωe-x-y2d2x-y·kyhεkydy（32）
2.3 離散化及迭代方案
對初始化
k + 1，0i，j = ki，j（33）
令
Qki，j=αgi，j--·k+1i，j（34）
使用迎風(fēng)差分格式離散求解，得
k+1，li，j+=max-x1k+1，li，j，02+min+x1k+1，li，j，02+max-x2k+1，li，j，02+min+x2k+1，li，j，02（35）
k+1，li，j-=min-x1k+1，li，j，02+max+x1k+1，li，j，02+min-x2k+1，li，j，02+max+x2k+1，li，j，02（36）
其中
+i，j=i+1，j-i，j（37）
-i，j=i，j-i-1，j（38）
式（32）的求解結(jié)果為
k+1，l+1i，j-k+1，li，jτ=maxQki，j，0k+1，li，j-+minQki，j，0k+1，li，j++maxF1，0+x1k+1，li，j+
minF1，0-x1k+1，li，j+maxF2，0+x2k+1，li，j+minF2，0-x2k+1，li，j+
4βd2hεk+1，li，jx∫Ωe-x-y2d2x-y·k+1，li，jyhεk+1，li，jydy（39）
整理，得
k+1，l+1i，j=k+1，li，j+τ＼[maxQki，j，0k+1，li，j-+minQki，j，0k+1，li，j++maxF1，0+x1k+1，li，j+
minF1，0-x1k+1，li，j+maxF2，0+x2k+1，li，j+minF2，0-x2k+1，li，j+
4βd2hεk+1，li，jx∫Ωe-x-y2d2x-y·k+1，li，jyhεk+1，li，jydy＼]（40）
ISRS模型的對偶算法為：
1）初始化。根據(jù)式（1）計(jì)算gx，初始化0x為符號距離函數(shù)并設(shè)置0=0;
2）設(shè)置懲罰參數(shù)的值;
3）設(shè)置時(shí)間步長和迭代次數(shù);
4）迭代過程：
for i =0，1，2，…，step
compute k+1 from （26）;
for k=0，1，2，…，K
for s=0，1，2，…，S
calculate k+1，s+1 from （40）
end for s
end for k
end for i
3 數(shù)值實(shí)驗(yàn)與分析
本部分針對本文提出的ISRS模型的對偶算法進(jìn)行分割實(shí)驗(yàn)，并與經(jīng)典SRS模型的AOS方法進(jìn)行比較，對比兩者的拓?fù)浔３值姆指钚Ч坝?jì)算效率。為了客觀比較實(shí)驗(yàn)結(jié)果，對每個(gè)實(shí)驗(yàn)中的兩種方法使用相同的初始化位置。實(shí)驗(yàn)中用黑色線條表示分割區(qū)域邊界。實(shí)驗(yàn)環(huán)境：CPU為Intel（R）Core（TM）i5-8500T，內(nèi)存8G，操作系統(tǒng)為Windows10，軟件平臺(tái)為Matlab2018b。
圖1為SRS模型的AOS方法分割合成手部圖像的效果圖，圖2為ISRS模型的對偶算法分割合成手部圖像的效果圖。
圖1 SRS模型的AOS方法分割合成手部圖像
（a）初始輪廓線;（b）迭代200次效果;（c）迭代500次效果;（d）迭代900次效果
圖2 ISRS模型的對偶算法分割合成手部圖像
（a）初始輪廓線;（b）迭代20次效果;（c）迭代40次效果;（d）迭代60次效果;（e）迭代100次效果;（f）迭代200次效果
圖1參數(shù)分別為1、0.5、1、2×2、-0.2、-0.2、1，其中h、τ、l均根據(jù)文獻(xiàn)[20]中參數(shù)設(shè)置;window表示搜索框的大小，α，β，ε均根據(jù)文獻(xiàn)[20]中參數(shù)做了適量調(diào)整，圖像大小127×150像素。圖2參數(shù)γ、α、β、ε、τ、τ2、τ3、l、d、window分別為1、10、0.01、0.1、0.5、0.2、0.125、1、4、3×3，其中，γ為Snake模型的參數(shù)，用來保持正常的分割作用;α為氣球力項(xiàng)參數(shù)，能使模型在平滑區(qū)域加速分割;β為排斥力項(xiàng)參數(shù)，用于保持0水平集的拓?fù)浣Y(jié)構(gòu);ε為域值，讓水平線函數(shù)保持0-1的形態(tài)，從而使函數(shù)變?yōu)榉侄纬Ｖ岛瘮?shù)，便于計(jì)算和保持良好的形態(tài);τ為以及k迭代過程中的時(shí)間步長;τ2為s迭代過程中的時(shí)間步長;τ3為迭代過程中的時(shí)間步長;以上三項(xiàng)均根據(jù)時(shí)間步長基本原則進(jìn)行取值，為長期實(shí)驗(yàn)所取得的經(jīng)驗(yàn)值;l為拓?fù)浔３志嚯x，l越大，拓?fù)浔３志嚯x越遠(yuǎn);d為高斯卷積核的標(biāo)準(zhǔn)差， window為搜索框大小;圖像大小300×354像素。
對比圖1圖2兩組圖像，可以看出兩組運(yùn)行結(jié)果基本一致，但本文改進(jìn)的模型及方法相比原來的模型及AOS方法迭代次數(shù)有了大幅度地減少。本文改進(jìn)的模型及方法在迭代到第100步時(shí)便達(dá)到收斂條件，完成分割，并且拓?fù)浣Y(jié)構(gòu)得到了保持，而傳統(tǒng)的模型及方法900步才完成分割。
圖3為SRS模型的AOS方法分割人為噪聲圖像的效果圖，圖4為ISRS模型的對偶算法分割人為噪聲圖像的效果圖。圖3參數(shù)分別為1、0.5、1、2×2、-0.2、-0.2、1，圖像大小200×138像素。圖4參數(shù)γ、α、β、ε、τ、τ2、τ3、l、d、window分別為8、10、0.01、0.1、0.5、0.2、0.125、1、4、3×3，圖像大小400×275像素。
圖3 SRS模型的AOS方法分割人為噪聲圖像
（a）初始輪廓線;（b）迭代200次效果;（c）迭代400次效果;（d）迭代800次效果
圖4 ISRS模型的對偶算法分割人為噪聲圖像
（a）初始輪廓線;（b）迭代10次效果;（c）迭代20次效果;（d）迭代40次效果;（e）迭代80次效果;（f）迭代200次效果
對比圖3圖4，兩種方法的運(yùn)行結(jié)果都較好地保持了圖像本身的拓?fù)浣Y(jié)構(gòu)，但本文改進(jìn)的模型及方法在迭代到第80步時(shí)便達(dá)到收斂條件，完成分割，并且拓?fù)浣Y(jié)構(gòu)得到了保持，而傳統(tǒng)模型及AOS方法在第800步才完成分割。
圖5為SRS模型的AOS方法分割腦部切片1的效果圖，圖6為ISRS模型的對偶算法分割腦部切片1的效果圖。圖5參數(shù)h、τ、l、window、α、β、ε分別為1、0.5、1、2×2、-0.2、-0.2、1，圖像大小125×125像素。圖6參數(shù)γ、α、β、ε、τ、τ2、τ3、l、d、window分別為1、50、0.01、1、0.3、0.12、0.125、1、4、5×5，圖像大小300×300像素。
圖5 SRS模型的AOS方法分割腦部切片1
（a）初始輪廓線;（b）迭代200次效果;（c）迭代400次效果;（d）迭代700次效果
圖6 ISRS模型的對偶算法分割腦部切片1
（a）初始輪廓線;（b）迭代10次效果;（c）迭代20次效果;（d）迭代30次效果;（e）迭代50次效果;（f）迭代100次效果
對比圖5圖6，本文改進(jìn)的模型及方法同原模型的AOS方法相比能更好地進(jìn)入到凹槽內(nèi)，并且在迭代次數(shù)上也有了大幅度地減少。本文改進(jìn)的模型及方法在迭代到第50步時(shí)便達(dá)到收斂條件，完成分割，并且拓?fù)浣Y(jié)構(gòu)得到了保持，而傳統(tǒng)模型及AOS方法在第700步才完成分割。
圖7為SRS模型的AOS方法分割腦部切片2的效果圖，圖8為ISRS模型的對偶算法分割腦部切片2的效果圖。圖7參數(shù)h、τ、l、window、α、β、ε分別為1、0.5、1、2×2、-0.2、-0.2、1，圖像大小125×125像素。圖8參數(shù)γ、α、β、ε、τ、τ2、τ3、l、d、window分別為1、50、0.01、10、0.3、0.12、0.125、1、4、5×5，圖像大小300×300像素。
（a）初始輪廓線;（b）迭代200次效果;（c）迭代400次效果;（d）迭代800次效果
（a）初始輪廓線;（b）迭代10次效果;（c）迭代20次效果;（d）迭代30次效果;（e）迭代50次效果;（f）迭代100次效果
對比圖7圖8，兩種方法都保持了對象本身的拓?fù)浣Y(jié)構(gòu)，本文改進(jìn)的模型及方法在迭代到第50步時(shí)便達(dá)到收斂條件，完成分割，并且拓?fù)浣Y(jié)構(gòu)得到了保持，而傳統(tǒng)模型及AOS方法在第800步才完成分割。
表1給出了以上四個(gè)對比實(shí)驗(yàn)的迭代次數(shù)、運(yùn)算時(shí)間以及運(yùn)算時(shí)間加速比的具體數(shù)值?？芍谟袇?shù)影響的情況下，本文也能較好地保持分割對象的拓?fù)浣Y(jié)構(gòu)，并在保持對象拓?fù)浣Y(jié)構(gòu)的基礎(chǔ)上大幅度減少了算法的迭代次數(shù)，模型的收斂速度也得到了大幅度提高。
4 結(jié)論
為了進(jìn)一步提高拓?fù)浔３址指钅Ｐ偷挠?jì)算效率，本文提出了ISRS模型及其對偶算法，簡化了計(jì)算過程，提高了運(yùn)算速度;GVF有向力場明顯加快了輪廓線在圖像狹窄區(qū)域的演化速度。對比實(shí)驗(yàn)表明，本文模型及其對偶算法的解類似于經(jīng)典SRS模型及其方法的解，可以有效防止分割輪廓的合并或分裂，保持原來的拓?fù)浣Y(jié)構(gòu)，并且在計(jì)算效率方面較經(jīng)典模型與算法有較大提高。下一步工作可以嘗試用高階模型來忽略分割過程中的微小細(xì)節(jié)，保持邊緣的紋理特征。
參考文獻(xiàn)
[1]國凱，王國棟，潘振寬，等.CT影像中的主動(dòng)脈血管鈣化檢測[J].青島大學(xué)學(xué)報(bào)（自然科學(xué)版），2013，26（1）：50-54.
[2]田佑仕，黃寶香，姜濤，等.融合圖像恢復(fù)的遙感數(shù)據(jù)變分分割方法研究[J].青島大學(xué)學(xué)報(bào)（自然科學(xué)版），2020，33（4）：1-8.
[3]吳則舉，王嘉琦，焦翠娟，等.基于深度學(xué)習(xí)網(wǎng)絡(luò)U-Net的輪胎帶束層分割算法研究[J].青島大學(xué)學(xué)報(bào)（自然科學(xué)版），2019，32（4）：22-29+35.
[4]MACDONALD D， KABANI N， AVIS D， et al. Automated 3-D extraction of inner and outer surfaces of cerebral cortex from MRI[J]. NeuroImage， 2000， 12 （3）：340-356.
[5]ZENG Y， SAMARAS D， CHEN W， et al. Topology cuts： A novel min-cut/max-flow algorithm for topology preserving segmentation in N-D images[J]. Computer Vision and Image Understanding， 2008， 112 （1）：81-90.
[6]CHEN M， CARASS A， OH J， et al. Automatic magnetic resonance spinal cord segmentation with topology constraints for variable fields of view[J]. NeuroImage， 2013，83： 1051-1062.
[7]WAGGONER J， ZHOU Y J， SIMMONS J， et al. Topology-preserving multi-label image segmentation[C]// IEEE Winter Conference on Applications of Computer Vision （WACV 2015）， Waikoloa， 2015： 1084-1091.
[8]LIU G Q， LI H F， YANG L. A topology preserving method of evolving contours based on sparsity constraint for object segmentation[J]. IEEE Access， 2017， 5：19971-19982.
[9]DEBROUX N， OZERE， S， LE GUYADER C. A non-local topology-preserving segmentation-guided registration model[J]. Journal of Mathematical Imaging and Vision， 2017， 59 （3）：432-455.
[10] DEBROUX N， LE GUYADER C. A joint segmentation/registration model based on a nonlocal characterization of weighted total variation and nonlocal shape descriptors[J]. SIAM Journal on Imaging Sciences， 2018， 11 （2）：957-990.
[11] CHAN H L， YAN S， LUI L M， et al. Topology-preserving image segmentation by beltrami representation of shapes[J]. Journal of Mathematical Imaging and Vision， 2018， 60 （3）：401-421.
[12] ZHANG D P， LUI L M. Topology-Preserving 3D image segmentation based on hyperelastic regularization[J]. Journal of Scientific Computing， 2021， 87（3）： 74.
[13] CASELLES V， KIMMEL R， SAPIRO G. Geodesic active contour. IJCV[J]. International Journal of Computer Vision， 1997， 22（1）：61-79.
[14] ZHAO H， CHAN T， MERRIMAN B， et al. A variational level set approach to multiphase motion[J]. Journal of Computational Physics， 1996， 127（1）： 179-195.
[15] 郭振波.基于變分水平集方法的多相圖像分割研究[D]. 青島： 中國海洋大學(xué)， 2008.
[16] HAN X， XU C Y， PRINCE J L. A topology preserving level set method for geometric deformable models[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence， 2003， 25（6）：755-768.
[17] ALEXANDROV O， SANTOSA F. A topology-preserving level set method for shape optimization[J]. Journal of Computational Physics， 2005， 204（1）：121-130.
[18] SUNDARAMOORTHI G， YEZZI A. Global regularizing flows with topology preservation for active contours and polygons[J]. IEEE Transactions on Image Processing， 2007， 16（3）：803-812.
[19] ROCHERY M， JERMYN I H， ZERUBIA J. Higher order active contours[J]. International Journal of Computer Vision， 2006， 69（1）：27-42.
[20] GUYADER C L， VESE L A. Self-repelling snakes for topology-preserving segmentation models[J]. IEEE Transactions on Image Processing， 2008， 17（5）：767-779.
[21] PAN H， SONG J， LIU W， et al. Using the split bregman algorithm to solve the self-repelling snake model＼[J＼]. Computer Vision and Pattern Recognition， 2021＼[2021-08-15＼]https：//doi.org/10.1007/s10851-021-01065-9.
[22] XU C， PRINCE J L. Snakes， shapes， and gradient vector flow[J]. IEEE Transactions on Image Processing， 1998：7（3）：359-369.
[23] CHAMBOLLE A. An algorithm for total variation minimization and applications[J]. Journal of Mathematical Imaging and Vision， 2004， 20（1-2）： 89-97.
[24] OSHER S， SETHIAN J A. Fronts propagating with curvature-dependent speed： Algorithms based on Hamilton-Jacobi formulations＼[J＼]. Journal of Computational Physics， 1988， 79 （1）： 12-49.
Topology Preserving Segmentation Model Based on
GVF Force and Its Dual Algorithm
SHEN Meng-jie， PAN Zhen-kuan， SONG Jin-tao， WEI Wei-bo
（College of Computer Science & Technology， Qingdao University， Qingdao 266071， China）
Abstract：
In view of the insufficient force of the self-repelling Snake model on narrow image regions， the calculation of traditional additive operator splitting method is complicated， and the memory usage will also increase rapidly with the increase of the image size. It is proposed that the Gradient Vector Flow directed force field based on the original model was added to accelerate the evolution speed of the contour line in the narrow area of the image， and a fast Dual algorithm is designed for the improved model to simplify the algorithm design and improve the efficiency of the solution. Numerous numerical experiments show that the proposed improved model and algorithm have a greater improvement in computational efficiency than the classic model and algorithm.
Keywords：
self-repelling snake model; topology-preserving segmentation; dual algorithm; gradient vector flow; variational method
收稿日期：2021-09-10
基金項(xiàng)目：
國家自然科學(xué)基金（批準(zhǔn)號：61772294，11472144） 資助;山東省聯(lián)合基金（批準(zhǔn)號：ZR2019LZH002） 資助。
通信作者：
潘振寬，男，博士，教授，主要研究方向?yàn)橛?jì)算機(jī)視覺、圖像處理，E-mail： zkpan@126.com
沈夢潔 潘振寬 宋金濤 魏偉波
摘要：針對自排斥Snake模型對于狹窄圖像區(qū)域作用力不足，傳統(tǒng)加性算子分裂方法計(jì)算復(fù)雜，內(nèi)存用量也會(huì)隨著圖像大小的增加而迅速增長等問題，提出了在原模型的基礎(chǔ)上增加梯度矢量流有向力場，以加快輪廓線在圖像狹窄區(qū)域的演化速度，并為改進(jìn)的模型設(shè)計(jì)快速對偶算法以簡化算法設(shè)計(jì)，提高求解效率。數(shù)值實(shí)驗(yàn)表明，改進(jìn)模型及算法在計(jì)算效率方面較經(jīng)典模型及算法有較大提高。
關(guān)鍵詞：自排斥Snake模型;拓?fù)浔３址指?對偶算法;梯度矢量流;變分法
中圖分類號：TP391 文獻(xiàn)標(biāo)志碼：A
圖像分割是把圖像分成各具特性的區(qū)域并提取感興趣目標(biāo)的技術(shù)和過程，在醫(yī)學(xué)影像[1]、遙感影像[2]、深度學(xué)習(xí)[3]等多個(gè)領(lǐng)域得到了廣泛應(yīng)用。在某些特定情況下，需要在圖像分割的基礎(chǔ)上保持分割目標(biāo)的拓?fù)浣Y(jié)構(gòu)來實(shí)現(xiàn)更好的分割效果，比如在醫(yī)學(xué)圖像中，大腦皮層表面的分割結(jié)果必須與真實(shí)世界大腦皮層結(jié)構(gòu)一致[4]，再比如具有噪聲或遮擋的對象的分割，可以使用拓?fù)浼s束來防止過度分割或欠分割等。學(xué)者在這一領(lǐng)域提出了許多模型及方法，如基于拓?fù)湎闰?yàn)的最小切割/最大流算法[5];具有可變視場拓?fù)浼s束的自動(dòng)磁共振脊髓分割方法[6];能在多標(biāo)簽圖像分割中保持拓?fù)涞姆椒╗7];約束輪廓演化的拓?fù)浔３址椒╗8];非局部的拓?fù)浔３址指钅Ｐ秃突诜蔷植啃螤蠲枋龇姆指钅Ｐ蚚9-10];使用形狀的Beltrami表示進(jìn)行拓?fù)浔３謭D像分割的方法[11];基于超彈性正則化的拓?fù)浔３秩S圖像分割模型[12]等。為使測地線活動(dòng)輪廓模型[13]具有拓?fù)浣Y(jié)構(gòu)的適應(yīng)性，模型求解時(shí)往往采用變分水平集方法[14]。變分水平集方法由于具有自適應(yīng)復(fù)雜拓?fù)浣Y(jié)構(gòu)變化、二維和三維圖像分割模型表達(dá)一致、數(shù)值計(jì)算方法穩(wěn)定以及具有集成多模型信息能力等特點(diǎn)，在圖像分割領(lǐng)域被廣泛應(yīng)用[15]。變分水平集方法可以方便地追蹤物體的拓?fù)浣Y(jié)構(gòu)改變，但對于有粘連的物體分割效果不理想。為了在輪廓演化過程中保持拓?fù)浣Y(jié)構(gòu)，通常在變分公式中添加一個(gè)約束項(xiàng)，以防止輪廓自相交，即合并或分裂，如隱式水平集框架中的簡單點(diǎn)檢測方案，該方案通過引入算法，在每次迭代時(shí)都監(jiān)視水平集函數(shù)的符號變化并防止水平集函數(shù)改變符號不簡單的網(wǎng)格點(diǎn)[16];或?qū)⑼負(fù)浔３謫栴}重鑄為水平集的形狀優(yōu)化問題，使用了符號距離函數(shù)，避免了演化輪廓的窄帶重疊[17];利用基于結(jié)能量最小化的方法來實(shí)現(xiàn)同樣的效果[18];或在引入非局部正則化項(xiàng)時(shí)使用了類似的思想，并應(yīng)用于遙感圖像中細(xì)長物體的跟蹤[19];而后出現(xiàn)了自排斥Snake模型（Self-repelling Snake，SRS）[20]，該模型使用隱式水平集方法表示曲線，并在經(jīng)典測地活動(dòng)輪廓線模型中添加了非局部排斥項(xiàng)[13]。SRS模型對于拓?fù)浔３址指钚Ч^好，但也存在一些問題，如對于狹窄圖像區(qū)域作用力不足，梯度下降方程的推導(dǎo)復(fù)雜，需要使用迎風(fēng)差分離散格式，內(nèi)存用量也會(huì)隨著圖像大小的增加而迅速增長。文獻(xiàn)[21]使用了Split Bregman算法求解原SRS模型，相對于SRS模型，減少了內(nèi)存用量，縮短了迭代時(shí)間。本文針對狹窄圖像區(qū)域作用力不足，設(shè)計(jì)一種加入梯度矢量流有向力場[22]的模型，通過引入對偶變量，利用凸優(yōu)化模型的快速對偶算法[23]，使用投影方法進(jìn)行約束，并對其迭代求解，與原模型及方法在多幅圖像上開展對比實(shí)驗(yàn)，通過對比分割效果、迭代次數(shù)以及運(yùn)行時(shí)間，驗(yàn)證了所提模型及方法的有效性和高效性。
1 相關(guān)研究工作基礎(chǔ)
1.1 SRS模型
SRS模型是在測地活動(dòng)輪廓線模型[10]基礎(chǔ)上加入氣球力項(xiàng)和排斥力項(xiàng)，使用隱式水平集方法來表示曲面，其中分割輪廓隱式表示為符號距離函數(shù)[24]的零水平線，模型的定義為：令fx：Ω→R為標(biāo)量值圖像，x∈Ω，Ω為圖像域，Ω為其邊界，標(biāo)準(zhǔn)邊緣檢測函數(shù)gx
gx=11+ρGσfs（1）
其中，s=1或2，ρ是縮放參數(shù)，Gσ表示圖像的高斯卷積，標(biāo)準(zhǔn)差為σ，對象邊界C用水平集函數(shù)的零水平集表示
C=x∈Ωx=0（2）
水平集函數(shù)被定義符號距離函數(shù)
x=-dx，C? x inside C0????? x∈Cd（x，C）?? x outside C （3）
其中，dx，C是點(diǎn)x與輪廓C之間的歐氏距離。作為一個(gè)符號距離函數(shù)，滿足以下約束條件
=1（4）
為了表示圖像區(qū)域和輪廓，引入了Heaviside函數(shù)H和Dirac函數(shù)δ，由于原始 Heaviside 函數(shù)是不連續(xù)的，因此不可微，所以采用以下正則化方案
Hε=121+ε+1πsinπε ≤ε1??????????? >ε0??????????? <-ε （5）
δε=12ε1+cosπε??? ≤ε0??????????? >ε （6）
SRS模型的能量泛函為
E（）=γ∫ΩgxHεxdx+α∫Ωgx1-Hεxdx-
β∫Ω∫Ωe-x-y2d2x·yhεxhεydxdy（7）
其中，等式右邊第一項(xiàng)為等高線的測地線長度[13]，改寫為γ∫Ωgxxδεxdx;第二項(xiàng)為氣球力項(xiàng)，能夠推動(dòng)分割輪廓越過弱邊緣;第三項(xiàng)為排斥力項(xiàng)，在分割輪廓線較近時(shí)能夠限制拓?fù)浣Y(jié)構(gòu)的改變。γ，α，β為平衡三個(gè)條件的懲罰參數(shù)，e-x-y2d2是用來測量兩個(gè)點(diǎn)x和y的接近度，距離越遠(yuǎn)的點(diǎn)排斥力越小，hεx和hεy表示x點(diǎn)和y點(diǎn)周圍的窄帶（l為實(shí)數(shù)）
hεx=Hεx+l1-Hεx-l（8）
hεy=Hεy+l1-Hεy-l（9）
SRS模型的能量極值問題
minE=γ∫Ωgxxδεxdx+α∫Ωgx1-Hεxdx-????? β∫Ω∫Ωe-x-y2d2x·yhεxhεydxdy s.t.=1 （10）
應(yīng)用變分方法解以上三個(gè)能量項(xiàng)，演化方程為
t=δεxγ·gx+αgx+4βd2hεx∫Ωe-x-y2d2x-y·yhεydy x∈Ωx，0=0x???? t=0=0???????? x∈Ωs.t.=1 （11）
關(guān)于=1的約束，采用以下動(dòng)態(tài)初始化方案
t+signx-1=0x，0=x （12）
式（12）是一個(gè)典型的Hamilton-Jacobi方程，可以通過迎風(fēng)差分格式[24]進(jìn)行離散求解。
1.2 SRS模型的AOS（Additive Operator Splitting）方法
文獻(xiàn)[20]中應(yīng)用的是傳統(tǒng)AOS方法。用半點(diǎn)差分格式和諧波平均近似對式（11）中的第一個(gè)式子等式右邊的第一項(xiàng)進(jìn)行離散，后兩項(xiàng)采用迎風(fēng)方案。通過將圖像的行和列分別連接起來，構(gòu)造了兩個(gè)半隱式格式
1-2τAx1kvk+1=k+τT2k+T3k（13）
1-2τAx2kwk+1=k+τT2k+T3k（14）
其中，Ax1和Ax2為兩個(gè)串聯(lián)矩陣，v和w為中間變量，T2和T3是演化方程中的第一個(gè)式子等式右邊的第二項(xiàng)和第三項(xiàng)的迎風(fēng)離散化，對于每個(gè)All∈x1，x2
Alijk=2γοkiοkigi+οkjgj???? j∈Nli-Σm∈Nl（i）2γοkiοkigi+οkmgm j=i0????????????? else （15）
其中，i，j為圖像中的兩點(diǎn)，Nli是矩陣Al中i的鄰域
οkj=i+1，j-i-1，j22+i，j+1-i，j-122（16）
最后，k+1計(jì)算為
k+1=12vk+1+wk+1（17）
總體上說，AOS方法解決了SRS模型的計(jì)算問題，但是引入了新的參數(shù)，且算法計(jì)算復(fù)雜，內(nèi)存用量也會(huì)隨著圖像大小的增加而迅速增長。
2 ISRS（Improved Self-repelling Snake）模型及對偶算法
2.1 對偶算法
通過使用對偶算法[23]，可以簡化算法設(shè)計(jì)，提高求解效率。對偶算法的主要解決思路是引入對偶變量
γ∫ΩgxHεxdx=max：≤γg（x）∫ΩHεx·dx（18）
對原變量和對偶變量進(jìn)行快速迭代計(jì)算，得到SRS模型的交替優(yōu)化形式為（k為迭代次數(shù)）
k+1，k+1=argmin：=1max：≤γgx∫ΩHεx·dx+α∫Ωgx1-Hεxdx-
β∫Ω∫Ωe-x-y2d2x·yhεxhεydxdy? （19）
交替優(yōu)化方法通過固定其他變量，從而求解一個(gè)變量的極值問題，式（19）可轉(zhuǎn)化為以下子優(yōu)化問題
k+1=argmax E1=Ek，s.t.≤γgx （20）
k+1=argmin E2=E，k+1（21）
首先，固定變量k，由式（19）得到關(guān)于k+1的能量泛函極值問題
k+1=max≤γgxEk，=max≤γgx∫ΩHεkx·xdx （22）
求解式（22），得到的梯度下降方程
t=-Hεkxs.t.≤γgx （23）
得
x=kx+τkpHεkx（24）
求對偶變量的解，得
k+1x=kx+τkpHεkxs.t. k+1≤γgx （25）
由式（25）得的投影公式
k+1xγgxk+1xmaxk+1（x），γgx（26）
固定k+1，求解k+1，由式（19）得到關(guān)于的能量泛函極值問題
k+1=min：=1E，k+1=min：=1∫ΩHεx·k+1dx+α∫Ωgx1-Hεxdx-
β∫Ω∫Ωe-x-y2d2x·yhεxhεydxdy （27）
對式（27）求偏導(dǎo)，得
E+εηε=∫Ω·k+1δεxηxdx-α∫Ωgxδεxηxdx-
4βd2∫Ωhεx∫Ωe-x-y2d2x-y·yhεydyηxdx（28）
其中，ε為一個(gè)極小增量，η為一個(gè)與性質(zhì)相同的任意函數(shù)。
通過梯度降方法來求解，得
t=αgx-·k+1δεkx+?? 4βd2hεkx∫Ωe-x-y2d2x-y·kyhεkydy=1 （29）
整理，得
t=αgx-·k+1δεkx+?? 4βd2hεkx∫Ωe-x-y2d2x-y·kyhεkydyt+signx-1=0 （30）
2.2 GVF模型
本文提出的方法添加了GVF（Gradient Vector Flow）模型[22]，該模型可以有效地克服模型不能收斂到輪廓凹陷的問題，且該模型是基于外力矢量的能量函數(shù)構(gòu)造的，利用該能量函數(shù)可以產(chǎn)生約束輪廓向邊界演化的力。其中GVF力的力場為：x，y=ux，y，vx，y，能量泛函定義為
E=α∫ 2dx+∫f 2-f2dx（31）
其中，f為與圖像數(shù)據(jù)有關(guān)的外部能量項(xiàng)，f為外部作用力。α為懲罰參數(shù)用來平衡能量項(xiàng)中的第一項(xiàng)與第二項(xiàng)，α越大，該力場就越平滑。在圖像邊緣附近f2較大，外部能量項(xiàng)（第二項(xiàng)）起主要作用，此時(shí)趨近于f。而在灰度比較均勻的區(qū)域，式中第二項(xiàng)為零，只有第一項(xiàng)起作用，此時(shí)具有平滑的作用。
令Qk=αgx-·k+1，用代替δεx，并將GVF力·加入到式（30）當(dāng)中，得
t=Qk+·+4βd2hεkx∫Ωe-x-y2d2x-y·kyhεkydy（32）
2.3 離散化及迭代方案
對初始化
k + 1，0i，j = ki，j（33）
令
Qki，j=αgi，j--·k+1i，j（34）
使用迎風(fēng)差分格式離散求解，得
k+1，li，j+=max-x1k+1，li，j，02+min+x1k+1，li，j，02+max-x2k+1，li，j，02+min+x2k+1，li，j，02（35）
k+1，li，j-=min-x1k+1，li，j，02+max+x1k+1，li，j，02+min-x2k+1，li，j，02+max+x2k+1，li，j，02（36）
其中
+i，j=i+1，j-i，j（37）
-i，j=i，j-i-1，j（38）
式（32）的求解結(jié)果為
k+1，l+1i，j-k+1，li，jτ=maxQki，j，0k+1，li，j-+minQki，j，0k+1，li，j++maxF1，0+x1k+1，li，j+
minF1，0-x1k+1，li，j+maxF2，0+x2k+1，li，j+minF2，0-x2k+1，li，j+
4βd2hεk+1，li，jx∫Ωe-x-y2d2x-y·k+1，li，jyhεk+1，li，jydy（39）
整理，得
k+1，l+1i，j=k+1，li，j+τ＼[maxQki，j，0k+1，li，j-+minQki，j，0k+1，li，j++maxF1，0+x1k+1，li，j+
minF1，0-x1k+1，li，j+maxF2，0+x2k+1，li，j+minF2，0-x2k+1，li，j+
4βd2hεk+1，li，jx∫Ωe-x-y2d2x-y·k+1，li，jyhεk+1，li，jydy＼]（40）
ISRS模型的對偶算法為：
1）初始化。根據(jù)式（1）計(jì)算gx，初始化0x為符號距離函數(shù)并設(shè)置0=0;
2）設(shè)置懲罰參數(shù)的值;
3）設(shè)置時(shí)間步長和迭代次數(shù);
4）迭代過程：
for i =0，1，2，…，step
compute k+1 from （26）;
for k=0，1，2，…，K
for s=0，1，2，…，S
calculate k+1，s+1 from （40）
end for s
end for k
end for i
3 數(shù)值實(shí)驗(yàn)與分析
本部分針對本文提出的ISRS模型的對偶算法進(jìn)行分割實(shí)驗(yàn)，并與經(jīng)典SRS模型的AOS方法進(jìn)行比較，對比兩者的拓?fù)浔３值姆指钚Ч坝?jì)算效率。為了客觀比較實(shí)驗(yàn)結(jié)果，對每個(gè)實(shí)驗(yàn)中的兩種方法使用相同的初始化位置。實(shí)驗(yàn)中用黑色線條表示分割區(qū)域邊界。實(shí)驗(yàn)環(huán)境：CPU為Intel（R）Core（TM）i5-8500T，內(nèi)存8G，操作系統(tǒng)為Windows10，軟件平臺(tái)為Matlab2018b。
圖1為SRS模型的AOS方法分割合成手部圖像的效果圖，圖2為ISRS模型的對偶算法分割合成手部圖像的效果圖。
圖1 SRS模型的AOS方法分割合成手部圖像
（a）初始輪廓線;（b）迭代200次效果;（c）迭代500次效果;（d）迭代900次效果
圖2 ISRS模型的對偶算法分割合成手部圖像
（a）初始輪廓線;（b）迭代20次效果;（c）迭代40次效果;（d）迭代60次效果;（e）迭代100次效果;（f）迭代200次效果
圖1參數(shù)分別為1、0.5、1、2×2、-0.2、-0.2、1，其中h、τ、l均根據(jù)文獻(xiàn)[20]中參數(shù)設(shè)置;window表示搜索框的大小，α，β，ε均根據(jù)文獻(xiàn)[20]中參數(shù)做了適量調(diào)整，圖像大小127×150像素。圖2參數(shù)γ、α、β、ε、τ、τ2、τ3、l、d、window分別為1、10、0.01、0.1、0.5、0.2、0.125、1、4、3×3，其中，γ為Snake模型的參數(shù)，用來保持正常的分割作用;α為氣球力項(xiàng)參數(shù)，能使模型在平滑區(qū)域加速分割;β為排斥力項(xiàng)參數(shù)，用于保持0水平集的拓?fù)浣Y(jié)構(gòu);ε為域值，讓水平線函數(shù)保持0-1的形態(tài)，從而使函數(shù)變?yōu)榉侄纬Ｖ岛瘮?shù)，便于計(jì)算和保持良好的形態(tài);τ為以及k迭代過程中的時(shí)間步長;τ2為s迭代過程中的時(shí)間步長;τ3為迭代過程中的時(shí)間步長;以上三項(xiàng)均根據(jù)時(shí)間步長基本原則進(jìn)行取值，為長期實(shí)驗(yàn)所取得的經(jīng)驗(yàn)值;l為拓?fù)浔３志嚯x，l越大，拓?fù)浔３志嚯x越遠(yuǎn);d為高斯卷積核的標(biāo)準(zhǔn)差， window為搜索框大小;圖像大小300×354像素。
對比圖1圖2兩組圖像，可以看出兩組運(yùn)行結(jié)果基本一致，但本文改進(jìn)的模型及方法相比原來的模型及AOS方法迭代次數(shù)有了大幅度地減少。本文改進(jìn)的模型及方法在迭代到第100步時(shí)便達(dá)到收斂條件，完成分割，并且拓?fù)浣Y(jié)構(gòu)得到了保持，而傳統(tǒng)的模型及方法900步才完成分割。
圖3為SRS模型的AOS方法分割人為噪聲圖像的效果圖，圖4為ISRS模型的對偶算法分割人為噪聲圖像的效果圖。圖3參數(shù)分別為1、0.5、1、2×2、-0.2、-0.2、1，圖像大小200×138像素。圖4參數(shù)γ、α、β、ε、τ、τ2、τ3、l、d、window分別為8、10、0.01、0.1、0.5、0.2、0.125、1、4、3×3，圖像大小400×275像素。
圖3 SRS模型的AOS方法分割人為噪聲圖像
（a）初始輪廓線;（b）迭代200次效果;（c）迭代400次效果;（d）迭代800次效果
圖4 ISRS模型的對偶算法分割人為噪聲圖像
（a）初始輪廓線;（b）迭代10次效果;（c）迭代20次效果;（d）迭代40次效果;（e）迭代80次效果;（f）迭代200次效果
對比圖3圖4，兩種方法的運(yùn)行結(jié)果都較好地保持了圖像本身的拓?fù)浣Y(jié)構(gòu)，但本文改進(jìn)的模型及方法在迭代到第80步時(shí)便達(dá)到收斂條件，完成分割，并且拓?fù)浣Y(jié)構(gòu)得到了保持，而傳統(tǒng)模型及AOS方法在第800步才完成分割。
圖5為SRS模型的AOS方法分割腦部切片1的效果圖，圖6為ISRS模型的對偶算法分割腦部切片1的效果圖。圖5參數(shù)h、τ、l、window、α、β、ε分別為1、0.5、1、2×2、-0.2、-0.2、1，圖像大小125×125像素。圖6參數(shù)γ、α、β、ε、τ、τ2、τ3、l、d、window分別為1、50、0.01、1、0.3、0.12、0.125、1、4、5×5，圖像大小300×300像素。
圖5 SRS模型的AOS方法分割腦部切片1
（a）初始輪廓線;（b）迭代200次效果;（c）迭代400次效果;（d）迭代700次效果
圖6 ISRS模型的對偶算法分割腦部切片1
（a）初始輪廓線;（b）迭代10次效果;（c）迭代20次效果;（d）迭代30次效果;（e）迭代50次效果;（f）迭代100次效果
對比圖5圖6，本文改進(jìn)的模型及方法同原模型的AOS方法相比能更好地進(jìn)入到凹槽內(nèi)，并且在迭代次數(shù)上也有了大幅度地減少。本文改進(jìn)的模型及方法在迭代到第50步時(shí)便達(dá)到收斂條件，完成分割，并且拓?fù)浣Y(jié)構(gòu)得到了保持，而傳統(tǒng)模型及AOS方法在第700步才完成分割。
圖7為SRS模型的AOS方法分割腦部切片2的效果圖，圖8為ISRS模型的對偶算法分割腦部切片2的效果圖。圖7參數(shù)h、τ、l、window、α、β、ε分別為1、0.5、1、2×2、-0.2、-0.2、1，圖像大小125×125像素。圖8參數(shù)γ、α、β、ε、τ、τ2、τ3、l、d、window分別為1、50、0.01、10、0.3、0.12、0.125、1、4、5×5，圖像大小300×300像素。
（a）初始輪廓線;（b）迭代200次效果;（c）迭代400次效果;（d）迭代800次效果
（a）初始輪廓線;（b）迭代10次效果;（c）迭代20次效果;（d）迭代30次效果;（e）迭代50次效果;（f）迭代100次效果
對比圖7圖8，兩種方法都保持了對象本身的拓?fù)浣Y(jié)構(gòu)，本文改進(jìn)的模型及方法在迭代到第50步時(shí)便達(dá)到收斂條件，完成分割，并且拓?fù)浣Y(jié)構(gòu)得到了保持，而傳統(tǒng)模型及AOS方法在第800步才完成分割。
表1給出了以上四個(gè)對比實(shí)驗(yàn)的迭代次數(shù)、運(yùn)算時(shí)間以及運(yùn)算時(shí)間加速比的具體數(shù)值?？芍谟袇?shù)影響的情況下，本文也能較好地保持分割對象的拓?fù)浣Y(jié)構(gòu)，并在保持對象拓?fù)浣Y(jié)構(gòu)的基礎(chǔ)上大幅度減少了算法的迭代次數(shù)，模型的收斂速度也得到了大幅度提高。
4 結(jié)論
為了進(jìn)一步提高拓?fù)浔３址指钅Ｐ偷挠?jì)算效率，本文提出了ISRS模型及其對偶算法，簡化了計(jì)算過程，提高了運(yùn)算速度;GVF有向力場明顯加快了輪廓線在圖像狹窄區(qū)域的演化速度。對比實(shí)驗(yàn)表明，本文模型及其對偶算法的解類似于經(jīng)典SRS模型及其方法的解，可以有效防止分割輪廓的合并或分裂，保持原來的拓?fù)浣Y(jié)構(gòu)，并且在計(jì)算效率方面較經(jīng)典模型與算法有較大提高。下一步工作可以嘗試用高階模型來忽略分割過程中的微小細(xì)節(jié)，保持邊緣的紋理特征。
參考文獻(xiàn)
[1]國凱，王國棟，潘振寬，等.CT影像中的主動(dòng)脈血管鈣化檢測[J].青島大學(xué)學(xué)報(bào)（自然科學(xué)版），2013，26（1）：50-54.
[2]田佑仕，黃寶香，姜濤，等.融合圖像恢復(fù)的遙感數(shù)據(jù)變分分割方法研究[J].青島大學(xué)學(xué)報(bào)（自然科學(xué)版），2020，33（4）：1-8.
[3]吳則舉，王嘉琦，焦翠娟，等.基于深度學(xué)習(xí)網(wǎng)絡(luò)U-Net的輪胎帶束層分割算法研究[J].青島大學(xué)學(xué)報(bào)（自然科學(xué)版），2019，32（4）：22-29+35.
[4]MACDONALD D， KABANI N， AVIS D， et al. Automated 3-D extraction of inner and outer surfaces of cerebral cortex from MRI[J]. NeuroImage， 2000， 12 （3）：340-356.
[5]ZENG Y， SAMARAS D， CHEN W， et al. Topology cuts： A novel min-cut/max-flow algorithm for topology preserving segmentation in N-D images[J]. Computer Vision and Image Understanding， 2008， 112 （1）：81-90.
[6]CHEN M， CARASS A， OH J， et al. Automatic magnetic resonance spinal cord segmentation with topology constraints for variable fields of view[J]. NeuroImage， 2013，83： 1051-1062.
[7]WAGGONER J， ZHOU Y J， SIMMONS J， et al. Topology-preserving multi-label image segmentation[C]// IEEE Winter Conference on Applications of Computer Vision （WACV 2015）， Waikoloa， 2015： 1084-1091.
[8]LIU G Q， LI H F， YANG L. A topology preserving method of evolving contours based on sparsity constraint for object segmentation[J]. IEEE Access， 2017， 5：19971-19982.
[9]DEBROUX N， OZERE， S， LE GUYADER C. A non-local topology-preserving segmentation-guided registration model[J]. Journal of Mathematical Imaging and Vision， 2017， 59 （3）：432-455.
[10] DEBROUX N， LE GUYADER C. A joint segmentation/registration model based on a nonlocal characterization of weighted total variation and nonlocal shape descriptors[J]. SIAM Journal on Imaging Sciences， 2018， 11 （2）：957-990.
[11] CHAN H L， YAN S， LUI L M， et al. Topology-preserving image segmentation by beltrami representation of shapes[J]. Journal of Mathematical Imaging and Vision， 2018， 60 （3）：401-421.
[12] ZHANG D P， LUI L M. Topology-Preserving 3D image segmentation based on hyperelastic regularization[J]. Journal of Scientific Computing， 2021， 87（3）： 74.
[13] CASELLES V， KIMMEL R， SAPIRO G. Geodesic active contour. IJCV[J]. International Journal of Computer Vision， 1997， 22（1）：61-79.
[14] ZHAO H， CHAN T， MERRIMAN B， et al. A variational level set approach to multiphase motion[J]. Journal of Computational Physics， 1996， 127（1）： 179-195.
[15] 郭振波.基于變分水平集方法的多相圖像分割研究[D]. 青島： 中國海洋大學(xué)， 2008.
[16] HAN X， XU C Y， PRINCE J L. A topology preserving level set method for geometric deformable models[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence， 2003， 25（6）：755-768.
[17] ALEXANDROV O， SANTOSA F. A topology-preserving level set method for shape optimization[J]. Journal of Computational Physics， 2005， 204（1）：121-130.
[18] SUNDARAMOORTHI G， YEZZI A. Global regularizing flows with topology preservation for active contours and polygons[J]. IEEE Transactions on Image Processing， 2007， 16（3）：803-812.
[19] ROCHERY M， JERMYN I H， ZERUBIA J. Higher order active contours[J]. International Journal of Computer Vision， 2006， 69（1）：27-42.
[20] GUYADER C L， VESE L A. Self-repelling snakes for topology-preserving segmentation models[J]. IEEE Transactions on Image Processing， 2008， 17（5）：767-779.
[21] PAN H， SONG J， LIU W， et al. Using the split bregman algorithm to solve the self-repelling snake model＼[J＼]. Computer Vision and Pattern Recognition， 2021＼[2021-08-15＼]https：//doi.org/10.1007/s10851-021-01065-9.
[22] XU C， PRINCE J L. Snakes， shapes， and gradient vector flow[J]. IEEE Transactions on Image Processing， 1998：7（3）：359-369.
[23] CHAMBOLLE A. An algorithm for total variation minimization and applications[J]. Journal of Mathematical Imaging and Vision， 2004， 20（1-2）： 89-97.
[24] OSHER S， SETHIAN J A. Fronts propagating with curvature-dependent speed： Algorithms based on Hamilton-Jacobi formulations＼[J＼]. Journal of Computational Physics， 1988， 79 （1）： 12-49.
Topology Preserving Segmentation Model Based on
GVF Force and Its Dual Algorithm
SHEN Meng-jie， PAN Zhen-kuan， SONG Jin-tao， WEI Wei-bo
（College of Computer Science & Technology， Qingdao University， Qingdao 266071， China）
Abstract：
In view of the insufficient force of the self-repelling Snake model on narrow image regions， the calculation of traditional additive operator splitting method is complicated， and the memory usage will also increase rapidly with the increase of the image size. It is proposed that the Gradient Vector Flow directed force field based on the original model was added to accelerate the evolution speed of the contour line in the narrow area of the image， and a fast Dual algorithm is designed for the improved model to simplify the algorithm design and improve the efficiency of the solution. Numerous numerical experiments show that the proposed improved model and algorithm have a greater improvement in computational efficiency than the classic model and algorithm.
Keywords：
self-repelling snake model; topology-preserving segmentation; dual algorithm; gradient vector flow; variational method
收稿日期：2021-09-10
基金項(xiàng)目：
國家自然科學(xué)基金（批準(zhǔn)號：61772294，11472144） 資助;山東省聯(lián)合基金（批準(zhǔn)號：ZR2019LZH002） 資助。
通信作者：
潘振寬，男，博士，教授，主要研究方向?yàn)橛?jì)算機(jī)視覺、圖像處理，E-mail： zkpan@126.com
