基于高階正則與非光滑數據擬合項的圖像邊緣檢測模型①

李 春,陳靜思,王鵬彥,李 健,羅 澤

1(中國科學院 計算機網絡信息中心,北京 100190)

2(中國科學院大學,北京 100049)

3(云南財經大學 云南省經濟社會大數據研究院,昆明 650221)

4(四川臥龍國家級自然保護區管理局,臥龍 623006)

1 圖像分割簡介

圖像分割和邊緣檢測是計算機視覺,特別是圖像分析中一個基礎而重要的問題(如:目標識別和圖像解析),其目的是把一幅給定圖像根據其圖像特征(如:邊緣,顏色,紋理,運動,角點特征等)劃分為不同區域[1].近年來,隨著計算機技術發展,使得圖像分割成文現代科學領域的研究熱點,圖像分割在醫學圖像分析,遙感,監測等領域得到廣泛應用[2–11].在諸如物體識別,分類和分割的應用中尋找物體邊緣是非常重要的一步.因此,使用和設計何種邊緣檢測算法直接影響這些應用性能.如:2019年Yuan等[12]提出了一個量子圖像邊緣檢測算法,將邊緣檢測算法靈活應用于量子計算.圖像邊緣檢測在雷達圖像檢測中也有重要應用,Ma等[13]提出通過稀疏表示進行SAR圖像邊緣檢測,作者提出了一種基于去噪算法的新型合成孔徑雷達(SAR)圖像檢測算法,該算法通過稀疏表達和一種新的形態學邊緣檢測器.首先,作者將Shearlet變換應用于SAR圖像以獲得圖像的稀疏表示.然后,將帶方向的形態學邊緣檢測器應用于Shearlet的方向子帶系數,其通過迭代去噪處理來恢復.2019年Sert等[14]提出了一個新的基于最大范數熵的中智學邊緣檢測方法.同時,邊緣檢測算法在醫學圖像領域也有重要應用,隨著醫學領域成像技術的快速發展,不同模態的醫學圖像具有不同的成像原理,反映了人體生理信息的不同重點和缺陷.邊緣檢測是醫學超聲圖像處理的關鍵步驟,檢測結果將直接影響醫生對疾病的診斷.圖像邊緣檢測可以看作邊緣點和非邊緣點的分類問題.基于此,Song等于2019年[15]提出了基于改進差分進化算法和Prewitt算子的醫學圖像邊緣檢測.Churchill等于2019年[16]提供了理論證據來解釋邊緣位置在CT影像重建中的潛在重要性.然后開發一種實用的方案來利用這一理論.最后,作者通過實驗證實了提出的使用邊緣遮蔽正則化可以提高CT重建的準確性和重建速度.

近年來,基于變分法和偏微分方程(PDES)的圖像分割方法得到了廣泛研究[17–24].早期圖像分割模型代表為活動輪廓模型,活動輪廓模型根據圖像邊緣信息和區域信息對給定圖像進行分割.基于邊緣信息和區域信息考慮,Kass等于1988年提出基于邊緣信息的能量最小化方法,即,Snake模型[25].其后,受到Snake活動輪廓模型啟發,Caselles等提出了綜合活動輪廓長度信息,提出了幾何活動輪廓模型[26],與此同時,Mumford和Shah提出了基于區域信息圖像分割變分框架,即,著名的Munmford-Shah (MS)模型[27].但是,由于MS模型能量泛函的非凸性,導致在當時對該模型直接求解十分困難.所以,Chan和Vese為了克服MS模型所存在的弊端,并整合了變分法和水平集方法[17],提出了Chan-Vese (CV)模型[28].

由于CV模型應用廣泛,所以在過去幾十年里,該模型一直受到研究學者青睞.同時,各研究學者開發各種對該模型求解的有效方法,例如:為克服初始化輪廓敏感度問題,Li等[29]修改CV模型,提出了利用懲罰函數作為約束且無需初始化過程的圖像分割模型.2014年Duan等[30]綜合了變量分離法,對偶方法,Bregman迭代和增廣的拉格朗日方法提出了更快更有效更的算法.這些算法在高計算精度和無需初始化方面取得重要成果.

隨著計算機技術的發展,深度學習在圖像分割中取得重要進展,最顯著的進展是深度卷積神經網絡(CNNS)不斷被應用到場景解析中.典型的網絡包括:DeepLab (V1-V3)[31,32],Refinenet[33],PSP-Net[34]和DANet[35],然而,這些模型的缺點是計算代價相當大,需要學習的參數特別多,為了降低這些模型在場景解析中的計算復雜度,過去幾年里,很多學者在降低計算復雜度方面做了大量工作.例如:ENet[36],ESPNet[37]大量減少了需要學習的參數.

但是,上述所提的變分優化模型中,基本上都是假設圖像所含噪聲為高斯噪聲,幾乎是用L2范數對高斯噪聲加以擬合.然而在現實生活中,情況相對較為復雜,往往圖像在成像時還有來自外加噪聲的污染,幸運的是,外加噪聲可以用L1范數得以很好擬合[38,39].利用L1范數作為數據擬合項作為分段常數圖像分割模型,并利用二階正則函數(TV2)對目標函數加以懲罰,使得所提出模型能更好地逼近目標函數的同時還能能更好地分割低對比度和含有外加噪聲的圖像.再者,本文通過變量分離方法把一復雜非凸問題轉換為若干個簡單凸子問題進行求解,從而巧妙地處理了不可導項和正則項.

2 預備知識

在本節中,因為Alternating Direction Method of Multipliers (ADMM)或Split Brgman迭代[40]方法是求解線性約束凸二次規劃問題的有效方法,所以在此對該方法做如下簡單介紹.考慮如下線性約束凸優化問題:

其中,通過上述目標函數,可得如下增廣拉格朗日泛函:

當上述系統方程滿足某些收斂條件時,ADMM方法在(2)中固定v更新u,固定u更 新v.因此,ADMM算法可以總結如算法1.

算法1.ADMM算法γ>0u0v0 b0=0初始化:設k=0,選擇 ,,and ;

for do(1)解問題argmin k=1,2,3,···,K u E1(u)+γ 2A1u+A2v?f?bk2 2;(2)求解問題argmin u E2(v)+γ 2A1u+A2v?f?bk2 2;(3)拉格朗日乘子更新bk+1=bk?(A1uk+1+A2vk+1?f);end

3 模型提出與數值求解

在本節中將介紹所提出的高階正則圖像輪廓檢測算法,隨后用變量分離法和ADMM算法對該模型詳細優化過程和部分收斂性分析加以描述.

設f:[0,1]2→R?L1(?)為觀測到信號(圖像),且滿足如下數學表達式,f=u+n,其中u為未受噪聲污染信號(圖像),n為外加噪聲(IN).為了重建u,可以考慮如下優化問題:

其中,?u為u的 分布導數,∥ ?u∥1為著名的total variation(TV)[41]有界變差空間中的半范數,λ>0為權參數.若用L2-數據擬合項 ∥f?u∥2替代∥f?u∥1數據擬合項,從而可得經典的Rudin-Osher-Fatemi (ROF)[41]模型.從而圖像分割模型數學表達式如下:

其中,Γ 為u的 不連續集合,H1為一維的Hausdorff測度,即,Γ的長度.

假設區域 [0,1]2被劃分為如下的N個區域,?1,?2,···,?N且同時滿足∪i=1,···,N?i=[0,1]2,?i∩?j=?,ij,且對于所有的i具 有ci∈R.特別地,當N=2可以得到除了數據擬合項之外,其他都類似于CV模型的模型.

由于 ∥f?u∥1不 可微,所以可以通過迭代得出c1c2近似解.在圖像恢復中,盡管基于TV[41]正則已經取,得了很好的效果,但是在圖像恢復過程中TV正則會引起階梯效應.而高階正則不僅有可以去除階梯效應的功能,而且高階正則還有光滑區域,很好地保持物體邊緣的效果,在本模型的實驗過程中,假設圖像都是帶有噪聲的,又因為高階正則比低階正則更具有光滑性和非線性性,所以為了更高的提取邊緣信息,綜合Jung等[38]提出的模型,我們提出如下的高階正則圖像邊緣檢測模型:

通過觀察上述目標函數可以發現,為了對該目標泛函進行求解,(c1,c2,u),即,兩步驟算法:第1步:先固定u求解c1,c2,第2步:固定c1,c2求解u.同時并對離散的散度算子做了如下定義,且散度算子div:(RM×N)2→RM×N具有如下的共軛性質.?div.u=p.?u,?u∈RN×M,p∈(RM×N)2.

故此,離散散度算子可以作如下定義:對

其中,Dx和Dy為如下的向前向后分算子.

以此類推,離散二階散度算子定義如下:div2:(RN×M)4→RN×M.且具有如下共軛性質:div2q.u=q.?u,?u∈RM×N,q∈(RN×M)4.對于,q=(q11,q12,q21,q22)∈(RN×M)4,我們定義,div2q(i,j)=Dxxq11(i,j)+Dxyq12(i,j)+Dyyq21(i,j)+Dyxq22(i,j).

在實際應用中,學者經常用算子分離法(split Bregman iteration和ADMM)對(6)進行求解.為達到求解目的,本文先把(6)無約束優化問題轉化成一個有約束優化問題,從而利用Bregman iteration進行求解,引入中間變量的好處在于:通過引入如下輔助變量,不僅把無約束優化問題變成一個有約束優化問題,引入輔助變量使用變量分離法把一個復雜問題分解成若干簡單子問題,使得分別對各子問題求解相對容易.

其中,v=(vx,vy),w=(wxx,wyy,wxy,wyx).這里需要進行說明的是該符號不是對v,w分別求導,而是表示坐標.顯而易見,無約束優化問題(6)和有約束優化問題(7)等價,從而可以利用ADMM或Split Bregman迭代進行求解.其中,問題(7)等價于:

Bregmen迭代算法總結如算法2.

算法2.Bregman迭代圖像邊緣檢測算法u1=u0,v1=?u0,w1=?u0 bk1=0,bk2=0初始化:設k=0,選擇 and ;for do 計算(8)的最優解;bk+1 k=1,2,3,···,K 1 =bk1+vk+1??uk+1 更新:;2 =bk2+wk+1??uk+1 更新:;bk+1 end

直接對(8)求解時十分困難,從而可把上式分離成多個單變量,即,把一個復雜問題分解成幾個子問題求解.Split Bregman iteration第一次被Goldstain和OSher[40]提出,它在很多寬松條件下等價于ADMM算法.該方法被廣泛地應用到圖像噪聲去除,圖像去模糊,圖像修補等相關領域.

3.1 u-子問題

對u-子問題進行求解時首先固定r=|f?c1|+|f?c2|.得到如下目標泛函:

此優化問題可以通過其優化條件加以求解,在連續情況下,其解為四階線性PDE,離散情況,其優化條件如下:

通過FFT (快速傅里葉變換)可以得到(10)如下近似解.其中,F 代表傅里葉變換,F?1代表傅里葉的逆變換.

3.2 v-子問題

v-子問題等價于求解如下優化問題:

上述優化問題的解可以使用收縮算子[42]表達,其解析表達如下:

其中,收縮算子shrink(x)[42]定義如下:

對于任意點α ∈[0,1]2,為了方便起見,令

3.3 w-子問題

w-子問題滿足如下優化條件:

類似于v-子問題,上述優化問題解為具有收縮性質的軟閾算子.

為了說明收縮算子的表達效果和設置依據我們考慮如下的一維最小值問題加以說明:

其中,x,h都為標量.令可得:

當h≥λα時 ,Q(x)的 最小值為h?λα≥0;當h≤?λα時,Q(x)的最小值為h+λα≤0 ;當? λα

收縮算子對稀疏解是有利的,最終大量的h值會趨近于區間[ ?λα,λα]的 某個值從而使得x變 得稀疏.當λ 值越大,區間[ ?λα,λα]就 越大,從而使得x的解就越稀疏.

對于本模型來說若圖像的大小為N×N,以式(13)為例,其線性時間復雜度為O (N2),所以,v,w-子問題可以被快速求解,綜上所述,對于v,w-子問題我們設置了如式(14)的收縮算子是合理的.

3.4 c 1,c2-子問題

現考慮如下最小值問題,從而得到c1,c2的近似解.先固定u,得到如下最小值問題:

其中,當i=1時 ,h=u,當i=2 時 ,h=1?u.

首先我們定義區域 [0,1]2上的特征函數,χ=從而常數值ci可以用如下的表達式來代替χci,從而可以重寫最小值(20):

其中,

為使用變量分離法,故此引入如下輔助變量,

式(23)可以通過ADMM或Split Bregman迭代進行求解.

從而,cki+1和eki+1的近似解解析表達式為:

故此ADMM圖像邊緣檢測算法總結如算法3.

算法3.ADMM圖像邊緣檢測算法λη1η2α β>0γ1,γ2>0,bk1=v0=0,bk2=w0=0,e0i=d0i=0 u0=1 ??[0,1]2初始化:設k=0,選擇,,,,,在某些區域 ,0其他;k=1,2,3,···,K for do hk=uk i=1hk=(1?uk),i=2 if ,,ck+1 i ck+1 i i=1,2 計算 :用式(27)計算 ,;ek+1 i ek+1 i i=1,2 計算 :用式(28)計算 ,;dk+1 i dk+1 i i=1,2f?ck+1 計算 :用 ;更新 :用式(26)計算 ,;rk+1 rk+1=f?ck+1 1 2 vk+1 vk+1 計算 :用式(13)計算 ;wk+1 wk+1 計算 :用式(16)計算 ;uk+1 uk+1 計算 :用式(11)計算 ;1 bk+11 =bk1+vk+1??uk+1 更新 :;bk+1 2 bk+12 =bk2+wk+1??uk+1 更新 :;bk+1 end

4 算法部分收斂性分析

本節中,受到文獻[38]啟發,本文將給出本算法的部分收斂性分析,最小值問題(6)可以通過引入變量分離辦法重寫成有約束條件:

其對應的拉格朗日泛函如下:

其中,μ1,μ2,μ3為對偶變量.

設X?=(c?1,c?2,u?,e?1,e?2,v?,w?,μ?1,μ?2,μ?3,μ?4)為 問 題(29)的KKT點,則點X?滿足如下的KKT條件:

其中,x⊙y為coordi nate-wise乘法.

定理.設Xk=(ck1,ck2,uk,ek1,ek2,vk,wk,d1k,d2k,bk1,bk2)為算法3中的迭代,并且,設:

假設,

證明:首先,對于u-子問題等價于求解最小化問題(9),從而得出如下優化條件:

其次,從優化條件和ADMM算法對各變量的表達式可得:

5 實驗結果

為了說明本文提出算法優越性,本文將給出該模型在灰度圖像,真實圖像,含噪聲圖像,CT圖像的分割實驗結果,在討論結果之前,本文先作以下說明,該算法的迭代終止條件應滿足如下的相對容忍度,即滿足下列條件:

此后,用該模型和Chan-Vese (CV)做比較.該算法與CV模型不同之處在于計算常數值ci,i=1,2,對于CV模型而言,ci,i=1,2在每次迭代過程中都有具體迭代表達式.該模型另外一個優點是很少依賴于初始化位置和微調參數,在本算法中,η1,η2,α ,β >0 ,γ1,γ2>0,所有參數都相對很固定.因為該模型是一個局部極小值問題,不同初始輪廓可能會收斂到不到能量局部極小值.我們之所以選擇和CV模型作比較,主要原因在于想體現高階正則模型與低階正則模型對圖像邊緣檢測不同效果影響.

圖1為CV模型和新模型對飛機圖像的分割效果圖,其中,第一,二行分別代表CV模型和高階正則模型對飛機圖像的分割效果圖;第一列代表原始圖像,第二列分別代表CV模型和高階正則模型對飛機圖像的初始化輪廓;第三列分別代表CV模型對飛機圖像迭代800次時的邊緣檢測情況和高階正則模型對飛機圖像迭代80次時的邊緣檢測情況;第四列分別代表CV模型對飛機圖像迭代800次時的實際分割效果圖和高階正則模型對飛機圖像迭代80次時的實際分割效果圖.

圖1 CV模型和高階正則模型對飛機圖像的檢測效果圖

圖2表示CV模型對自然風景圖像的分割效果圖,圖3表示高階正則模型對自然風景圖像的分割效果圖.其中圖2的4幅圖分別為:第一幅代表原始圖像,第二幅代表CV模型對自然風景圖像的初始化輪廓,第三幅代表CV模型對自然風景這幅圖像迭代800次時的邊緣檢測圖,第四幅代表CV模型對自然風景這幅圖迭代800次時的實際分割效果圖.圖3的4幅圖分別代表的意義為:第一幅為原始圖像,第二幅為高階正則模型對自然風景圖像的初始化輪廓圖,第三幅圖像表示高階正則模型對自然風景圖迭代200次時的邊緣檢測圖,第三幅圖像表示高階正則模型對自然風景圖迭代200次時的實際分割效果圖.

圖2 CV模型對自然風景圖像的分割效果圖

從圖4可以看出,CV模型對自然風景這一幅圖像迭代30次后基本收斂,而高階正則模型需要迭代到40次左右才收斂,但是值得注意的是,從兩個模型對同一幅圖像Loss函數曲線圖可知,CV模型基本收斂后Loss函數曲線圖還有小幅度跳躍情況,然而,高階正則模型收斂后Loss函數曲線圖相對更平滑,說明其分割效果也很好,更能精確地逼近物體邊界.

下面我們給出一個例子說明該高階正則模型在醫學圖像領域的應用,圖像分割在醫學圖像領域的應用就是通過發展大量的自動或者半自動圖像分割方法,精準分割醫學圖像,從而替代醫生的大量手工標注,從而使得醫生從大量繁重的體力勞動中解放出來.

我們用本模型在病人肺部CT圖上做分割實驗,如圖5所示,從其分割效果圖得知,高模型不僅能對病人肺部做精準分割,CV模型和高階正則模型Loss函數曲線收斂圖如圖6,可知,新模型Loss函數曲線圖相對CV模型而言相對較光滑,這更能說明該模型比CV模型來更能精確地逼近物體邊界.

需要說明的是,在分割過程中,CV模型的模型參數設置均為原始論文的參數,考慮到CV模型的數據項數擬合高斯噪聲,而本文所設計的模型數據保真項模擬的是外加噪聲,所以,在實驗過程中,對于CV模型而言,并沒有人工地加入任何噪聲.但是,對于高階正則模型而言,由于本文假設其成像過程中可能含有外加噪聲(如:椒鹽噪聲),所以在實驗過程中我們用Matlab中的加噪聲函數imnoise(Im0,'salt ＆ pepper' 0.6);每幅原始圖像都人為的加了參數為0.6的椒鹽噪聲.新模型對所有圖像分割的迭代總次數均設置為200,CV迭代總次數均設為800,迭代停止條件為stoppling_tol=–10,λ =20,α =β=1,η1=η2=5 0,γ1=γ1=50.

圖3 高階正則模型對自然風景圖像的分割效果圖

圖4 CV模型和高階正則模型對自然風景圖像分割的Loss函數曲線收斂圖

圖5 高階正則模型對病人肺部CT圖像的分割效果圖

圖6 CV模型和高階正則模型對病人肺部CT圖像分割的Loss函數曲線收斂圖

6 總結與展望

本文提出一個修改的Chan-Vese模型,并引入高階正則函數對目標函數進行懲罰,然后對新模型用ADMM或者Split Bregmen迭代算法進行數值求解.再者,本文還分析了該模型的數學性質,并給出該模型的部分收斂性分析,實驗結果表明,通過引入高階正則函數后,該模型不僅能分割對比對低的物體,而且還可以探測帶噪聲的物體.大量實驗表明,該算法在各領域具有廣泛應用價值,在未來的工作中,我們會綜合一些深度學習的模塊,去探測物體的landmarks點,從而使得算法更自動,更魯棒.